This section outlines the main methods used in this investigation. A test case of particulate organic carbon (POC) in the North Atlantic down to 1000 m is used. The observed data, explained in detail in Sect. 2.1, are obtained from autonomous ocean Argo floats. The model tested is the one-dimensional (depth) configuration of the ocean biogeochemical model PISCES-v2_RC (Aumont et al., 2015, 2017), a component of NEMO4 v4.0.1 (Nucleus for European Modelling of the Ocean version 4), as outlined in Sect. 2.2.

The type of GA used is BRKGA (Goncalves and Resende, 2011). The outline of this method, including the crossover, is described in Sect. 2.3. We use the workflow manager Autosubmit (Manubens-Gil et al., 2016; Uruchi et al., 2021) to create a workflow that facilitates the various steps of the algorithm, as outlined in Sect. 2.4.

This paper outlines two test case experiments in which the reference data are an output of a simulation with default parameters, another three in which the reference data are observed data from three locations in the North Atlantic, and, last, a set of cross-experiments. Section 2.5 outlines the details of these experiments.

2.2 PISCES 1D and parameters

PISCES-v2 (Aumont et al., 2015) is an OBGCM of intermediate complexity that represents the cycles of the main inorganic nutrients (N, P, Si, and Fe), carbonate chemistry, and organic matter compartments, including phytoplankton and zooplankton organisms (with two size classes each), dissolved organic matter, and particulate organic matter, making up 24 prognostic variables or tracers in total. Here we use a model version, PISCES-v2_RC, that incorporates the POC reactivity continuum parameterisation (Aumont et al., 2017). This model version is included as the OBGCM component of NEMO4 v4.0.1 (Madec et al., 2022) and will hereafter be referred to as PISCES.

In PISCES, detrital POC is represented by two tracers, i.e. POC for detritus smaller than 100 µm and GOC for detritus larger than 100 µm. To avoid confusion between PISCES tracers and the term POC, used here as a generic concept and to refer to observations, PISCES tracer names are italicised. It is important to note that total POC as sampled in situ is made up of detrital matter and living biomass. Therefore, the correspondence between PISCES tracers and observations must be established. Here we define SPOC as the sum of the PISCES tracers for nanophytoplankton (PHY), microphytoplankton (PHY2), microzooplankton (ZOO) and small detritus (POC), and LPOC as the sum of large detritus (GOC) and mesozooplankton (ZOO2). These idealised fractions show good correspondence with those determined from BGC–Argo data (Galí et al., 2022). It is important to keep in mind that detrital POC is a variable proportion of total POC, which generally increases with depth. In the mesopelagic, detrital POC represents around 70 % of total POC globally with the default PISCES parameterisation (Table 3 in Galí et al., 2022).

Our study focuses on nine PISCES parameters (Table 1) expected to strongly influence mesopelagic POC dynamics according to model equations (Aumont et al., 2015, 2017) and preliminary analyses (Appendix A and B). These parameters control POC formation in the surface productive layer through microphytoplankton mortality, gravitational POC fluxes, POC degradation rates, and interception and fragmentation of sinking POC by mesopelagic zooplankton. Preliminary tests also included the parameters unass and unass2, which determine POC production from the unassimilated fraction of phytoplankton biomass ingested by zooplankton. However, they were eventually excluded because these parameters have a strong impact on upper ocean (epipelagic) ecosystem dynamics, which are beyond the scope of our study.

Table 1Definitions of the PISCES parameters included in the optimisation experiments, along with their default values, optimisation ranges, and units.

Download Print Version | Download XLSX

This investigation uses PISCES configured for one spatial dimension (1D) and to run offline (Galí et al., 2022). The 1D configuration has the same vertical levels as the 3D configuration (in our setup, 75 levels of gradually increasing thickness – L75 vertical grid), but the horizontal grid is reduced to an idealised domain of 3×3 cells. In this configuration, tracer concentrations change over the temporal and vertical dimensions as a result of local sources and sinks, vertical diffusion, particle sinking through the water column, and fluxes at the ocean–atmosphere boundary. PISCES computes the sources and sinks and the gravitational sinking of detrital particles at each biological time step (here set to 45 min, which is one-fourth of the NEMO4 v4.0.1 time step). Then, the NEMO component TOP (Tracers in the Ocean Paradigm) calculates vertical diffusion using dynamical fields, which are precalculated in a previous NEMO run, with a time step of 3 h. The 1D configuration does not allow for the advection of biogeochemical tracers. Simulations are spun up by repeating the same annual forcing over 4 years, and simulation year 5 is used for the comparison against observations.

Being one-dimensional, the model only requires one computational core and runs at a speed of roughly 1 simulation year per minute on a supercomputer, which allows for multiple simulations to be run in parallel. The numerical parameters that will be constrained are stored in text files called name lists and can be easily modified prior to each simulation without requiring recompilation. In the experiments (Sect. 2.5), parameters were allowed to vary between the lower and upper bounds based on what we considered physically or biologically reasonable according to the experimental and modelling literature.

2.3 Genetic algorithm (GA)

A GA is a type of evolutionary algorithm used for optimisation that, in general, is analogous to natural selection in the sense that a population of p individuals are tested for their strength (or fitness) using a cost function. At each generation, weaker individuals are eliminated, while stronger individuals pass on their characteristics by pairing with other individuals to produce λ offspring. In most applications, including this one, p=λ. A GA is considered a stochastic optimisation method, which is well balanced between elitist and exploratory behaviours. Being elitist in this sense is the property of reaching an optimal solution with efficiency, and being exploratory refers to increasing the range of possible solutions. Being exploratory is particularly important to ensure that the algorithm does not reach a local minimum of the cost function by leaving some regions of the search space unexplored. The usual method of recombination in the GA is the crossover, which is the action of two individuals from a generation producing offspring for the next. This is the primary discovery force of the GA. In our case, an individual is a vector of floating point numbers that represent the values of the parameters. A crossover occurs when two individuals are selected, and a new individual vector is created by taking a random combination of components from the two parent individuals. In general, crossovers are intended to be elitist by ensuring that individuals with higher strength are more likely to be chosen. This process is known as selective pressure.

Another feature inspired by genetics is the concept of mutations. The purpose of mutations is to make the algorithm more exploratory by randomly changing or perturbing parts in individual members or adding randomly generated individuals to the population. This is usually done with a very small probability, emulating transcription errors that occur within natural gene passing.

Once the crossovers are completed and the new generation is made, their strength is again measured and the process is repeated. This continues until a certain condition is met. This can be whenever the value of the cost function of the strongest member reaches a certain value, or if no change is noted after a certain number of generations, or simply after a predetermined number of generations.

2.3.1 Biased random key genetic algorithm (BRKGA)

A BRKGA is a particular type of GA in which each gene is a vector of floats rather than a bitstring, which is typical of traditional GAs (De Jong et al., 1993). This is useful for addressing the issue of uneven distance between solutions, inherent to bitstrings, and appropriate for this problem because the set of parameters to be optimised can be treated as a vector. The behaviour of the BRKGA can be adjusted by changing the so-called metaparameters (Fig. 1) that are described below. Initially, p sets of parameters are generated at random using a uniform distribution with appropriate bounds (Sect. 2.2). At each generation, the p_e individuals with the best score, known as the elite subpopulation, are selected, where $p_{\mathrm{e}}<p/\mathrm{2}$. These are passed directly to the next generation. The remainder of the vectors are placed into the non-elite subpopulation. Next, a set of p_m randomly generated vectors is introduced into the population as mutants and passed directly onto the next generation in order to make the algorithm more exploratory, performing the same role as mutations in standard GAs. The set of vectors of the next generation is completed by generating $p-(p_{\mathrm{e}}+p_{\mathrm{m}})$ vectors by crossover. A crossover in this case is a method used to generate a new vector by selecting two parents at random, and then each element of the new vector is randomly picked from one of the two parents. In a normal random key GA, the parents are selected completely at random from the whole of the previous set of parameters, with a 0.5 probability of an element coming from either parent. However, in a BRKGA, one parent vector comes from the elite set and the other from the non-elite set. In addition, the probability of an element coming from an elite parent is determined by ρ, where ρ>0.5. This has shown, in previous investigations, to cause faster convergence to an optimal solution (Goncalves and Resende, 2011). Finally, to make the algorithm more exploratory, after the crossover is completed, all values are slightly perturbed to allow the exploration of values close to those of the elite vectors. It is worth noting that this slight perturbation may allow the parameters to evolve beyond their initial range. Given that the parameter ranges are also not well constrained, this allows the algorithm to explore the possibility of finding optimal values outside the given range; however, the feasibility of the values is at the discretion of the user.

Figure 1A visualisation of the BRKGA's process from one generation to another (Júnior et al., 2020).

2.3.2 Cost function

Deciding on an ideal cost function to measure the misfit between the results of each simulation and the observed data requires a number of considerations. In this case, the limitations of the model itself and the particular properties of the data need to be taken into account. An important model limitation is that there exist inherent physical biases and, in some cases, uncertainties in the conversion factor between the model variable and its observed counterpart. In addition, we wish to compare trends, in particular the seasonality of the data. For this, simply calculating the difference between observed data and simulated outputs, or bias, is not sufficient.

To ensure sensible fitting, in addition to bias, the correlation and the normalised standard deviation need to be considered. The root mean square error, RMSE, is a widely used parameter in this type of investigation; however, in certain cases, it has been found to reward reductions in model variability, for example, over the seasonal cycle (Jolliff et al., 2009). An alternative metric known as the ST score is used. This is defined as follows:

$$\begin{array}{}\text{(1)}& \mathrm{ST}=\sqrt{{\text{Bias}}_{\mathrm{m}}^{\mathrm{2}}+S_{\mathrm{3}}^{\mathrm{2}}},\end{array}$$

where Bias_m of an individual simulation is defined as its mean bias (over all data points) divided by the mean bias of the individual with the highest bias in the particular generation, that is,

$$\begin{array}{}\text{(2)}& {\text{Bias}}_{\mathrm{m}}={\displaystyle \frac{{\mathrm{Bias}}_i}{{\mathrm{Bias}}_{max}}},\end{array}$$

and S₃ is a function of normalised standard deviation, σ, and correlation, R. Jolliff et al. (2009) test this particular cost function using bio-optical data, generally characterised by lognormal or similarly right-skewed distributions that reflect the exponential growth and decay of plankton organisms. For this reason, a normal logarithmic scale is used, a choice that is supported by preliminary experiments where the BRKGA performance with linear vs. logspace statistics was evaluated. Jolliff et al. (2009) state various possible formulae. Since it is of high importance to correctly determine seasonality in this investigation and in this field in general, it is most sensible to choose a cost function that prevents situations in which normalised standard deviation and bias are rewarded at the expense of correlation. Considering the three described options, preliminary tests indicated that S₃ served this purpose most appropriately, as follows:

$$\begin{array}{}\text{(3)}& S_{\mathrm{3}}=\mathrm{1.0}-\left(e^{-}{\displaystyle \frac{(\mathit{\sigma }-\mathrm{1.0}{)}^{\mathrm{2}}}{\mathrm{0.18}}}\right)\left({\displaystyle \frac{(\mathrm{1}+R)}{\mathrm{2}}}\right).\end{array}$$

2.4 Workflow

Running a BRKGA requires performing a number of iterations until a termination condition is achieved. This does not represent a technical challenge if the fitness function can be calculated directly from the generation members. However, in some cases, such as the one presented in this work, an external model is responsible for calculating the result that will be the input to the cost function. As a consequence, the need for the parallel execution and management of many different and interdependent tasks requires using tools called workflow managers or metaschedulers, which are commonly used to run ensemble experiments with climate models. Here we use a state-of-the-art workflow manager called Autosubmit (Manubens-Gil et al., 2016). Autosubmit is developed with ESMs in mind, and is typically used to run complex simulations composed of multiple different tasks executed in one or multiple clusters via SSH (secure shell) connection. Autosubmit can automatically handle the submission of these tasks, respecting their dependencies and managing failures with minimal user intervention, thereby providing tools to monitor (Uruchi et al., 2021) the experiment execution. In addition, it allows multiple jobs to run simultaneously in parallel or packed in macrojobs (wrappers) by automatically allocating the required computing resources.

Autosubmit experiments are hierarchically composed of start dates, members, and chunks. A single experiment can run different start dates that can be divided into members in which each member contains an individual simulation. This feature was added to facilitate ensemble forecasts. In addition, each member is usually divided into different sequential chunks in order to save checkpoints of the model state at regular intervals. With these features, Autosubmit has the ability to run multiple members in parallel and therefore is suitable to run a GA in which there are different individuals in the same generation. This allows the size of the experiment to be adjusted easily and many different quantities of population and generations to be tested with ease. The use of Autosubmit to facilitate multiple instances of a computational model in a BRKGA is a novel one. One shortcoming of this method, however, is that the workflow size is static, and there is no feature to terminate the experiment after a certain condition is met. This means that the only viable stopping condition of the BRKGA is after a predetermined number of generations; otherwise, the stopping condition would have been if no evolution is observed after a certain number of generations.

Our particular workflow consists of three different types of job. The first is the initialisation of the experiment and is only run once at the very beginning of the experiment. The second is the simulation, which is run once per individual in parallel in each generation. Finally, the postprocessing, which includes the crossover, is run once per generation. An example of a workflow for a toy experiment of four populations and four generations is shown in Fig. 2.

Figure 2An example of the Autosubmit workflow.

Download

2.4.1 Initialisation

The initialisation script starts by setting up the directory in which the simulations are run by copying the executable of the model and the necessary input files into it. Included within the initialisation is a simulation run with a vector of the default parameters, and certain statistical measurements between its output and the observed data are taken that are necessary for postprocessing and calculation of the cost function. Finally, the script generates the initial set of vectors at random.

2.4.2 Simulation

The second script, which runs p times at each generation in parallel, starts by setting up the environment for each simulation. It then reads its corresponding vector from the generated set and edits the name lists to contain the updated parameters. Afterwards, the simulation is run and the cost function calculated.

2.4.3 Crossover

The final script runs once per generation after all simulations of the respective generation are completed. First, it reads the cost function statistics calculated after each simulation and uses them to calculate the ST score of SPOC and LPOC. It then ranks each of the simulations according to the sum of the two ST scores. Then it performs the crossover, as described in Sect. 2.3.1, to produce a new set of parameters in the same format so that it can be read by the following generation's simulation scripts.

2.5 Experiments

To investigate the potential of the BRKGA, different sets of experiments are run. Each set contains five experiments (to test consistency and robustness) with distinct and randomly generated initial populations, with 100 individual simulations over 100 iterations. Their details are summarised in Table 2.

Table 2A summary of the experiments run using the workflow.

Download Print Version | Download XLSX

Initially, we determine the capabilities of the BRKGA by testing how well it can find a known set of parameters. To do this, experiment sets D9 and D5 are run using the output of a simulation with the default parameters being the reference data at location LAB1. In set D9, nine parameters are tested to check which ones can be constrained from SPOC and LPOC data. This leads us to select five parameters, which are tested in set D5 and, additionally, give us an indication of how the method behaves when different sizes of vectors of parameters are used.

Experiment set O5_LAB1 uses the BRKGA as intended, where the reference data are observed data from LAB1, and the outputs are analysed. This is further compared with experiment sets O5_LAB2 and O5_LAB3, which are run in LAB2 and LAB3, respectively. This is to investigate how the results obtained reflect the wider region.

Finally, cross-simulations are run, whereby a representative vector of parameters from each of the experiment sets O5_LAB1, O5_LAB2, and O5_LAB3 is selected to run a single simulation in the other two locations. This is to further check how robust the BRKGA is and if the vectors produced are representative of the region. In fact, a certain homogeneity is expected across the three locations because of their similar physical and biogeochemical properties. The BRKGA not capturing this homogeneity would suggest that the tool is compensating for other errors in the attempt to minimise the cost function resulting in an overfitting of the optimal vector of parameters.

A set of experiments was designed to test the potential of a newly developed BRKGA. As a validation, the BRKGA was first tested against the output of a simulation produced with known default parameter settings. For the first set of experiments, we chose nine parameters, expressed as a vector, insuring a broad selection. This guided our choosing of parameters that could be constrained with confidence from the evaluated variables (in this case, SPOC and LPOC). In addition, in this set (and all others in the paper) five identical experiments were run at a time, and all results were similar to each other; this indicates that this method behaves consistently and reliably. The next set of experiments was identical to the previous set, except that there were only five parameters selected from the initial nine-parameter vector. This particular set of experiments produced results that were closer to the result of the default parameter vector with less computation. This leads us to believe that the size of the experiment required is dependent on the size of the parameter vector. One of the main contributions of this work is to use a state-of-the-art ocean model as a prior step to the calculation of the fitness function, with all the complexity that this option entails. This is only possible because of the aforementioned availability of computing power, and it is also highly facilitated by the usage of advanced scientific workflow solutions that allow the integration of the model executions in the evolutionary workflow (Oana and Spataru, 2016; Dueben and Bauer, 2018; Rueda-Bayona et al., 2020).

After the experiments against default data, the BRKGA was then tested by using observed data from ocean floats in the North Atlantic as the reference data. A set of five experiments was run for each BGC–Argo float annual time series, using the same settings as in the previous set. From Fig. 9, we can see that the results show a similar level of consistency to those with the default data. There is a visible improvement in the outputs of the simulations that use a set of parameters that have been optimised by the BRKGA compared to the outputs with the default parameter simulations (Figs. 12 and 13). However, most of the optimised parameter values tended rapidly towards the optimisation bounds (Table 1), and some even exceeded them thereafter because parameters were allowed to exceed the bounds by a small percentage in each generation. This behaviour makes us question whether the optimised values are realistic, although it is also possible that we imposed too-strict bounds in some cases, given the wide range of plausible ranges that characterises some parameters (see below). The problem of obtaining a right answer for the wrong reasons is common to all PO methods when applied to complex and heavily parameterised systems (Löptien and Dietze, 2019; Kriest et al., 2020). Therefore, PO must always be followed by a critical evaluation of the results. If a parameter converges repeatedly to unrealistic values, regardless of the value of other parameters, this may indicate that a process is poorly represented by the model equations. In such cases, PO can prompt further model development.

Another concern that arises from the results is the need to carefully evaluate the behaviour of the cost function. This is well illustrated by Figs. 4 and 5, which show that, on occasion, some statistics were improved at the expense of others, for example, bias at the expense of correlation. Correctly balancing bias, variability, and pattern (correlation) statistics in the cost function is critical to obtain meaningful PO results. Traditional cost functions based solely on the RMSE tend to reward solutions with too-low variability, whereby the positive biases cancel out the negatives (Jolliff et al., 2009, and references therein). Cost functions such as the ST score used here (Jolliff et al., 2009) were designed to avoid this problem. However, their behaviour is also sensitive to the overall variability and the signal-to-noise ratio in the data. Our preliminary tests suggested slightly better BRKGA performance after the logtransformation of the data. This procedure reduced the weight of the very high POC concentrations present only in the surface layer in spring–summer, favouring the representation of the portions of the water column with lower POC (i.e. the mesopelagic). Unlike the model outputs used to test the BRKGA in the first set of experiments, the BGC–Argo POC estimates are noisy. Therefore, the cost function may have been less effective when faced with the observed data. The performance of the cost function could also be improved by applying different weights to each variable, in this case SPOC and LPOC, which is common practice when the reference variables exhibit very different variability ranges (Friedrichs et al., 2007; Ayata et al., 2013).

Further work quantifying the effectiveness of the cost functions across different situations would probably improve the efficacy of the BRKGA. Yet, it must be highlighted that the test case chosen to evaluate the BRKGA is an exigent one because model skill was already very good with the default parameters, even though PISCES was not originally tuned to fit these particular observations. Ongoing work with a different optimisation case indicates that the BRKGA can produce larger and simultaneous improvements in all skill metrics when starting from a state of very poor model performance, which is, in this case, the seasonal cycle of sea surface chlorophyll a in the Tasman Sea ( Joan Llort,, personal communication, 2022; data not shown). Therefore, the tradeoffs between skill metrics observed here during the evolution of the experiments may indicate that the optimisation was operating close to the best skill attainable with a given set of model equations and considering observational uncertainty.

As a further test of our approach, two more sets of experiments were carried out in different locations in the Labrador Sea, resulting in a reasonable consistency of the optimal parameter set across the region. To confirm this, the optimal parameter sets for the three locations were cross-referenced by using each parameter set in each of the other two locations in single simulations. Results from this cross-testing suggest that the parameters produced have the potential to be representative of the region or even exchangeable among multiple locations (Table 3), meaning that the BRKGA is not compensating for other biases (e.g. physics) by overfitting (Löptien and Dietze, 2019; Kriest et al., 2020). This is an important aspect because it means the BRKGA could be used to investigate the large-scale spatial variability in key biogeochemical parameters. In particular, in the past decade, several authors have investigated the spatial variability in the transfer efficiency of POC from the surface ocean to its interior using different approaches and arriving at contrasting conclusions (Henson et al., 2011; Marsay et al., 2015; Guidi et al., 2015; Weber et al., 2016; Schlitzer, 2004; Wilson et al., 2015). Such spatial variability would be very difficult to establish in a three-dimensional framework because of the high computational cost required. This is an example of the still-open scientific questions that could be tackled with our approach.

The optimisation of PISCES parameters against BGC–Argo presented in our study illustrates how PO can help us understand a dynamical system better. Here we will briefly discuss the lessons learnt from the O5 experiments, while keeping in mind that a detailed review of PISCES parameter values and their biogeochemical implications is beyond the scope of this paper. It is also noteworthy that the interpretation provided here draws only from the analysis of the best-performing parameter set in each BRKGA experiment. Full exploitation of the results, with thousands of alternative model realisations, could yield further insights on how parameters interact in a space constrained by optimal model performance.

In the three O5 experiments, wsbio converged to values between 0 and 1 m d⁻¹. The decrease in the POC sinking speed improved the fit to observations by reducing the plumes of sinking POC that formed below intense phytoplankton blooms in the simulations (Fig. 12). Galí et al. (2022) identified these plumes as being the main reason for SPOC model–data misfit in the upper mesopelagic in several subpolar locations in the Northern and Southern hemispheres. The decrease in wsbio effectively turned the SPOC fraction into suspended POC, which is plausible according to field studies that sorted POC fractions according to their sinking speed (Riley et al., 2012; Baker et al., 2017). The evolution of the remaining parameters acted to adjust the magnitude and shape of POC vertical profiles. Increased xremip implies a steeper vertical decrease in both POC and GOC, with a stronger effect on POC given its much longer residence time in the mesopelagic. The xremip parameter represents the maximal specific remineralisation rate attainable in the model, corresponding to freshly produced detritus, normalised to a temperature of 0 ^∘C with a power law temperature dependence (Aumont et al., 2015). Our optimised xremip range, 0.078–0.114 d⁻¹, is consistent with the median of the highest values found across six field and laboratory studies when normalised to 0^∘ in the same way, i.e. 0.10 d⁻¹ (Belcher et al., 2016, and references therein), with an absolute maximum of 0.135 d⁻¹ (Iversen and Ploug, 2010). Decreased POC sinking speed and increased remineralisation would deplete mesopelagic POC during the productive season (central panel of Fig. B1) and, hence, SPOC if they were not compensated by other processes. In our PO experiments, this deficit was compensated by slightly increased surface microphytoplankton mortality and aggregation (wchld and wchldm), which supply SPOC and LPOC.

Interpretation of the evolution of the grazflux parameter is more complex. The flux-feeding rate depends on the product of mesozooplankton biomass, grazflux, and particle sinking speed. In PISCES, a fraction of the intercepted GOC is fragmented into POC. Therefore, flux feeding acts to remove POC and GOC (preferentially the fast-sinking GOC) and simultaneously produce POC; this process is becoming an important POC source in the lower mesopelagic (Figs. A1 and A2). Although a decrease in grazflux provided the best fit to observations, an increase in grazflux could also improve model skill, as shown by the relatively good skill of experiment a3nj during the first few generations (Figs. 9 and 10). This dual behaviour is confirmed by sensitivity analyses (Fig. B3) which show that, for a given xremip, a decrease in POC (which improves the model–data fit in the upper mesopelagic) can be achieved by either increasing or decreasing grazflux. Overall, these findings provide further evidence for the difficulty of constraining this important parameter (Jackson, 1993; Stemmann et al., 2004; Gehlen et al., 2006; Stukel et al., 2019).

The selection of a subset of model parameters is a common limitation of PO experiments, and although we based it on objective criteria, we acknowledge that it remains somewhat arbitrary. The stepwise reduction in the number of parameters from nine to five obeys the need to assess the GA performance with a varying number of parameters and also to reduce the degrees of freedom, given that only two variables were used as reference observations. Among the excluded parameters, wsbio2 certainly deserves examination in future experiments, given its primary control on the fate of GOC (large detritus). There are the following three main reasons that led us to exclude wsbio2 from this work: (1) our optimisation exercise focused on POC stocks, which are largely dominated by the SPOC fraction that typically represents around 85 % of total POC (Galí et al., 2022, and references therein), (2) experiment set D9 suggested that, unlike wsbio, wsbio2 might be difficult to constrain even from error-free data provided by the default model run, and (3) estimation of wsbio2 relies, more than any other parameter, on observational LPOC estimates which suffer from larger uncertainty than SPOC estimates (Galí et al., 2022).

Mesopelagic POC dynamics provide a relevant optimisation case because of their role in oceanic carbon sequestration (Martin et al., 2020; Henson et al., 2022). Beyond the specific results presented here, our approach has the additional value of showing that similarly good fits can be obtained with different combinations of parameter values (Ward et al., 2010), pointing to the many degrees of freedom in PISCES and similar OBGCMs. Hence, we call for a continuous reassessment of parametric uncertainty as new types of observations (e.g. from BGC–Argo) become available. Future work is granted to study the sensitivities, interdependencies, and optimal values of PISCES parameters through more comprehensive experiments.

The code of NEMO v4.0.1 and PISCES-v2_RC are publicly available at https://www.nemo-ocean.eu/ (last access: April 2021), https://doi.org/10.5281/zenodo.3878122 (Madec et al., 2022). The PISCES 1D configuration used in this study is available at https://earth.bsc.es/gitlab/mgalitap/p1d_share/-/tree/gapoc/ (last access: May 2021) and https://doi.org/10.5281/zenodo.5243343 (Galí et al., 2021b). The code for the workflow of the genetic algorithm is readily available at https://doi.org/10.5281/zenodo.5205760 (Falls et al., 2021).

These data were collected and made freely available by the international Argo programme and the national programmes that contributes to it (https://argo.ucsd.edu, last access: April 2021; https://www.ocean-ops.org, last access: May 2021). The Argo programme is part of the Global Ocean Observing System. The repository that contains the processed BGC–Argo data and matching PISCES 1D simulations, titled “Datasets for the comparison between POC estimated from BGC–Argo floats and PISCES model simulations”, is readily available at https://doi.org/10.5281/zenodo.5139602 (Galí et al., 20211).

MF wrote the paper, with contributions from all co-authors. MC and MA contributed to the topics of computing and GAs. MG, RB, and JL contributed to the topics of ocean biogeochemistry. MF developed the code of the workflow of the GA (including the implementation of the GA itself and the configuration of the workflow manager) and developed the code to produce the figures and the data in the tables. MC developed the configuration of NEMO-PISCES for the MareNostrum 4 HPC platform and provided guidance to maximise the workflow efficiency. MG configured the PISCES 1D offline simulations and processed the observed data from the BGC–Argo floats. MG and RB conceived the study.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This research has been supported by the Fundación Bancaria Caixa d'Estalvis i Pensions de Barcelona (grant no. LCF/BQ/PI18/11630009) and the Ministerio de Ciencia e Innovación (grant nos. PID2019-107952GA-I00 and CGL2017-84493-R).

This paper was edited by Olivier Marti and reviewed by Urmas Raudsepp and one anonymous referee.