1 INTRODUCTION

The deep chlorophyll maximum (DCM) is a subsurface peak in chlorophyll a, and refers to the equilibrium of phytoplankton growth with opposing resource gradients, including light from the surface and nutrients from below (Huisman et al., 2006; Harrison and Smith, 2011). The DCM is a widespread ecological phenomenon that exists in most oceans and lakes (Fee, 1976; Barbiero and Tuchman, 2001). Previous studies have demonstrated that the DCM has a close relationship with climate processes (Boyce et al., 2010; Doney et al., 2012), biogeochemical cycles (Brandini et al., 2014; Laufkötter et al., 2015), and particularly the carbon cycle (Paquette and Beisner, 2018). To understand the dynamics of DCM formation, numerical models describing marine ecosystems have been widely used (e.g., Christian et al., 2002; Heinle and Slawig, 2013).

Ocean ecosystem models usually involve multiple parameters that represent various physical and biochemical processes (Fulton et al., 2003; Kantha, 2004). The uncertainties of these parameters are primary sources of discrepancies between modeled DCM and observations (Losa et al., 2004). However, accurate estimates of these parameters are challenging due to the following reasons. First, parameters that contribute to ocean ecosystem simulations tend to change under different environmental conditions (Yoshiyama et al., 2009). For example, although there are many different phytoplankton species in the ocean ecosystem, the parameters associated with these species are commonly set as constants. Furthermore, it is unrealistic to expect to improve the accuracy of all parameters by field observations because of the large number of uncertain parameters and the high cost of ocean observations. Therefore, it is necessary to identify important parameters (i.e. sensitive parameters) that cause relatively large simulation or prediction errors. Improving these priority parameters can lead to better simulations or predictions of the ocean ecosystem.

To-date, there have been many proposed approaches to identify sensitive parameters. Specifically, some of them employ the data assimilation technique (e.g., Fiechter, 2012; Weir et al., 2013; Schartau et al., 2017). This strategy seeks to minimize the difference between temporally varying modeled and observed data by adjusting the value of control parameters. For instance, Evensen et al. (1998) constructed a penalty function for estimating poorly known parameters and found this method was only feasible for problems involving weak nonlinearity. Fennel et al. (2001) estimated parameter sensitivity based on posterior errors built with a Hessian matrix. Their study emphasized the importance of parameter optimization in testing the reliability of dynamic models. In contrast to data assimilation, some other approaches employ a variance-based methodology (Saltelli et al., 2005; Chu et al., 2007). This strategy tries to identify sensitive parameters by assessing their contribution to output variability. By utilizing special metrics of warping functions and irreversible predictability, Chu et al. (2007) explored the nonlinear sensitivity of model outputs to parameters and found this strategy feasible for parameters with an extremely wide range. Nonetheless, this approach requires a large quantity of computing resources, especially to investigate the sensitivities of multiple parameters. In recent years, Mu et al. (2010) proposed a new approach, conditional nonlinear optimal perturbation related to parameters (CNOP-P), to investigate the impact of parameter uncertainties with a nonlinear regime. The CNOP-P approach seeks to find the most sensitive parameter or parameter combination that induces the largest effect on simulating/predicting investigated events. This approach saves computing power and simultaneously considers nonlinear interactions among different parameters. To-date, the CNOP-P approach has been successfully utilized in parameter sensitivity analyses of net primary productivity in eastern China (Sun and Mu, 2014), predictions of soil carbon (Sun and Mu, 2017), and double-gyre variations (Yuan et al., 2019).

In this study, we investigated parameter sensitivity on simulating the DCM using the CNOP-P approach. By employing a nutrients-phytoplankton ecosystem model, we first identified the sensitivity of each parameter and then investigated the sensitivity of combinations that involves different parameters. The rest of this paper is organized as follows. A brief description of the nutrients-phytoplankton ecosystem model and the CNOP-P approach are presented in Section 2. Sensitivity studies of single parameter and parameter combination are discussed in Section 3. The improvements in the DCM simulation caused by reducing parameter uncertainties are presented in Section 4. Finally, Section 5 presents the summary and discussion.

2 MODEL AND APPROACH 2.1 A nutrients-phytoplankton ecosystem model

In this paper, the DCM was simulated by a onedimensional ecosystem model, which is able to reproduce most ecological features in the upper ocean (e.g., Huisman et al., 2006; Löptien, 2011; Liccardo et al., 2013). This model describes phytoplankton and nutrient dynamics limited by the supply of light and nutrients. The governing equations of phytoplankton Ph (cells/m³) and nutrients N (mmol/m³) are given by:

where z is the water depth, t is the integration time, m is the loss rate of phytoplankton, κ is the coefficient of vertical turbulent diffusivity, α is the nutrient content of phytoplankton, ε is the nutrient recycling coefficient, ν is the phytoplankton sinking velocity; I is light intensity; and H_N and H_I are the half-saturation constants for nutrient-limited and light-limited growth, respectively.

Due to attenuation from concentrated phytoplankton and background turbidity, incident light intensity decreases with water depth. Using Lambert-Beer's equation (Kirk, 1994), the expression for the incident light intensity can be written as

where I_in is the incident light intensity at the top of the water column; k and K_bg are the attenuation coefficients of background turbidity and concentrated phytoplankton, respectively. As shown in Table 1, the values of all parameters are consistent with previous studies (Liccardo et al., 2013; Wang and Mu, 2015).

Regarding boundary conditions, the fluxes of phytoplankton and nutrients are

where Z_B denotes the bottom of the water column.

By employing a forward difference scheme, we discretized Eqs.1 & 2 in resolution of 0.5 m (600 vertical grids in total) and integrated the model for 20 years with a time step of 0.01 h. As the depleted nutrients and attenuated light for phytoplankton growth achieve a balance below the surface (Fig. 1a & b), the system gradually reaches a stable DCM equilibrium. We present the vertical nutrient and phytoplankton profiles of the 20^th year in Fig. 1c & d. The maximum phytoplankton biomass is located at ~85 m, consistent with observational results and model simulations (Navarro and Ruiz, 2013; Wang and Mu, 2015).

2.2 CNOP-P approach and its application 2.2.1 CNOP-P approach

Mu et al. (2010) proposed the CNOP-P to investigate the effects of uncertainties of single parameter and parameter combination consisting of multiple parameters. In this study, the optimal single parameter or parameter combination revealed by the CNOP-Ps is regarded as the most sensitive parameter(s), which leads to maximal uncertainty in numerical simulations. The detailed derivation of the CNOP-P approach is introduced next.

A general nonlinear evolution equation can be expressed as

where t represents time and ^F is a nonlinear propagator. U∈Rⁿ is the state variable, U₀ is its initial value, and p=(p₁, p₂, p₃, ∙∙∙, p_n) is the vector consisting of all model parameters. Let M_t be the nonlinear propagator of the differential equations from the initial time 0 to t, then U(t)= M_t(U₀, p) is the solution of Eq.6 at time t. U(t; U₀, p) and U(t; U₀, p)+u(t; U₀, p′) are the solutions for vectors p and p+p′, where u(t; U₀, p′) represents the departure from the reference state U(t; U₀, p) caused by parameter perturbation p′. They satisfy the following equations

For a given norm ||•||, a parameter perturbation p_Ω′ is called a CNOP-P if and only if

where J is the cost function and Ω is the constraint of parameter perturbation. Essentially, CNOP-P (p_Ω′) is the optimal parameter perturbation causing the largest variation of the cost function.

2.2.2 Settings for calculating the CNOP-P

In this study, the parameter vector consists of 10 parameters (Table 1). These parameters are related to major biological and physical processes, including nutrient (p₁) and light (p₂, p₈, p₉, p₁₀) limitations that affect the growth of phytoplankton, loss of phytoplankton from death (p3) and gravity (p₄), turbulent mixing of water (_p5), and internal nutrient supply (p₆, p₇). We denote them in a vector p=(p₁, p₂, p₃, ∙∙∙, p₁₀). Previous studies have shown that phytoplankton biomass is a good proxy for DCM intensity (Klausmeier and Litchman, 2001; Wang and Mu, 2015). Thus, the objective function which is employed to seek the maximal change in DCM intensity is defined as

where ||•|| is the L₂ norm and Ph′ is the change in phytoplankton caused by parameter perturbation p′. In addition, the equilibrium of nutrients and phytoplankton presented in Fig. 1c & d is taken as the reference state to calculate the CNOP-P. The optimization time (t=1 080 d) was selected based on the diffusive timescale (Beckmann and Hense, 2007; Wang and Mu, 2015). This optimization time can ensure the ecosystem reaches a new equilibrium after the model is superimposed with parameter perturbations. We also verified that the parameter sensitivity results remain unchanged when the optimization time is changed to 720 d or 1 440 d. In addition, the differential evolution optimization algorithm, which displays strong global convergence and robustness (Storn and Price, 1997), is used to obtain the CNOP-P.

3 SENSITIVITY ANALYSIS OF SINGLE PARAMETER AND MULTIPLE PARAMETERS 3.1 Sensitivity analysis of single parameter

For each of the 10 selected parameters, we use the CNOP-P approach to calculate its optimal perturbation and examine its impacts on DCM simulations. Specifically, we employ p_i′ to denote the perturbation of p_i. To examine whether the sensitivity of single parameter depends on the constraint of parameter perturbation, different constraint values are selected with the following equation

where λ is the factor modulating the amplitude of perturbation constraint. To ensure a reasonable DCM position of 50–100 m (e.g., Huisman et al., 2006; Liccardo et al., 2013), λ≤0.4 is selected. Moreover, different values of λ, which ranges from 0.1 to 0.4 at an interval of 0.1, are used to explore the impacts of different perturbation constraints on the CNOP-P calculations. Hence, to study the sensitivity under different constraints, four different optimization experiments are conducted for each parameter.

Figure 2 demonstrates the sensitivity of the 10 parameters according to the cost function values caused by the CNOP-Ps. For different constraints, the parameter sensitivity rankings are similar. Moreover, the CNOP-P for each parameter is located at the boundary of the perturbation constraint. Here, we focus on the four most sensitive parameters since they induce relatively large changes in the DCM. The sensitivity ranking among the four most sensitive parameters basically remains unchanged for different parameter perturbation constraints. The most sensitive parameter is always p₉ (K_bg, background turbidity), which limits light supply in the model. The corresponding cost function value is more than four times as large as the average caused by the CNOP-Ps related to the other parameters. Furthermore, p₆ (α, nutrient content of phytoplankton), p₇ (ε, nutrient recycling coefficient), and p₅ (κ, vertical turbulent diffusivity) are the other three most sensitive parameters. The cost function values caused by the CNOP-Ps of these three parameters are about three times as large as the average caused by the six nonsensitive parameters.

Parameter p₉ affects the light supply from the surface, whereas the other three sensitive parameters are related to nutrient supply. Therefore, superimposing these parameter perturbations affects the DCM equilibrium between light and nutrient supplies. To investigate how these four sensitive parameters specifically affect the DCM, the DCM changes caused by superimposing the CNOP-Ps of the sensitive parameters were examined. For brevity, only the CNOP-Ps obtained with λ=0.3 are discussed here, which induce corresponding perturbations of p₉′=-0.015, p₆′=-3.0×10¹⁰, p₇′=0.15, and p₅′=0.12.

Reduced background turbidity (p₉′) allows light to reach deeper into the water column, leading to the rapid phytoplankton growth and biomass increase in the lower euphotic zone. At the same time, phytoplankton consumes more nutrients and thereby limits the phytoplankton growth in the upper layer (Fig. 3a & e). Eventually, a one-third increase in phytoplankton biomass and a deeper DCM are found after integration of 360 d (Fig. 3i). Additionally, the enhanced vertical turbulent diffusivity (p₅′) and nutrient recycling coefficient (p₇′) both increase the nutrient supply (Fig. 3g & h). The decline in the nutrient content of phytoplankton (p₆′) indirectly results in a greater nutrient supply because each cell requires less nutrients under this condition (Fig. 3f). The abundant nutrient supply first benefits the phytoplankton growth in the upper euphotic zone (Fig. 3b–d). The concentrated phytoplankton in the upper water column subsequently prevents light penetration, further leading to a decrease in phytoplankton in the lower euphotic zone. This results in intensified and shallower DCMs (Fig. 3j–l).

3.2 Sensitivity analysis of multiple parameters

Despite the analysis of single parameter sensitivities, it is important to remember that in ocean ecosystem models, the physical and biological processes represented by various parameters do not exist independently. Interactions among parameters play an important role in the identification of sensitive parameters (Chu et al., 2007; Sun and Mu, 2016; Wang et al., 2017). Thus, in this study, the impacts of multiple parameter perturbations on DCM simulations are further investigated by calculating the CNOP-Ps of parameter combinations.

Studies on parameter combinations can reflect the interactions among multiple parameters. Particularly, we aim to figure out whether the optimal parameter combination is a simple combination of the top four sensitive parameters (p₅, p₆, p₇, p₉). If yes, are their CNOP-Ps the same to those obtained from the single parameter analysis? If not, how do the interactions among various parameters affect the identification of sensitive parameters? To answer these questions, we examined the sensitivity of 210 (C₁₀⁴) groups of parameter combinations. Note that only the perturbation constraint of λ=0.3 is investigated for brevity. In addition, we have verified that the most sensitive parameter combination does not change at λ=0.1, 0.2, and 0.4. Furthermore, 16 CNOP-Ps of parameter combinations do not reach the perturbation constraint boundaries. This may be attributed to the nonlinear interactions among different parameters in a combination.

Figure 4 demonstrates the cost function values caused by the CNOP-Ps of the 210 groups of parameter combinations. The maximum cost function value occurs when the parameter combination of (p₃, p₅, p₆, p₇) is perturbed. It is nearly 2.5 times as large as the average caused by the CNOP-Ps of all parameter combinations. It is noted that the optimal parameter combination is different with (p₅, p₆, p₇, p₉). The cost function caused by the perturbation of (p₅, p₆, p₇, p₉) only shows the 13^th largest value. Except for p₃ (m, rate of phytoplankton loss) and p₉, the other three parameters (p₆, p₇, p₉) are the same in these two combinations.

To further validate the importance of p₃, Table 2 lists the frequency of each parameter in the 10 most sensitive combinations. We found that p₃ occurs eight times, whereas p₉ only occurs twice. Thus, p₃ plays an important role in the identification of sensitive parameter combinations. However, with regard to single parameter sensitivity, p₉ is the most sensitive and p₃ is not sensitive at all (Fig. 2). This raised the question why p₃ displays such importance in the optimal parameter combination. Are the nonlinear interactions among different parameters responsible for this?

To answer these questions, we superimposed the CNOP-P of the optimal parameter combination (p₃′=-0.003, p₅′=0.12, p₆′=-3.0×10^-10, p₇′=0.15) where p₃′ represents a decrease in the rate of phytoplankton loss, leading to an increase in phytoplankton. The nutrient supply is essential to maintain this significant increase. In this case, p₅′, p₆′, and p₇′ provide abundant nutrients, as discussed in Section 3.1. Overall, these four types of perturbations eventually lead to nutrient increases (Fig. 5c), which results in increased phytoplankton biomass and a shallower DCM (Fig. 5a & e). However, no additional nutrients are provided when only p₃ is perturbed. Therefore, p₃ is not sensitive as a single parameter.

In contrast, p₉ shows less sensitivity in the sensitivity analysis of multiple parameters. The optimal perturbation (p₅′=0.12, p₆′=-3.0×10^-10, p₇′=0.15, and p₉′=-0.005 66) was superimposed on the reference state. The influence of p9′ on the DCM is opposite to the other three parameters (Fig. 3a). Specifically, p₉′ induces a deeper DCM, whereas the DCM moves to a shallower depth with p₅′, p₆′, and p₇′. Therefore, there is an opposite effect when the model is simultaneously perturbed with p₅′, p₆′, p₇′, and p₉′. As such, phytoplankton showed less of an increase compared with that caused by the optimal perturbation combination (p₃′=-0.003, p₅′=0.12, p₆′=-3.0×10^-10, p₇′=0.15) and p9 shows less sensitivity (Fig. 5b & f).

The interactions among different parameters present a remarkable influence on the identification of sensitive parameters. For the ocean ecosystem, there are effects that enhance and cancel the physical and biological processes represented by different parameters. The sensitivity analysis of single parameter only evaluates the importance of one parameter and so cannot consider enhancement or cancellation effects. In contrast, the sensitivity analysis of multiple parameters with the CNOP-P approach can reflect different physical or biological processes (especially important processes) at the same time. Hence, sensitive parameters identified by the CNOP-Ps of parameter combinations should be more effective at improving the performance of DCM simulations.

4 IMPACT OF REDUCING PARAMETER UNCERTAINTIES

One of the most important objectives in identifying sensitive parameters is to reduce the uncertainties in numerical simulations or predictions. Errors in sensitive parameters are expected to be corrected through observations or laboratory experiments. Thus, we investigated whether DCM simulations can be significantly improved by eliminating uncertainties in sensitive parameters. Specifically, the sensitive parameters identified in the sensitivity analyses of single parameter and multiple parameters are examined, respectively.

In most cases, we are most concerned about parameter perturbations that can lead to large uncertainties in numerical simulations (Sun and Mu, 2016). Thus, using the same perturbation constraint (λ=0.3), the CNOP-P of the 10 parameters was calculated. As a result, p_all′=(-0.007 5, -12, -0.003, -0.012 6, 0.12, -3.0×10¹⁰, 0.15, -2.4×10¹⁰, 0.006 652, 60) was obtained, which lead to the largest DCM simulation errors. The uncertainties of the sensitive parameters, identified by the CNOP-Ps of single and multiple parameters, were reduced based on p_all′. Specifically, the uncertainties of the parameters p₃, p₅, p₆, and p₇, identified with the optimal parameter combination, were first reduced by utilizing the controlling factor δ₁ defined in Eq.11.2 (denoted as p_case1′). For comparison, we also reduced the uncertainties in the other two parameter combinations (denoted as p_case2′ and p_case3′). According to δ₂ (Eq.11.3), p_case2′ reduces uncertainties of the four sensitive parameters obtained from the single parameter analysis. Additionally, p_case3′ reduces the uncertainties of the remaining parameters, except for the optimal parameter combination, controlled by δ₃ (Eq.11.4).

The improvement extent τ for uncertainties in DCM simulations is quantified by

where μ is the ratio modulating the extent of reductions of parameter perturbation (varying from 0.1 to 0.4 at an interval of 0.1) and integration of t=360 d.

Figure 6 illustrates the improvements in DCM simulations after reducing the uncertainties of the three groups of parameters. Better DCM simulations were consistently achieved with larger values of μ. However, for different values of μ, reducing perturbations of p₁, p₂, p₄, p₈, p₉, and p₁₀ produced no significant improvements. This indirectly reflects the importance of identifying sensitive parameters by the CNOP-P approach because most of these six parameters are not included in the CNOP-Ps of either sensitive single or multiple parameters. More importantly, for cases with the same value of μ, reducing parameter errors from the most sensitive parameter combination (p₃, p₅, p₆, p₇) always produces the best model performance. For different values of μ, improvements in reducing perturbations from the most sensitive parameter combination are 10%–25% better than those caused by the sensitive parameters p₅, p₆, p₇, and p₉. In particular, when μ is 0.4, the improvement induced by this group of parameters achieves a high value (~76%). Therefore, focusing on sensitive parameters identified by CNOP-P of single and, especially, multiple parameters is very effective in improving DCM simulations.

5 SUMMARY AND DISCUSSION

Using a nutrients-phytoplankton ecosystem model and the CNOP-P approach, this study explored the impacts of parameter uncertainties on simulating the DCM. First, the DCM was successfully simulated by the ocean ecosystem model. Then, the sensitivity of single parameter was investigated by calculating the CNOP-Ps related to 10 individual parameters. A sensitivity ranking of the 10 parameters was obtained according to the cost function values caused by the CNOP-Ps. The top four sensitive parameters were further analyzed because they were shown to induce relatively large changes in the DCM. Specifically, the most sensitive single parameter is p₉ (K_bg, background turbidity). The CNOP-P of p₉ enhances the light supply for the ecosystem, and thereby leads to increased phytoplankton and a deeper DCM. The other three sensitive parameters are p₆ (α, nutrient content of phytoplankton), p₇ (ε, nutrient recycling coefficient), and p₅ (κ, vertical turbulent diffusivity). By enhancing the nutrient supply, the CNOP-Ps of these three parameters all lead to increases in phytoplankton and upward moves of the DCM.

Considering the importance of interactions among different parameters, the sensitivity of multiple parameters was further investigated by calculating the CNOP-Ps related to combinations of four parameters. As a result, the optimal parameter combination (p₃, p₅, p₆, p₇) was identified. Interestingly, the parameters in the optimal parameter combination are not the same as the top four sensitive parameters obtained from the single parameter analysis. In the optimal parameter combination, p₉ is replaced by p₃ (m, loss rate). Further investigation indicates p₃ becomes more important because of the abundant nutrient supply caused by perturbations of p₅, p₆, and p₇. In contrast, there is a cancellation effect between the perturbation of p₉ and the other three parameters. This cancellation effect also limits p₉′ to a boundary value.

Finally, we quantified the improvement in DCM simulations by reducing the uncertainties of the identified sensitive parameters. Compared with nonsensitive parameters, reducing uncertainties of the sensitive parameters always exhibit better performances. Moreover, the sensitive parameters in the optimal parameter combination are more effective than those obtained from the single parameter analysis. This again emphasizes the importance and necessity of considering the interactions of various parameters for parameter sensitivity analysis.

The CNOP-P approach is effective in identifying sensitive single parameters and parameter combinations. Our results highlight the importance of considering interactions among different parameters for the identification of sensitive parameters. The CNOP-P approach related to parameter combination can identify sensitive parameters more effectively, especially for models that involve a large number of parameters. Reducing the uncertainties of these sensitive parameters is an effective and economical way to improve numerical simulations or predictions. CNOP-P approach related to parameter combination can identify sensitive parameters more effectively, especially for models that involve a large number of parameters. Reducing the uncertainties of these sensitive parameters is an effective and economical way to improve numerical simulations or predictions. Moreover, the CNOP-Ps provide a potential way to understand the dynamics of the concerned ecological phenomenon, in this case of the DCM.

It should be also noted that nonlinear interactions of other numbers of parameters are worth investigating. More experiments would be required to investigate other CNOP combinations with different numbers of parameters. One of our future goals is to determine the best parameter combination that optimally improves DCM simulations with regard to effectiveness and efficiency. Besides, a detailed analysis on the nonlinearity is required next. Overall, our results are preliminary but encouraging. It is expected that future ecosystem simulations could greatly benefit from the CNOP-P approach. To achieve this goal in an operational framework, substantial work will be required.

6 DATA AVAILABILITY STATEMENT

All datasets generated and analyzed during this study are available from the corresponding author on reasonable request.