Institute of Oceanology, Chinese Academy of Sciences
Article Information
 YUAN Shijin, LI Mi, WANG Qiang, ZHANG Kun, ZHANG Huazhen, MU Bin
 Optimal precursors of doublegyre regime transitions with an adjointfree method
 Journal of Oceanology and Limnology, 37(4): 11371153
 http://dx.doi.org/10.1007/s0034301972359
Article History
 Received Dec. 17, 2017
 accepted in principle Feb. 12, 2018
 accepted for publication Aug. 29, 2018
2 Key Laboratory of Ocean Circulation and Waves, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China
The doublegyre circulation includes a subpolar gyre and a subtropical gyre, which is a typical occurrence in the northern midlatitude ocean basins (Shen et al., 1999). Doublegyre circulation can capture the dynamic behaviors of the western boundary currents (WBC) system and reveal dynamic instability. The study for doublegyre circulation contributes to the study of the nonlinear dynamics of winddriven ocean circulations (Jiang et al., 1995; Dijkstra, 2005).
Moore (1999) examined the stochastic wind forcing of ocean gyre circulation using the generalized linear stability theory applied to the barotropic vorticity equation. Sura and Penland (Sura and Penland, 2002) showed how "minor" details of stochastic forcing could have a significant effect on the temporal and spatial response of the forced system. Primeau (Primeau, 2002; Primeau and Newman, 2008) used a quasigeostrophic model and a shallowwater model to study the impact of different viscosity coefficients and windforcing parameters on the multiple equilibria states and lowfrequency variability of doublegyre circulation. Nauw et al. (2004) indicated the effects of lateral friction on the asymmetries of doublegyre circulation. Pierini (2010) showed that a red noise wind forcing could induce lowfrequency variability in a shallowwater model through a coherence resonance mechanism. Sapsis and Dijkstra (2013) studied the interactions of additive wind stress noise and nonlinear dynamics in a quasigeostrophic model of doublegyre circulation.
Doublegyre regime transitions are one of lowfrequency variability phenomenon (Nauw and Dijkstra, 2001) of doublegyre circulation, the study on which can aid in understanding the dynamic mechanism and enhance the predictability of doublegyre circulation (Mahadevan et al., 2001).
Initial perturbations have an important effect on predictability (Mu et al., 2002). One focus of initial perturbations is to address the nature of them that is most likely to develop into atmospheric or climate events we are interested in. The initial perturbations under some conditions, are called as the optimal precursors (OPRs) (Mu et al., 2014) to the atmospheric and oceanic event. Finding the OPRs that cause doublegyre regime transitions will help us understand the dynamic mechanism and the process, improve the prediction skills.
CNOP method proposed by Mu et al. (2003) has been successfully applied to the predictability study of atmospheric and oceanic events induced by initial perturbations, such as El NiñoSouthern Oscillation (ENSO) (Duan et al., 2004; Mu et al., 2004, 2014), targeted observations of typhoons (Qin and Mu, 2011; Qin et al., 2013), Indian Ocean dipole (IOD) (Feng and Meyers, 2003), the Kuroshio Large Meander (Wang et al., 2012, 2013a, b) and the Kuroshio transport variation (Wang et al., 2017; Zhang et al, 2017a). The initial perturbations corresponding to CNOP leads to the largest nonlinear evolution at the prediction time under a given constraint. Among these researches, Terwisscha et al. (2008), Zhang et al. (2016) investigated the roles of initial perturbations on doublegyre circulation with a quasigeostrophic model and a 1.5layer shallowwater model respectively. ROMS can simulate doublegyre circulation and apparently considers the effect of baroclinic, in which there are barotropic and baroclinic modes (Shchepetkin and Mcwilliams, 2005). Studying CNOPtype initial perturbations with ROMS is a promising work. In this paper, we find the optimal precursors (OPRs) with ROMS, a kind of CNOPtype initial perturbations, which cause doublegyre regime transitions. We observe the evolution of the OPRs with time, superpose the OPRs to initial state to verify the doublegyre regime transitions. We analysis the dynamic mechanism on how barotropy and barocline have effects on doublegyre regime transitions from the point of energetics.
In addition, we propose a disjointfree method, called Principal Component Analysisbased Simulated Annealing (PCASA) algorithm, to solve CNOPtype initial perturbations of ROMS. PCASA is a member of the family of dimension reduction intelligent algorithms which have been applied successfully to solve CNOP of ZC model (Mu et al., 2015; Wen et al., 2015), MM5 model (Zhang et al., 2017b).
The remainder of this paper is organized as follows. The ROMS and CNOP method are briefly introduced in Section 2, and the solving CNOP of ROMS with PCASA algorithm is described in Section 3. Experiment settings, including ROMS initialization, cost function setting, equilibrium state acquisition, PCASA parameter setting, are introduced in Section 4. In Section 5, we show the experiment results and give the corresponding analyses. At last, the conclusions are given in Section 6.
2 ROMS AND CNOP METHOD 2.1 ROMSROMS can run standalone or be coupled to atmospheric and/or wave models. It follows the Earth System Modeling Framework (ESMF) (Hill et al., 2004) conventions for model coupling: initialize, run and finalize. The dynamical kernel of ROMS is comprised of four separate models including the nonlinear (NLM), tangent linear (TLM), representer tangent linear (RPM), and adjoint (ADM).
ROMS is a freesurface, terrainfollowing, primitive equations ocean model widely used by the scientific community for a diverse range of applications (e.g., Haidvogel et al., 2000; Di Lorenzo, 2003; Dinniman et al., 2003; Marchesiello et al., 2003; Peliz et al., 2003; Budgell, 2005; Warner et al., 2005a, 2005b; Wilkin et al., 2005). The algorithms that comprise ROMS computational nonlinear kernel are described in detail in Shchepetkin and McWilliams (2003, 2005), and the tangent linear and adjoint kernels and platforms are described in Moore et al. (2004). ROMS includes accurate and efficient physical and numerical algorithms and several coupled models for biogeochemical, biooptical, sediment, and sea ice applications. The sea ice model is described in Budgell (2005). It also includes several vertical mixing schemes (Warner et al., 2005a), multiple levels of nesting and composed grids.
ROMS can be configured for any region of the world's oceans, ranging from the local to the basin scale (Moore et al., 2004).
2.2 CNOP methodThe definition of CNOP is the initial perturbation x_{0}^{*} leading to the largest nonlinear evolution at prediction time t under a given constraint δ, which can be written as follows:
where
where M is the nonlinear model propagator from time 0 to t, x is the initial state, x_{0} is the perturbation, x_{t} is the evolution of x_{0} from time 0 to t.
3 SOLVING CNOP OF ROMS WITH PCASA ALGORITHMPCASA algorithm proposed by us is a member of the family of dimension reduction intelligent algorithms. It is an adjointfree method, so solving CNOP of ROMS with PCASA algorithm calls only the nonlinear (NLM) of ROMS, adjoint (ADM) is not called. The process for solving CNOP of ROMS with PCASA algorithm is depicted in Fig. 1, which contains two steps:
STEP 1: dimension reduction with Principal Component Analysis (PCA).
STEP 2: finding CNOPtype perturbations with Simulated Annealing (SA).
The input of STEP 1 is the data samples, the output of that is feature vectors of the data samples. This step is only done once unless the data samples are changed. The detailed process of this step is as follows:
(1) Obtaining data samples of perturbation by integrating ROMS.
(2) Preprocessing data samples from (1), including normalizing and centering.
(3) Carrying out eigenvalue decomposition of the covariance matrix of dataset from (2).
(4) Extracting feature vectors.
The input of STEP 2 is a solution generated randomly, the output of that is the optimal solution. Of course, for converging quickly of Simulated Annealing, we use some tricks to generate the above solution from ROMS outputs. This step is usually performed many times considering the stochasticity of Simulated Annealing. Among these optimal solutions, the optimal solutions corresponding the maximum of cost function is treated as CNOPtype perturbations. CNOPtype perturbations that cause doublegyre regime transitions are the optimal precursors (OPRs).
The detailed process of this step is as follows:
(1) Integrating ROMS after projecting the solution to original high dimensions space by feature vectors of the data samples from STEP 1.
(2) Computing cost function from (1).
(3) Updating the solution by metropolis criterion.
(4) Judging if the iteration termination condition is satisfied, if YES to (5), otherwise to (1).
(5) Output the solution, this is the optimal solution.
The Pseudocode for solving CNOP of ROMS with PCASA algorithm is in Table 1.
4 EXPERIMENT SETTING 4.1 ROMS SettingROMS is configured for simulating doublegyre circulation, whose setting follows Moore et al. (2004). The model domain is configured as a flatbottom rectangular ocean basin with dimensions of 1 000 km zonally, 2 000 km meridionally, and 500 m vertically, divided into four vertical levels with a thickness of 125 m. The model is solved on a midlatitude βplan centered at 45°N with a horizontal resolution of 18.5 km. The model is forced by a zonally uniform wind stress, which can be described as follows:
where τ_{0}= 0.05 N/m^{2}, y is meridional coordinate and L_{y} is the meridional extent of the basin.
Weakly relaxing the temperature of the water column to a simple analytical profile which is a function of both latitude and depth, which can be described as follows:
where a=0.95℃, T_{0}=5℃, y and z are meridional coordinate and vertical coordinate.
In addition, thirdorder upstream horizontal advection, secondorder horizontal diffusion along zsurfaces, linear bottom drag, and secondorder vertical viscosity and diffusion are used. Some other parameters are shown in Table 2.
4.2 Equilibrium states acquisitionWe focus on the doublegyre circulation transitions in the multipleequilibria regime. Viscosity coefficient is one of the important model parameters of ROMS, which is denoted as V_{c}. With V_{c}=580 m^{2}/s, we get nonperiod oscillations of doublegyre circulation states, the kinetic energy time sequence diagram for 10– 30 years is shown in Fig. 2, regarding the first 10 years as the model initialing phase. Wave peaks are jet up states and jet down states, wave troughs are symmetric states.
Wave peaks for t=4 790 days and t=5 040 days are selected, which are denoted by the red crosses in Fig. 2. From contour plots of sea surface height, it is a jet down state at t=4 790 days (similar to Fig. 4b), it is a jet up state at t=5 040 days (similar to Fig. 4a).
Taking the jet down state for t=4 790 days as the initial state, integrating ROMS for 30 years with different V_{c}, we find that only under the condition 860 m^{2}/s≤V_{c}≤1 140 m^{2}/s, doublegyre circulation keeps at the jet down state; 1 140 m^{2}/s≤V_{c}≤1 280 m^{2}/s, doublegyre circulation changes from the jet down state to the symmetric state. Taking the jet up state for t=5 040 days as the initial state, changing V_{c}, integrating ROMS for 30 years, we also find that only under the condition 860 m^{2}/s≤V_{c}≤1 140 m^{2}/s, doublegyre circulation keeps at the jet up state; 1 140 m^{2}/s≤V_{c}≤1 280 m^{2}/s, doublegyre circulation changes from the jet up state to the symmetric state. So for 860 m^{2}/s≤V_{c}≤1 140 m^{2}/s, multiequilibrium states can be acquired, whereas for 1 140 m^{2}/s≤V_{c}≤1 280 m^{2}/s, only one single equilibrium state is obtained.
To characterize the multiequilibrium states, on the basis of Simonnet et al. (2003), we define a modified transport difference (TD) variable, which is written as follows:
where ζ_{max} and ζ_{min} are the maximum and minimum heights of the sea surface, respectively.
From the definition, TD=0 indicates a perfect symmetric state, with two counterrotating gyres of equal strength and spatial extent. TD>0 corresponds to the state in which the subpolar gyre is more intense than the subtropical gyre, but is smaller in spatial extent than the latter. Conversely, TD < 0 refers to the state in which the subtropical gyre has a smaller extent, but is more intense than the subpolar gyre. Figure 3 gives the bifurcation diagram of these data, which shows the TDs of the equilibrium states with respect to V_{c}. The bifurcation diagram shows a perfect pitchfork structure, which agrees with the previous studies of Speich et al. (1995) and Terwisscha et al.(2008).
According to the shown bifurcation diagram, V_{c} =860 m^{2}/s is used to evaluate the oscillation between multiequilibrium states of the doublegyre circulation. The contour plots of sea surface height (SSH) for the multiequilibrium states in this circumstance are illustrated in Fig. 4. Figure 4a depicts the jetup equilibrium state, and Fig. 4b depicts the jetdown state, and these two states can coexist with V_{c} =860 m^{2}/s. Based on the multiequilibrium states, we intend to find the initial perturbations that can induce the transitions between them and investigate the transition process to better understand the dynamic mechanism of doublegyre circulation. These initial perturbations is viewed as the optimal precursors that might trigger the transitions The CNOP method is used to find these initial perturbations.
4.3 Cost function and constraint settingIn the study of the optimal precursors of doublegyre circulation regime transitions, we mainly focus on u (the zonal velocity), v (the meridional velocity) and ζ (SSH), the state vector incorporates u, v and ζ. So, x_{0} and x_{t} in Eq.7 can be written as follows:
where u_{t}, v_{t} and ζ_{t} denote the evolutions of initial perturbations u_{0}, v_{0} and ζ_{0} at time t, respectively.
Previous research shows that energy norm can distinguish among different conditions of doublegyre circulation (Jiang et al., 1995; Pierini, 2006). The most intensive effects of the jet and vortex occur in area Ʌ (0 km≤x≤600 km, 750 km≤y≤1250 km, marked by the white box in Fig. 3). Therefore, the objective function J(x_{0}) is defined as the energy norm of x_{t} over area Ʌ. Thus, J(x_{0}) is written as follows:
where h is the thickness of each vertical layer, and g is the gravitational acceleration.
The constraint is chosen as the total energy norm of initial perturbations in the entire domain, which can be written as follows:
In the following experiments, δ is set to 4.0×10^{11} m^{5}/ s^{2} for the jetup state as the reference state to study the transition from the jetup state to the jetdown state, 4.5×10^{11} m^{5}/s^{2} for the jetdown state as the reference state to study the transition from the jetdown state to the jetup state. The value of constraint δ is nearly 10%–12% of the energy norm of the corresponding reference states. Other constraints are also experimented. The CNOPtype perturbations for smaller constraints cannot cause the transition between doublegyre states at the specified prediction time, and that for larger constraints can cause the transition between doublegyre states. When constraints are greater than or equal to the above given constraints, patterns of the CNOPtype perturbations are similar for different constraints. The larger the constraint is, the faster doublegyre regime transitions occurs. So, the above given constraint is the minimum for which CNOPtype perturbations can cause the transition between doublegyre states.
In the experiment, the optimized time t is set to 1 000 days. We have also made experiments on CNOPtype perturbations obtained under other optimized times. It is found that the spatial patterns of the CNOPtype perturbations and their evolution are similar to those of 1 000 days when the optimization time is more than 1 000 days, and the doublegyre circulation will remain in the state of 1 000 days after 1 000 days. On the contrary, when the optimization time is less than 1 000 days, only the precursor of the transition can be observed, but the transition cannot be completed within the optimization time. However, due to the limitation of computational efficiency, when computing cost function, ROMS nonlinear model is only integrated for 100 days. The experimental results show that CNOPtype initial perturbations obtained by the objective function, which is the deviation degree between the evolution of perturbations for 100 days and the reference state, can also lead to the transition of doublegyre circulations.
Consequently, only the CNOPtype perturbation for δ=4.0×10^{11} m^{5}/s^{2} (when taking the jetup state as the reference state) and δ=4.0×10^{11} m^{5}/s^{2} (when taking the jetdown state as the reference state) with t=1 000 days are investigated in this paper. According to the definition of CNOP, the CNOPtype perturbation has the largest nonlinear evolution at prediction time t under the given constraint, and lead to the largest prediction error. So, the CNOPtype perturbation can be regarded as the optimal precursor of doublegyre regime transition.
4.4 PCASA settingThe model domain is configured as a flatbottom rectangular ocean basin with dimensions of 1 000 km zonally, 2 000 km meridionally, and 500 m vertically, divided into four vertical levels with a thickness of 125 m, the resolution of doublegyre simulation is 18.5 km. According to the horizontal grid setting and the vertical grid setting in Yu (1998), the dimension of u, v, ζ is 55×110×4, 56×109×4, 56×110 respectively, totally 54 776.
Before finding CNOPtype perturbations with SA, we reduce 54 776 dimensions with PCA to a low dimension. We obtain data samples of perturbation by integrating ROMS. We integrate the ROMS to simulate the evolution of the doublegyre circulation for 30 years with two different viscosity coefficients, 640 and 320 (units of m^{2}/s). The reason that we use these two coefficients is that the doublegyre circulation is period doubling with V_{c}=640 m^{2}/s and chaotic with V_{c}=320 m^{2}/s. Thus, we can obtain various perturbation samples by subtraction between the two outputs. During the 30 years integration, we save the corresponding records every 10 days, and thus both datasets both have sizes of 1 096. Hence, the subtraction operation between the two datasets can supply 1 096 samples of perturbations for the step of dimension reduction with PCA. The output of this step is feature vectors of the data samples.
In Eq.7, x_{0} is the initial perturbations with 54 776 dimensions, the aim for the step of dimension reduction with PCA is to represent x_{0} with linear combinations of feature vectors of the data samples, i.e. representing x_{0} with L·ω. Of course, there is the same constraint, i.e. L·ω≤δ.
L is the feature matrix from PCA, ω is the projection of x_{0} in feature space. The dimension of ω is far less than that of x_{0}.
We need to determine the number of eigenvectors in L, denoted as l. If a larger value is chosen, the feature space can absorb more information from the original space, but additional computation time is required. In contrast, a smaller l can speed up the computation, but a portion of information from the original space will be lost. To ascertain a proper value for l, we introduce an accumulative eigenvalue ratio R_{l} which is the percentage of energy conserved by the feature space, written as follows:
where λ_{i} is the i_{th} largest eigenvalues from PCA.
As shown in Fig. 5, the solid line represents the accumulative eigenvalue ratio, and the dashed line denotes its gradient. A proper value of l should make the accumulative eigenvalue ratio comparatively high but allow for a relatively low gradient. At l=70, the accumulative eigenvalue ratio is 81.62%, and the increase rate of the ratio is slow after l=70, according to the gradient. In other words, l greater than 70 results in limited accuracy but increases the computing time. Therefore, l is finally set to 70, and the dimension ω of is 70.
As discussed above, x_{0} with 54 776 dimensions is now represented with linear combinations of 70 feature vectors of the data samples, then SA randomly search the optimal initial perturbations in feature space with 70 dimensions. This accelerate the search process of SA.
The parameters setting of SA are shown in Table 3, which are achieved by repeated experiments. Generally, T_{0} is 10% of order of magnitude of energy norm of reference states.
5 EXPERIMENT RESULTS AND ANALYSESWe study both the transition from the jetup state to the jetdown state and the converse transition. PCASA is used to calculate the CNOPtype perturbations and the perturbations located on the boundary of the initial constraints, i.e., δ=4.0×10^{11} m^{5}/s^{2} for the jetup state as the reference state and δ=4.5×10^{11} m^{5}/s^{2} for the jetdown state as the reference state, which agrees with the conclusion of Liu (2008).
5.1 From the jetup state to the jetdown stateFor the transition from the jetup state to the jetdown state, the SSH patterns of the CNOPtype perturbation and its evolution are shown in Fig. 6. The perturbation first shows a twocell structure located in the area (0 km≤X≤600 km, 750 km≤Y≤1250 km). An anticyclonic eddy appears below this structure, and therefore the perturbation evolves into a threecell structure.
To investigate whether this type of perturbation can induce the transition between doublegyre states, the perturbation is superimposed on the jetup reference state, and the nonlinear model is integrated for 1 000 days. Figure 7 shows the SSH evolution of doublegyre circulation after superimposing the CNOPtype perturbation. It is noted that the doublegyre circulation presents a jetdown state at t=1 000 days, which means that the CNOPtype perturbation definitely induces the transition between doublegyre states.
Of further interest, we investigate the evolution of the CNOPtype perturbation and the transition process between doublegyre states to figure out the dynamic mechanism of the variability of doublegyre circulation. After superimposing the CNOPtype perturbation, the SSH pattern is clearly still in a jetup state at the initial time, and the perturbation mostly shows a twocell structure with one anticyclonic eddy on the top of a larger cyclonic eddy (Fig. 6a), and this type of structure makes the jetup extension elongated towards the north (Fig. 7b). Subsequently, the cyclonic eddy of the perturbation gradually extends towards the north (Fig. 6c–d), which causes the jetup extension to weaken and shed (Fig. 7c–d). At t=50 days, a strong anticyclonic eddy occurs below the cyclonic eddy, while the strength of the upper anticyclonic eddy gradually decreases (Fig. 6e–f). This pattern is similar to the socalled onebump gyre mode computed by Simonnet et al. (2003). Primeau and Newman (2008) also found a similar pattern during bifurcation analysis and revealed that such a type of gyre mode has a possible relationship to doublegyre regime transitions. This type of spatial structure makes the doublegyre circulation reach a symmetry state (Fig. 7f). The cyclonic eddy at the top gradually becomes stronger (Fig. 6g–i), and the doublegyre circulation finally reaches a jetdown state at t=1 000 days (Fig. 7i).
We divide this process into two stages. In the first stage from 0 days to 50 days, the perturbation primarily shows a twocell structure. With the extension of the cyclonic eddy towards the north, the jetup extension is therefore elongated and finally shed, which causes the doublegyre circulation to reach a symmetry state. The second stage occurs from 50 days to 1 000 days, during which an anticyclonic eddy appears below the cyclonic eddy, and the perturbation evolves into a threecell structure. The growth of the upper anticyclonic eddy drives the doublegyre circulation to finally reach a jetdown state.
5.2 From the jetdown state to the jetup stateSimilar results can be found for the transition from the jetdown state to the jetup state. Figure 8 shows the CNOPtype perturbation and its evolution, and Fig. 9 demonstrates the SSH evolution of doublegyre circulation after superimposing the CNOPtype perturbation. It can be observed that the variation occurs due to the CNOPtype perturbation and that the process is symmetrical to the converse transition discussed above.
5.3 The dynamic mechanism analysis from the point of energeticsWe can conclude that the transition between doublegyre states can be induced by the rapid growth of CNOP. Therefore, it is important to determine the contributors to the CNOP evolution. We explore the evolution of CNOPs with reference to the energetics analysis method (Chang and Oey, 2014; Oey, 2008). The energy evolution equation can be written as follows:
where u and v are the velocities in the zonal and meridional directions, respectively; g is the gravitational constant; ρ is the density; and N^{2} is the buoyancy frequency. Variables with primes represent the perturbation, and those without primes denote the reference state. The left terms in the bracket represent the total energy (including kinetic and potential energy) of the perturbation. The first and second terms on the righthand side are the barotropic conversation rate term (BT) and the baroclinic conversation rate term (BC), respectively. BT represents the energy conversation rate from the background kinetic energy to the perturbation kinetic energy, which can describe the barotropic instability, and BC represents the energy conversation rate from the background potential energy to the perturbation potential energy, which indicates the baroclinic instability. Figure 10 shows the potential vorticity of the jetup reference state at different layers. The value of potential vorticity at the first layer is larger than that of the third layer. Figure 11 shows the BTs and BCs for the transition from the jetup state to the jetdown state at t=50 and 1 000 days, respectively. It can be observed that larger BTs and BCs exist in the region where the CNOPtype perturbation is located, and this phenomenon means that CNOP obtains energy from the background state through both the barotropic and baroclinic instabilities. Moreover, compared with BC, BT is relatively larger, which implies that barotropic instability contributes more significantly to the rapid growth of CNOP.
Figures 12 and 13 show how the converse case transit from the jetdown state to the jetup state. The same conclusion can be obtained that BTs are always larger than BCs, which means that barotropic instability plays a more important role than baroclinic instability in the rapid growth of CNOP.
5.4 Random perturbationsTo verify whether CNOPtype perturbations are the optimal precursors, we generate 30 random initial perturbations for the jetup and jetdown reference state respectively, then superimpose these perturbations on the corresponding reference states, and the transitions don't occur after integrating ROMS. This means that CNOPtype perturbations are the optimal precursors of doublegyre regime transitions.
6 CONCLUSIONIn this paper, we propose PCASA algorithm to solve CNOP of ROMS aiming at studying the transitions between doublegyre states induced by CNOPtype perturbations.
We first explore the multiple linear steady states. The bifurcation diagram is found to have a pitchfork structure. To examine the nonlinear instability of doublegyre circulations, the jetup and jetdown states are separately taken as the reference state, and the CNOP approach is applied to find the optimal initial perturbations that might induce the transitions between these two states. Instead of using the adjointbased method to solve the CNOP, we apply the PCASA algorithm, which primarily consists of two steps: dimension reduction and simulated annealing. The dimension reduction process reduces the dimension of the solution space from 54 776 to 70, which makes it possible for simulated annealing to search for the global optimum.
To explore whether the calculated CNOPtype perturbations can trigger the transitions between doublegyre states, we superimpose the perturbations on the corresponding reference states, and the transitions occur after integrating ROMS. Taking the transition from the jetup state to the jetdown state as an example, the evolution process of CNOP can be divided into two stages. In the first stage, the perturbation primarily exhibits a twocell structure with one anticyclonic eddy on top of a larger cyclonic eddy. The northern extension of the cyclonic eddy makes the jetup extension elongate and shed, which causes the doublegyre circulation reach a symmetry state. In the second stage, an anticyclonic eddy appears below the cyclonic eddy, and the perturbation evolves into a threecell structure. The growth of the upper anticyclonic eddy drives the doublegyre circulation to finally reach a jetdown state. With the energetics analysis method, we find that the CNOP obtains energy from the background state through both the barotropic and baroclinic instabilities, and barotropic instability contributes more significantly to the rapid growth of CNOP. Similar results can be found in the converse process of transition from the jetdown state to the jetup state.
Random initial perturbations cannot cause the transitions of doublegyre circulations. CNOPtype perturbations are the optimal precursors of doublegyre regime transitions.
In the future, we will apply the CNOP approach to investigate the parameter errors of ROMS when simulating doublegyre circulation. After obtaining the optimal parameter errors combination of ROMS, we construct ensemble forecasting system obtaining samples by reducing parameter errors to improve the prediction skills of doublegyre circulation.
7 DATA AVAILABILITY STATEMENTThe data supporting the conclusions of this study are included in this article.
Budgell W P. 2005. Numerical simulation of iceocean variability in the Barents Sea region. Ocean Dynamics, 55(34): 370387.
DOI:10.1007/s1023600500083 
Chang Y L, Oey L Y. 2014. Analysis of STCC eddies using the OkuboWeiss parameter on model and satellite data. Ocean Dynamics, 64(2): 259271.
DOI:10.1007/s1023601306807 
Di Lorenzo E. 2003. Seasonal dynamics of the surface circulation in the Southern California current system. Deep Sea Research Part Ⅱ:Topical Studies in Oceanography, 50(1416): 23712388.
DOI:10.1016/S09670645(03)001255 
Dijkstra H A. 2005. Nonlinear Physical Oceanography:A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Nino. Springer Science & Business Media, New York.

Dinniman M S, Klinck J M, Jr W O S. 2003. Crossshelf exchange in a model of the Ross Sea circulation and biogeochemistry. DeepSea Research Part Ⅱ, 50(22): 31033120.

Duan W S, Mu M, Wang B. 2004. Conditional nonlinear optimal perturbations as the optimal precursors for El NinoSouthern Oscillation events. Journal of Geophysical Research:Atmospheres, 109(D23): D23105.

Feng M, Meyers G. 2003. Interannual variability in the tropical Indian Ocean:a twoyear timescale of Indian Ocean Dipole. Deep Sea Research Part Ⅱ:Topical Studies in Oceanography, 50(1213): 22632284.
DOI:10.1016/S09670645(03)000560 
Haidvogel D B, Arango H G, Hedstrom K, Beckmann A, MalanotteRizzoli P, Shchepetkin A F. 2000. Model evaluation experiments in the North Atlantic Basin:simulations in nonlinear terrainfollowing coordinates. Dynamics of Atmospheres and Oceans, 32(34): 239281.
DOI:10.1016/S03770265(00)00049X 
Hill C, DeLuca C, Balaji, Suarez M, Da Silva A. 2004. The architecture of the Earth System modeling framework. Computing in Science & Engineering, 6(1): 1828.

Jiang S, Jin F F, Ghil M. 1995. Multiple equilibria, periodic, and aperiodic solutions in a winddriven, doublegyre, shallowwater model. Journal of Physical Oceanography, 25(5): 764786.
DOI:10.1175/15200485(1995)025<0764:MEPAAS>2.0.CO;2 
Liu Y M. 2008. Maximum principle of conditional optimal nonlinear perturbation. JournalEast China Normal University (Natural Science), (2): 131134.
(in Chinese with English abstract) 
Mahadevan A, Lu J, Meacham S, MalanotteRizzoli P. 2001. The predictability of largescale winddriven flows. Nonlinear Processes in Geophysics, 8(6): 449465.
DOI:10.5194/npg84492001 
Marchesiello P, Mcwilliams J C, Shchepetkin A. 2003. Equilibrium structure and dynamics of the California current system. Journal of Physical Oceanography, 33(4): 753783.
DOI:10.1175/15200485(2003)33<753:ESADOT>2.0.CO;2 
Moore A M, Arango H G, Di Lorenzo E, Di Cornuelle B, Miller A J, Neilson D J. 2004. A comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model. Ocean Modelling, 7(12): 227258.
DOI:10.1016/j.ocemod.2003.11.001 
Moore A M. 1999. Windinduced variability of ocean gyres. Dynamics of Atmospheres and Oceans, 29(24): 335364.
DOI:10.1016/S03770265(99)00010X 
Mu B, Wen S C, Yuan S J, Li H Y. 2015. PPSO:PCA based particle swarm optimization for solving conditional nonlinear optimal perturbation. Computers & Geosciences, 83: 6571.

Mu M, Duan W S, Wang B. 2003. Conditional nonlinear optimal perturbation and its applications. Nonlinear Processes in Geophysics, 10(6): 493501.
DOI:10.5194/npg104932003 
Mu M, Duan W S, Wang J C. 2002. The predictability problems in numerical weather and climate prediction. Advances in Atmospheric Sciences, 19(2): 191204.
DOI:10.1007/s003760020016x 
Mu M, Sun L, Dijkstra H A. 2004. The sensitivity and stability of the ocean's thermohaline circulation to finiteamplitude perturbations. Journal of Physical Oceanography, 34(10): 2305.
DOI:10.1175/15200485(2004)034<2305:TSASOT>2.0.CO;2 
Mu M, Yu Y S, Xu H, Gong T T. 2014. Similarities between optimal precursors for ENSO events and optimally growing initial errors in El Niño predictions. Theoretical and Applied Climatology, 115(34): 461469.
DOI:10.1007/s007040130909x 
Nauw J J, Dijkstra H A, Chassignet E P. 2004. Frictionally induced asymmetries in winddriven flows. Journal of Physical Oceanography, 34(9): 20572072.
DOI:10.1175/15200485(2004)034<2057:FIAIWF>2.0.CO;2 
Nauw J J, Dijkstra H A. 2001. The origin of lowfrequency variability of doublegyre winddriven flows. Journal of Marine Research, 59(4): 567597.
DOI:10.1357/002224001762842190 
Oey L Y. 2008. Loop Current and deep eddies. Journal of Physical Oceanography, 38(7): 14261449.
DOI:10.1175/2007JPO3818.1 
Peliz A, Dubert J, Haidvogel D B, et al. 2003. Generation and unstable evolution of a densitydriven Eastern Poleward Current:The Iberian Poleward Current. Journal of Geophysical Research Oceans, 108(C8): 3268.
DOI:10.1029/2002JC001443 
Pierini S. 2006. A kuroshio extension system model study:decadal chaotic selfsustained oscillations. Journal of Physical Oceanography, 36(8): 16051625.
DOI:10.1175/JPO2931.1 
Pierini S. 2010. Coherence resonance in a doublegyre model of the Kuroshio Extension. Journal of Physical Oceanography, 40(1): 238248.
DOI:10.1175/2009JPO4229.1 
Primeau F, Newman D. 2008. Elongation and contraction of the western boundary current extension in a shallowwater model:a bifurcation analysis. Journal of Physical Oceanography, 38(7): 14691485.
DOI:10.1175/2007JPO3658.1 
Primeau F. 2002. Multiple equilibria and lowfrequency variability of the winddriven ocean circulation. Journal of Physical Oceanography, 32(8): 22362256.
DOI:10.1175/15200485(2002)032<2236:MEALFV>2.0.CO;2 
Qin X H, Duan W S, Mu M. 2013. Conditions under which CNOP sensitivity is valid for tropical cyclone adaptive observations. Quarterly Journal of the Royal Meteorological Society, 139(675): 15441554.
DOI:10.1002/qj.v139.675 
Qin X H, Mu M. 2011. A study on the reduction of forecast error variance by three adaptive observation approaches for tropical cyclone prediction. Monthly Weather Review, 139(7): 22182232.
DOI:10.1175/2010MWR3327.1 
Sapsis T P, Dijkstra H A. 2013. Interaction of additive noise and nonlinear dynamics in the doublegyre winddriven ocean circulation. Journal of Physical Oceanography, 43(2): 366381.
DOI:10.1175/JPOD12047.1 
Shchepetkin A F, Mcwilliams J C. 2003. A method for computing horizontal pressuregradient force in an oceanic model with a nonaligned vertical coordinate. Journal of Geophysical Research Oceans, 108(C3): 3090.
DOI:10.1029/2001JC001047 
Shchepetkin A F, Mcwilliams J C. 2005. The regional oceanic modeling system (ROMS):a splitexplicit, freesurface, topographyfollowingcoordinate oceanic model. Ocean Modelling, 9(4): 347404.
DOI:10.1016/j.ocemod.2004.08.002 
Shen J, Medjo T T, Wang S. 1999. On a winddriven, doublegyre, quasigeostrophic ocean model:numerical simulations and structural analysis. Journal of Computational Physics, 155(2): 387409.
DOI:10.1006/jcph.1999.6344 
Simonnet E, Ghil M, Ide K, Temam R, Wang S H. 2003. Lowfrequency variability in shallowwater models of the winddriven ocean circulation. Part Ⅱ:timedependent solutions. Journal of Physical Oceanography, 33(4): 729752.

Speich S, Dijkstra H, Ghil M. 1995. Successive bifurcations in a shallowwater model applied to the winddriven ocean circulation. Nonlinear Processes in Geophysics, 2(34): 241268.

Sura P, Penland C. 2002. Sensitivity of a doublegyre ocean model to details of stochastic forcing. Ocean Modelling, 4(3): 327345.

Van Scheltinga A D T, Dijkstra H A. 2008. Conditional nonlinear optimal perturbations of the doublegyre ocean circulation. Nonlinear Processes in Geophysics, 15(5): 727734.
DOI:10.5194/npg157272008 
Wang Q, Ma L B, Xu Q Q. 2013a. Optimal precursor of the transition from Kuroshio large meander to straight path. Chinese Journal of Oceanology and Limnology, 31(5): 11531161.
DOI:10.1007/s0034301323011 
Wang Q, Mu M, Dijkstra H A. 2012. Application of the Conditional Nonlinear Optimal Perturbation Method to the Predictability Study of the Kuroshio Large Meander. Advances in Atmospheric Sciences, 29(1): 118134.
DOI:10.1007/s0037601101990 
Wang Q, Mu M, Dijkstra H A. 2013b. The similarity between optimal precursor and optimally growing initial error in prediction of Kuroshio large meander and its application to targeted observation. Journal of Geophysical Research:Oceans, 118(2): 869884.
DOI:10.1002/jgrc.20084 
Wang Q, Tang Y M, Pierini S, Mu M. 2017. Effects of singularvectortype initial errors on the shortrange prediction of kuroshio extension transition processes. Journal of Climate, 30(15): 59615983.
DOI:10.1175/JCLID160305.1 
Warner J C, Geyer W R, Lerczak J A. 2005b. Numerical modeling of an estuary:A comprehensive skill assessment. Journal of Geophysical Research Oceans, 110(C5): C05001.

Warner J C, Sherwood C R, Arango H G, et al. 2005a. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Modelling, 8(1): 81113.

Wen S C, Yuan S J, Mu B, Li H Y, Ren J H. 2015. PCGD: Principal componentsbased great deluge method for solving CNOP. In: 2015 IEEE Congress on Evolutionary Computation (CEC). IEEE, Sendai, Japan.

Wilkin J L, Arango H G, Haidvogel D B, et al. 2005. A regional ocean modeling system for the Longterm Ecosystem Observatory. Journal of Geophysical Research Oceans, 110(C6): C06S91.

Yu J S. 1998. Setting up and calibration analysis of three dimensional ocean current forecasting mode (2/4). http://www.doc88.com/p704860833471.html. Accessed on 20121230. (in Chinese)

Zhang K, Wang Q, Mu M, Liang P. 2016. Effects of optimal initial errors on predicting the seasonal reduction of the upstream Kuroshio transport. Deep Sea Research Part Ⅰ:Oceanographic Research Papers, 116: 220235.
DOI:10.1016/j.dsr.2016.08.008 
Zhang L L, Yuan S J, Mu B, Zhou F F. 2017b. CNOPbased sensitive areas identification for tropical cyclone adaptive observations with PCAGA method. AsiaPacific Journal of Atmospheric Sciences, 53(1): 6373.
DOI:10.1007/s1314301700058 
Zhang X, Mu M U, Pierini S. 2017a. Optimal precursors triggering the kuroshio extension state transition obtained by the conditional nonlinear optimal perturbation approach. Advances in Atmospheric Sciences, 34(6): 685699.
DOI:10.1007/s0037601762637 