A new optimization strategy is proposed to identify the sensitivities of simulations of atmospheric and oceanic models to uncertain parameters. The strategy is based on a nonlinear optimization method that is able to estimate the maximum values of specific parameter sensitivity measures; meanwhile, it takes into account interactions among uncertain parameters. It is tested using the Lorenz’63 model and an intermediate complexity 2.5-layer shallow-water model of the North Pacific Ocean. For the Lorenz’63 model, it is shown that the parameter sensitivities of the model results depend on the initial conditions. For the 2.5-layer shallow-water model used to simulate the Kuroshio large meander (KLM) south of Japan, the optimization strategy reveals that the prediction of the KLM path is insensitive to the uncertainties in the bottom friction coefficient, the interfacial friction coefficient, and the lateral friction coefficient. Rather, the KLM prediction is relatively sensitive to the uncertainties of the reduced gravity representing ocean stratification and the wind stress coefficient.
Numerical models have become an important tool for studying and predicting atmospheric and oceanic states. Generally, there are many uncertain parameters used to parameterize unresolved physical processes in atmospheric and oceanic models. The uncertainties of these parameters can add significant uncertainties in the simulation of atmospheric and oceanic models (e.g., Murphy et al. 2004; Edwards and Marsh 2005; Brierley et al. 2010; Levine-Moolenaar et al. 2012). Furthermore, the effects of uncertainties in different parameters on the model simulation and prediction are very diverse (Chu 1999; Orrell 2003; Yin et al. 2014). Hence, it is of importance to perform sensitivity analyses on model parameters and then identify which parameters have larger influences on the model results (the sensitive parameters) compared to others (the insensitive parameters). By excluding the insensitive parameters, the number of model parameters that need accurate estimation is reduced. This offers significant benefits in the assessment of model responses to parameter variations, data assimilation, model tuning, and so on (e.g., Cacuci 2003).
Various approaches have been developed to analyze the sensitivities of model outputs to parameters. The one-at-a-time (OAT) method is a direct and common method (Hamby 1994; Zaehle et al. 2005; Saltelli et al. 2010). The procedure varies one of model parameters at a time while fixing other parameters, then rerunning the model, and monitoring the changes in the model results studied. The parameters inducing the largest changes are then considered to be the most sensitive. The OAT method, however, ignores the interactions among uncertain parameters. Another option is the adjoint sensitivity method (Errico 1997), which has often been applied to analyze the parameter sensitivities in atmospheric and oceanic models, for instance, by Hall et al. (1982), Janisková and Morcrette (2005), and Daescu and Todling (2010). However, this method is based on linear approximations, which may be only valid for sufficiently small parameter perturbations and the development of model trajectories over a very short period. In this situation, this approach may not consider nonlinear interactions among model parameters.
Realizing the limitations of the above two approaches, Chu et al. (2007), Marzban (2013), and Marzban et al. (2014) applied the variance-based sensitivity analysis (VSA) method to investigate the model parameter sensitivity problem. Their results showed that the VSA method is able to capture the interactive role among model parameters. However, the VSA method only samples a limited number of points in the parameter perturbation space to examine parameter sensitivities. In addition, for multiple parameters in complicated atmospheric and oceanic models, the VSA method needs too many model simulations, which leads to computational feasibility issues.
The conditional nonlinear optimal parameter perturbation (CNOP-P) method proposed by Mu et al. (2010) can be used to assess the largest effects of model parameter uncertainties. The CNOP-P method is a constraint optimization method, which can find the parameter perturbations that cause the largest uncertainties among the model results. Based on the CNOP-P approach, Yin et al. (2014) and Sun and Mu (2016) explored model parameter sensitivities. However, the strategy for identifying the sensitive parameters in Yin et al. (2014) depends on the defined L2 constraint norm for parameter perturbations. Yin et al. (2014) judged the sensitivity of parameters according to the values of relative optimal parameter perturbations: the bigger the parameter perturbation is, the more sensitive the model output is to the parameter. However, for another constraint norm, such as a box constraint norm (i.e., L∞), the amplitudes of the optimal perturbations for different parameters may be the same. In this situation, it is impossible to determine parameter sensitivities using the strategy of Yin et al. (2014). For this reason, Sun and Mu (2016) calculated the CNOP-P with a box constraint norm to determine sensitive parameters. To reduce the computational cost, they first excluded the insensitive parameters from all the model parameters by calculating the single parameter CNOP-P (like the OAT method). During this process, simultaneous variations of parameters are not considered. In addition, to determine the combination of five sensitive parameters from 10 model parameters, Sun and Mu (2016) had to perform the nonlinear optimization calculations times. This reflects that the computational cost for their method is still too high. Hence, it is necessary to develop a complementary method of identifying the sensitivities of model parameters.
In this study, we propose a new optimization strategy to identify the parameter sensitivities in atmospheric and oceanic models. The new strategy is able to consider the nonlinear interaction role among model parameters and at the same time is computationally cheaper. Interestingly, we note that a method with the same aim was described by Oakley and O’Hagan (2004), but it was based on the Bayesian statistical analysis. In the next section, we will describe this new optimization strategy. Section 3 presents the test results of the new strategy when applied to the Lorenz’63 model (Lorenz 1963) and a 2.5-layer shallow-water ocean model used to simulate the Kuroshio path variations south of Japan. A summary and discussion are given in section 4.
2. The new optimization strategy
Generally, a numerical model can be formally written as
where X(t) is the state vector of model at time t, t(P) denotes a nonlinear propagator with a parameter vector P = (P1, P2, …, Pm), and X0 is the initial state vector.
This study mainly focuses on the uncertainties of model parameters. We therefore suppose that there exist no uncertainties in the initial conditions. In this situation, the amplitude of state perturbations in model output (simulation or prediction) at final time t caused by the perturbations of all parameters p = (p1, p2, …, pm) can be denoted as
where ||⋅|| indicates a chosen norm that measures the perturbation amplitude and P = (P1, P2, …, Pm) is regarded as the reference parameter vector. Similarly, the amplitude of the state perturbations induced by all the parameter perturbations except the ith perturbation pi, 1 ≤ i ≤ m (i.e., setting pi = 0), can be written as
Then, the magnitude of the state perturbations yielded by the parameter perturbation pi itself and its interaction roles with perturbations of other parameters is defined as
In general, the parameter uncertainties are limited, which means that they fall within some intervals. Therefore, the parameter perturbations are constrained by the specific conditions, that is, (p1, p2, …, pm) ∈ Cε. Of course, it is impractical to assess the effect of the parameter perturbation pi through taking different parameter perturbations and integrating the model repeatedly. Instead, we search for the maximum for the function Ji over the whole feasible region Cε by solving the following constraint optimization problem:
The above equation means that when the parameter perturbations take optimal values (p1,ε, p2,ε, …, pm,ε) that satisfy the given constraint condition (p1,ε, p2,ε, …, pm,ε) ∈ Cε, the uncertainty in the ith parameter and its interaction with uncertainties in the other parameters will lead to the maximal functional Ji,ε of model state perturbations at time t. This is to say, objective function value Ji,ε measures the maximal effect of the parameter perturbation pi, given the perturbations of the other parameters. Hence, for each model parameter perturbation, we can obtain the maximal effect by solving the optimization problem (5).
The quantities in (5) can be proved to satisfy the following inequalities:
The first inequality can be satisfied by setting all the parameter perturbations to zero except the ith parameter perturbation. This inequality indicates that the maximal effect of a single parameter perturbation is less than the objective function value Ji,ε, implying that the interaction among the model parameter uncertainties has been considered in (5). The second inequality is rather trivial considering (3) and (4). It means that the maximal effect of a single parameter perturbation pi and its interaction with the other parameter perturbations is always less than that caused by all of the model parameter perturbations.
The above discussion considers the maximal effects of uncertainty in a single parameter on the model outputs. Of course, we can also consider the situation for multiple parameters. Taking two parameters as an example, the objective function can be written as the following formula through setting pi = 0 and pj = 0:
Obviously, we can maximize Jij to explore the influence of the combination of pi and pj. In this case, however, the computational cost cannot be significantly reduced. In fact, we do not need to maximize Jij because we can easily prove the following inequality:
where Ji and Jj are defined by (4). The proof of this inequality is shown in the appendix. The inequality (8) indicates that when considering the interactions among parameters, the maximal joint effect of multiple parameter uncertainties is less than the sum of the maximal effect of each parameter uncertainty, which guarantees that when judging whether the parameters Pi and Pj are insensitive, we only need to calculate each term on the right-hand side of (8). If the value of each term is very small and their sum is less than a given value (that is defined for specific experiments in section 3), the parameters Pi and Pj can be considered to be insensitive. For more parameters, similar inequalities can also be obtained. This is the fundamental reason why our method can reduce the computational cost.
According to the above discussion, the optimization strategy for identifying the sensitivities of the model parameters can be described as follows: the optimization problem (5) for each model parameter is first solved and the maximum Ji,ε for each parameter is obtained; then the maxima are ranked from smallest to largest and we assume that the rank can be denoted as J1,ε ≤ J2,ε ≤ ⋅⋅⋅ ≤ Jm,ε; third, a criterion value denoted by Cr is defined according to a concrete physical problem; fourth, if J1,ε + J2,ε + ⋅⋅⋅ + Jk,ε ≤ Cr and J1,ε + J2,ε + ⋅⋅⋅ + Jk,ε + J(k+1),ε > Cr, the parameters corresponding to J1,ε, J2,ε, …, Jk,ε are regarded as being insensitive [because the maximal influence caused by the uncertainties in these parameters is still very small; see (5) and (8)] while the remaining parameters are relatively sensitive. For convenience, this optimization parameter sensitivity analysis method is abbreviated as OPSA. It is worth stressing that Cr is important in identifying the parameter sensitivity in the OPSA method. In realistic physical problems, Cr should be determined based on the understanding for the specific physical phenomena to be investigated.
The key for the OPSA method is to solve the optimization problem (5). The nonlinear optimization algorithm—the spectral projected gradient (SPG) algorithm (Birgin et al. 2000)—is employed to solve this problem. To make the SPG algorithm work, the information about the constraint condition, the objective function, and the gradient of the objective function with respect to each model parameter perturbation needs to be provided. The specific constraint condition and objective function norm for each model used will be given in section 3. In this study, for simplicity, the gradient information is calculated using the definition of the gradient, i.e., [Ji(p1, p2, …, pk + , …, pm) − Ji(p1, p2, …, pm)]/, k = 1, 2, …, m, where is a small perturbation. Of course, the gradient information can also be obtained using other methods, such as the adjoint method (Errico 1997). Furthermore, when performing the optimization calculation, we need to provide initial-guess values to the SPG algorithm. In the numerical experiments in section 3, 10 different random initial guess values are used in each optimization computation. In this way, we try our best to avoid having local maximal values of the objective function regarded as global ones. The specific details on the optimization calculations have been shown in Wang et al. (2012). It is worth pointing out that, to determine the sensitivities of m parameters in a model, the OPSA method only needs to perform the optimization calculation m times. This computational cost is thus much less than that in Sun and Mu (2016), in which the optimization calculations need to be performed m + times if one would like to determine n (1 < n < m) sensitive parameters.
3. Experimental design and results
In this section, we will identify parameter sensitivities using the OPSA method in section 2 for the simple Lorenz’63 model and an intermediate-complexity shallow-water ocean model.
a. Lorenz’63 model
To test the OPSA method, the Lorenz’63 model (Lorenz 1963) is considered to be a suitable tool. The model is frequently applied to test some new principles and methods because of its simplicity. For example, Chu (1999) employed this model to examine the effects of the uncertainties in the initial conditions and the model parameters on the model results. Marzban (2013) used this model to investigate the sensitivities of the model outputs versus parameters using the VSA approach. Yin et al. (2015) tested a nonlinear ensemble parameter perturbation approach for climate prediction based on this model. Furthermore, the Lorenz’63 model can simulate the atmospheric chaotic characteristic, which is qualitatively similar to that of the large-scale atmosphere circulation (Palmer 1999).
As is well known, in the Lorenz’63 model there are three state variables: the intensity of the convective motion x, the horizontal temperature variation y, and the vertical temperature variation z, as well as three additional parameters: the Prandtl number σ, the Rayleigh number r, and a nondimensional wavenumber b. The model is discretized using a fourth-order Runge–Kutta scheme with a dimensionless time step of 0.01. In this study, the reference parameter vector in the optimization problem (5) is set to P = (P1, P2. P3) = (σ, r, b) = (10, 28, 8/3). The parameter perturbation is denoted as p = (p1, p2, p3) = (σ′, r′, b′). The objective function norm is defined as the L2 norm; that is,
where the vector xt = () represents the perturbations of the model state variables at time t induced by parameter perturbations, which can be obtained through taking the difference between the model integration results at time t with and without parameter perturbations. Furthermore, the constraint condition for the parameter perturbations in (5) is defined as follows:
This constraint means that the amplitude of each parameter perturbation is located in a box of width ε times the absolute value of the corresponding reference parameter. So far, we have completed all settings of the optimization problem for the Lorenz’63 model. In the following, we will discuss the model parameter sensitivities by solving the optimization problem using the SPG algorithm. Prior to that, according to the OPSA method, the criterion value Cr for the Lorenz’63 model should be determined. Here, it is roughly defined as 30% of the maximal effects caused by uncertainties in all parameters [i.e., the last term in (6)]. This definition is simply because the maximal effects are induced by the uncertainties of three parameters in the model. We consider that on average each parameter uncertainty will approximately result in about 30% of the maximal effects. If the maximal effect caused by the uncertainty of some parameter is less than such an average level, this parameter is thought to be insensitive, relative to the other parameters. Of course, this definition is very rough, but by doing so the obtained results are consistent with those shown by Marzban (2013) in which parameter sensitivities were analyzed using the VSA approach (as mentioned below). To show the robustness of the results, different initial conditions, optimization (or prediction) times, and constraint conditions will be considered.
At first, the initial conditions are fixed at (x0, y0, z0) = (−14, −13, 47) which is the same as in Marzban (2013). The optimization time t is taken as 20 time steps. Different constraint conditions are taken. Table 1 shows the constraint conditions and the optimization results. The corresponding criterion values are also given in Table 1. It is seen that for different constraint values ε the objective function values for the parameter σ are always smallest and much less than those for the other two parameters. In addition, the objective functions for σ are less than the corresponding criterion values Cr, while the sums of the objective functions for σ and b are greater than Cr. This implies that the model outputs are insensitive to the variations in the parameter σ. This result is consistent with that obtained by Marzban (2013), in which the VSA approach was used.
Next the constraint condition is fixed at ε = 0.2, which represents that the amplitude of each parameter perturbation does not exceed 20% of the absolute value of the corresponding reference parameter value. The optimization times are taken as 20, 40, and 60 time steps. The results in Table 2 also show that the insensitive parameter is still the parameter σ and the result is robust for different optimization times. Furthermore, the objective function values induced by the perturbations of the parameters r and b are close and much greater than the criterion values, indicating that the model results are relatively sensitive to the r and b parameters.
To examine the influences of the model state on the parameter sensitivities under different initial conditions, the constraint condition and the optimization time are fixed at ε = 0.2 and t = 20 time steps, respectively. Table 3 shows the optimization results with respect to different initial conditions. It is found that the insensitive parameters vary for different initial conditions. Specifically, the objective function value for the parameter σ is the smallest for the initial conditions (−14, −13, and 47). However, for the initial conditions (0, 1.0, and 0) and (5.0, 5.0, and 5.0), the objective function values caused by the perturbations of the parameter b are the smallest. Simultaneously, the objective functions for b are less than the values Cr, while the sums of the objective functions for σ and b are greater than Cr. Hence, the sensitivities of the model parameters depend on the initial conditions. For the initial conditions (−14, −13, and 47), the insensitive parameter is σ, while the insensitive parameter is b for both sets of initial conditions [(0, 1.0, and 0) and (5.0, 5.0, and 5.0)].
The first inequality in (6) illustrates that the influences of parameter uncertainties are larger when the interactions among the parameters are taken into account. To quantify the differences, we obtain the optimization results of a single parameter by solving the maximization problem defined at the first line of (6). The other settings are the same as those when solving the optimization problem (5). Table 4 lists the results of the single-parameter optimization for different initial conditions. Comparing Tables 3 and 4, it is clear that the objective function values obtained from the single-parameter optimization are significantly less than those obtained when considering the parameter interaction. This demonstrates the importance of the interaction role between the model parameters.
In summary, the OPSA method was applied to explore the relative sensitivities of the short-term simulations of the Lorenz’63 model to different parameters based on three different sets of initial conditions. For the initial conditions (−14, −13, and 47), the model results are insensitive to the uncertainty in the parameter σ. Interestingly, we also find that the sensitivities of the Lorenz’63 model parameters are dependent on the initial state. For the initial states (0, 1.0, and 0) and (5.0, 5.0, and 5.0), the model results are almost unaffected by the uncertainty in the parameter b. Of course, we should note that the above results are dependent on the choice of Cr.
b. Intermediate-complexity shallow-water ocean model
In this section, we will further test the OPSA method in an intermediate-complexity ocean model: a 2.5-layer shallow-water model that is used to simulate the Kuroshio path variations south of Japan.
Before introducing the model, we first review the Kuroshio path variations. Previous study suggested that the Kuroshio path south of Japan exhibits an interannual variation between a typical large meander state and a non-large meander state (Taft 1972). In addition, Xu et al. (2010) and Nakamura et al. (2012) pointed out that the occurrence of the Kuroshio large meander (KLM) path has significant impacts on the local weather and climate. Hence, prediction of the KLM has attracted the attention of many researchers. Over the past decade, the studies examining KLM prediction mainly focus on the effects of the uncertainties in the initial conditions on the prediction (e.g., Kamachi et al. 2004; Ishikawa et al. 2004; Fujii et al. 2008; Wang et al. 2013). There are hardly any studies on the influences of the uncertainties among the model parameters. So, it is of significance to investigate the sensitivities of the model parameters and their effects in the KLM prediction.
Considering that shallow-water models can well reproduce the interannual variation of the Kuroshio path south of Japan (Qiu and Miao 2000; Ishikawa et al. 2004; Pierini 2006; Wang et al. 2012), we use a 2.5-layer shallow-water model to examine the sensitivity of the KLM prediction to model parameters applying the OPSA method. It is well known that the shallow-water model contains six parameters. To simulate the Kuroshio path variations, the values of these six parameters are, respectively, set to the two reduced gravity values, g1 = 0.029 and g2 = 0.048 m s−2, the lateral friction coefficient AH = 600 m2 s−1, the wind stress coefficient α = 0.1 Pa, the interfacial friction coefficient between two active layers RI = 1.17 × 10−7 s−1, and the bottom friction coefficient RB = 1.17 × 10−7 s−1. These values are often estimated from observations and, therefore, have relatively large uncertainties. As such, in this study we mainly focus on the relative sensitivities of these six parameters.
The model is discretized using the implicit Crank–Nicholson scheme, which allows us to take a relatively large time step (Dijkstra 2005). Here, a time step of 10 days is used. To simulate the Kuroshio path variations south of Japan, the model domain is set to a part of the western North Pacific basin V: (15°–55°N, 122°E–158°W) with a horizontal resolution of 0.2° × 0.2°. To keep the modeled Kuroshio from entering the East China Sea, the 200-m-deep contour is taken as the continental boundary. No-slip boundary conditions are used in the model. The monthly climatological wind stress of Hellerman and Rosenstein (1983) is used to drive the model.
The model is initialized from a state of rest and is integrated for 15 yr. After integrating for 9 yr, the solution of the model reaches a relatively steady state. So, the first 9 yr are regarded as the spinup period. We take a KLM process from the last 6 yr to show the modeled Kuroshio path variations. Figure 1 indicates that the Kuroshio takes a nonlarge meander path at 11.92 yr; after 40 days, a small meander path emerges around 33°N, 138°E; as this small meander path evolves, a typical KLM forms at 12.17 yr and the southernmost position of the Kuroshio axis is defined as the 480-m contour of the upper-layer thickness that is located at less than 31.5°N; subsequently, the KLM path is maintained for about 8 months (Figs. 1c–e) and the Kuroshio path eventually returns the nonlarge meander state at 12.94 yr. This process is similar to those obtained from other ocean models (Qiu and Miao 2000; Tsujino et al. 2006).
For simplicity, we will test the OPSA method for identifying the model parameter sensitivities only for the prediction of the KLM formation process. To solve the optimization problem (5), two KLM formation processes are chosen as the reference states, denoted as events 1 and 2. Event 1 is shown in Figs. 1a–c. The KLM formation process in event 2 is similar to that in event 1 (figure not shown) and the starting time is 13 yr. For these two events, the KLM formation periods are both 90 days (e.g., Figs. 1a–c). Hence, the optimization time in (5) is set to t = 90 days. Furthermore, the reference parameter vector in (5) is P = (P1, P2, P3, P4, P5, P6) = (g1, g2, AH, α, RI, RB). The parameter perturbation vector is denoted as p = (p1, p2, p3, p4, p5, p6) = (). The objective function norm is defined as the kinetic energy of the perturbation at the final time of optimization over the region Λ: 25°–35°N, 132°–140°E, where the KLM occurs. As such, the norm is written as
where xt = () is the perturbation vector of the model variables at time t caused by the model parameter perturbations; and represent the perturbations of the upper- and lower-layer zonal velocities, respectively; and are the perturbations of the meridional velocities; and are the upper- and lower-layer thickness perturbations; and and denote the upper- and lower-layer thicknesses, respectively, at time t in the reference state. This kinetic energy norm of the perturbation is defined because it is able to distinguish different Kuroshio path states. That is to say, the value of this norm is large when the difference between the Kuroshio path obtained after superimposing parameter perturbations and the reference Kuroshio path is large, and vice versa (Wang et al. 2012). We also define another norm such as the quadratic sum of the sea surface height (SSH) perturbation over the region Λ, but the results obtained were qualitatively similar to those acquired by the above norm. Hence, we only show the results based on the perturbation kinetic energy norm [see (11)] here.
In addition, similar to the constraint condition (10) for the Lorenz’63 model, here the amplitude of each parameter perturbation is also constrained within a box having a width ε times the absolute value of the corresponding reference parameter. After completing these settings, we solve the optimization problem (5) using the SPG algorithm.
The optimization experiments are performed for ε = 0.2. The reason for choosing this constraint condition is as follows: for the wind stress amplitude α in the North Pacific Ocean, Kutsuwada (1982) estimated its uncertainty using observational data and found that it is not beyond 0.02 Pa, which is equivalent to 20% of the reference value; for consistency, other parameters are also constrained within 20% of the corresponding reference values. In addition, and similar to section 3a, the criterion value is also defined as 30% of the maximal effects caused by uncertainties in all model parameters. Hence, the criterion values Cr for events 1 and 2 are 1.46 × 1012 and 1.39 × 1012 m5 s–2, respectively, which approximately corresponds to the SSH error of ~13 cm in the southern region of Japan. Interestingly, Usui et al. (2006) suggested, based on a data assimilation analysis, that the prediction for the Kuroshio path is unacceptable when the SSH error is beyond 13.5 cm. This is broadly consistent with the definition of Cr. That is to say, if the SSH error caused by uncertainties in some parameters is less than 13 cm, these parameters are considered to be insensitive because the prediction for the Kuroshio path is still acceptable even when these parameters contain uncertainties.
Figure 2 indicates the objective function values obtained by optimizing different parameters. For both events 1 and 2, the three smallest objective function values are caused by uncertainties in the bottom friction coefficient RB, the interfacial friction coefficient RI, and the lateral friction coefficient AH; their sums are 1.27 and 1.07 × 1012 m5 s−2, which are less than the corresponding criterion values, Cr. This implies that the KLM prediction is insensitive to the parameters RB, RI, and AH. The objective function values for the reduced gravity values g1 and g2, as well as the wind stress coefficient α, are relatively large and larger than the corresponding values Cr, especially for g1 and α. Physically, the reduced gravities represent the ocean stratification. Hence, the above results reflect that the changes in the ocean stratification and the amplitude of wind stress may have important effects on the KLM prediction.
We also perform the optimization calculation of a single parameter [see the first line of (6)] and the results are shown in Fig. 3. We find that the objective function values obtained from the single-parameter optimization are much less than those obtained by considering the parameter interaction role (cf. ordinates of Figs. 2 and 3). This implies that the interaction role between the parameter uncertainties is very strong in the KLM prediction. Similar to Fig. 2, Fig. 3 also indicates that the maximal effects of the uncertainties in the parameters RB and RI on the KLM prediction are very weak for both events 1 and 2. Additionally, if we pay attention to the rank of the parameter sensitivities, as in Zaehle et al. (2005) and Sun and Mu (2016), we find that the ranks are different in Figs. 3a and 3b. In Fig. 3a, the wind stress coefficient α is more sensitive than the parameter g2, while the situation is reversed in Fig. 3b. Hence, for the single-parameter optimization, the parameter sensitivities vary as reference states change, although the reference states are very similar and both correspond to the formation process of the KLM. However, for the OPSA method considering parameter interaction, the ranks of the parameter sensitivities for different reference states are almost consistent (see Fig. 2). Another interesting phenomenon is that for event 2, the objective function value for g2 is slightly less than that for α in Fig. 2b, while the former is significantly larger than the latter in Fig. 3b. This stresses the importance of considering the interaction among parameter uncertainties. In fact, Kurogi and Akitomo (2006) have demonstrated the importance of the interaction between the ocean stratification and the wind stress amplitude to the shift of the Kuroshio path.
4. Summary and discussion
In this study, a new optimization strategy called OPSA is proposed to determine the sensitivities of parameters in atmospheric and oceanic models. This strategy applies a nonlinear optimization method for estimating the maximum values of model parameter sensitivity measures. Relative to other methods, this approach is cheaper to compute. On the other hand, the OPSA method can consider the role of interaction among the parameter uncertainties.
The method is tested using the Lorenz’63 model and an intermediate-complexity 2.5-layer shallow-water model. For the Lorenz’63 model, it is found that the parameter sensitivities are unaffected by the different constraint conditions and optimization times, but they depend on the initial conditions (i.e., the reference state). For the initial conditions (−14, −13, and 47), the model output is insensitive to the parameter σ, while it is insensitive to the parameter b for both (0, 1.0, and 0) and (5.0, 5.0, and 5.0).
For the 2.5-layer shallow-water model used to simulate the KLM event south of Japan, the OPSA method reveals that the prediction of the KLM is insensitive to the uncertainties in the parameters RB, RI, and AH. On the contrary, the prediction is mostly affected by the uncertainties of the reduced gravity values g1 and g2 and the wind stress coefficient α. Interestingly, Kurogi and Akitomo (2006) and Akitomo (2008) also found, based on dynamics analysis, that the changes in the stratification in the upper ocean and the wind stress play important roles in the shift of the Kuroshio bimodal path. In addition, Sun et al. (2013) suggested that the variability of the ocean stratification can change the kinetic energy transfer from the atmosphere into the ocean’s upper layer induced by wind stress. For example, increasing stratification can reduce the kinetic energy transport from the ocean’s upper layer to the lower layer, which will induce an increase in the kinetic energy in the ocean’s upper layer. They further pointed out that the kinetic energy variations will induce changes in the flow velocity of the western boundary current. Hence, the uncertainties in the ocean stratification and the wind stress amplitude may have marked influences on the KLM prediction. Of course, further investigation is needed to clarify how the stratification and the wind stress amplitude affect the prediction, but this is beyond the scope of this paper.
The optimization of a single parameter without considering the interaction among parameter uncertainties is also calculated in this study. Sun and Mu (2016) confirmed that this single-parameter optimization is similar to the OAT method. The results here show that the objective function values obtained by the single-parameter optimization are much less than those obtained by our optimization strategy for both the Lorenz’63 model and the shallow-water model, which is consistent with the inequality (6). This implies that the uncertainties in the model outputs caused by parameter uncertainties are seriously underestimated by the single-parameter optimization method (or the OAT method). Hence, it is important to consider the interaction among parameter uncertainties when identifying model parameter sensitivities.
In the present study, the criterion value Cr is simply defined as 30% of the maximal effects caused by uncertainties in all model parameters. Of course, as mentioned in section 2, for the physical problems considered in the future study, the value Cr should be determined according to our understanding of the specific physical phenomena and processes. In addition, this study mainly uses two relatively simple models to examine the usefulness of the OPSA method. However, for a realistic atmospheric or oceanic general circulation model, the number of degrees of freedom is very high. In this situation, although the method is computationally cheaper than for other existing methods, it is still expensive to comprehensively investigate the sensitivities of model results to the uncertainties of the parameters. This holds especially for the dependence of the parameter sensitivity on the system state, because the model integrations will require considerable time and computational resources. On the other hand, for a complex model, the objective function defined in (5) may be nonsmooth. This causes great challenges to the realistic optimization calculation. Fortunately, some nonsmooth optimization algorithms have been developed that may be of use (e.g., Steward et al. 2012; Zheng et al. 2012). However, how to analyze the parameter sensitivity of a realistic atmospheric or oceanic model using the OPSA method based on the nonsmooth optimization algorithms still needs to be investigated, which is our next step of work.
The authors thank two anonymous reviewers for their useful suggestions. This study was supported by the National Natural Scientific Foundation of China (41306023, 41576015 and 41490644), the Qingdao National Laboratory for Marine Science and Technology (QNLM2016ORP0107), the NSFC Innovative Group (Grant 41421005), the National Programme on Global Change and Air-Sea interaction (GASI-IPOVAI-06), the NSFC-Shandong Joint Fund for Marine Science Research Centers (U1406401), and the Natural Sciences and Engineering Research Council (NSERC) of Canada Discovery Grant to YT. It was also partially supported by the Netherlands Organization for Scientific Research (NWO) under the Complexity project PreKurs.
Proof of the Inequality (8)
In (7), the objective function Jij can be written as
According to the expression of Jj, for the term in the second set of parentheses we have
The above inequality holds for arbitrary (p1, p2, …, pm) ∈ Cε. This results in the inequality (8):