Journal of Oceanology and Limnology   2019, Vol. 37 issue(4): 1137-1153     PDF       
http://dx.doi.org/10.1007/s00343-019-7235-9
Institute of Oceanology, Chinese Academy of Sciences
0

Article Information

YUAN Shijin, LI Mi, WANG Qiang, ZHANG Kun, ZHANG Huazhen, MU Bin
Optimal precursors of double-gyre regime transitions with an adjoint-free method
Journal of Oceanology and Limnology, 37(4): 1137-1153
http://dx.doi.org/10.1007/s00343-019-7235-9

Article History

Received Dec. 17, 2017
accepted in principle Feb. 12, 2018
accepted for publication Aug. 29, 2018
Optimal precursors of double-gyre regime transitions with an adjoint-free method
YUAN Shijin1, LI Mi1, WANG Qiang2, ZHANG Kun2, ZHANG Huazhen1, MU Bin1     
1 School of Software Engineering, Tongji University, Shanghai 201804, China;
2 Key Laboratory of Ocean Circulation and Waves, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China
Abstract: In this paper, we find the optimal precursors which can cause double-gyre regime transitions based on conditional nonlinear optimal perturbation (CNOP) method with Regional Ocean Modeling System (ROMS). Firstly, we simulate the multiple-equilibria regimes of double-gyre circulation under different viscosity coefficient and obtain the bifurcation diagram, then choose two equilibrium states (called jet-up state and jet-down state) as reference states respectively, propose Principal Component Analysis-based Simulated Annealing (PCASA) algorithm to solve CNOP-type initial perturbations which can induce double-gyre regime transitions between jet-up state and jet-down state. PCASA algorithm is an adjoint-free method which searches optimal solution randomly in the whole solution space. In addition, we investigate CNOP-type initial perturbations how to evolve with time. The results show:(1) the CNOP-type perturbations present a two-cell structure, and gradually evolves into a three-cell structure at predictive time; (2) by superimposing CNOP-type perturbations on the jet-up state and integrating ROMS, double-gyre circulation transfers from jet-up state to jet-down state, and vice versa, and random initial perturbations don't cause the transitions, which means CNOP-type perturbations are the optimal precursors of double-gyre regime transitions; (3) by analyzing the transition process of double-gyre regime transitions, we find that CNOP-type initial perturbations obtain energy from the background state through both barotropic and baroclinic instabilities, and barotropic instability contributes more significantly to the fast-growth of the perturbations. The optimal precursors and the dynamic mechanism of double-gyre regime transitions revealed in this paper have an important significance to enhance the predictability of double-gyre circulation.
Keywords: optimal precursors    double-gyre regime transitions    conditional nonlinear optimal perturbation (CNOP)    Principal Component Analysis-based Simulated Annealing (PCASA)    multipleequilibria regimes    
1 INTRODUCTION

The double-gyre circulation includes a sub-polar gyre and a sub-tropical gyre, which is a typical occurrence in the northern mid-latitude ocean basins (Shen et al., 1999). Double-gyre circulation can capture the dynamic behaviors of the western boundary currents (WBC) system and reveal dynamic instability. The study for double-gyre circulation contributes to the study of the nonlinear dynamics of wind-driven 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 quasi-geostrophic model and a shallow-water model to study the impact of different viscosity coefficients and wind-forcing parameters on the multiple equilibria states and low-frequency variability of double-gyre circulation. Nauw et al. (2004) indicated the effects of lateral friction on the asymmetries of double-gyre circulation. Pierini (2010) showed that a red noise wind forcing could induce low-frequency variability in a shallow-water model through a coherence resonance mechanism. Sapsis and Dijkstra (2013) studied the interactions of additive wind stress noise and nonlinear dynamics in a quasi-geostrophic model of double-gyre circulation.

Double-gyre regime transitions are one of lowfrequency variability phenomenon (Nauw and Dijkstra, 2001) of double-gyre 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 double-gyre 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ño-Southern 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 quasi-geostrophic model and a 1.5-layer shallow-water model respectively. ROMS can simulate double-gyre circulation and apparently considers the effect of baroclinic, in which there are barotropic and baroclinic modes (Shchepetkin and Mcwilliams, 2005). Studying CNOP-type initial perturbations with ROMS is a promising work. In this paper, we find the optimal precursors (OPRs) with ROMS, a kind of CNOP-type initial perturbations, which cause double-gyre regime transitions. We observe the evolution of the OPRs with time, superpose the OPRs to initial state to verify the double-gyre 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 disjoint-free method, called Principal Component Analysis-based Simulated Annealing (PCASA) algorithm, to solve CNOP-type 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 ROMS

ROMS 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 free-surface, terrain-following, 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, bio-optical, 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 method

The definition of CNOP is the initial perturbation x0* leading to the largest nonlinear evolution at prediction time t under a given constraint δ, which can be written as follows:

    (1)

where

    (2)

where M is the nonlinear model propagator from time 0 to t, x is the initial state, x0 is the perturbation, xt is the evolution of x0 from time 0 to t.

3 SOLVING CNOP OF ROMS WITH PCASA ALGORITHM

PCASA algorithm proposed by us is a member of the family of dimension reduction intelligent algorithms. It is an adjoint-free 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:

Fig.1 The process for solving CNOP of ROMS with PCASA algorithm

STEP 1: dimension reduction with Principal Component Analysis (PCA).

STEP 2: finding CNOP-type 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) Pre-processing 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 CNOP-type perturbations. CNOP-type perturbations that cause double-gyre 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 Pseudo-code for solving CNOP of ROMS with PCASA algorithm is in Table 1.

Table 1 The Pseudo-code for solving CNOP of ROMS with PCASA algorithm
4 EXPERIMENT SETTING 4.1 ROMS Setting

ROMS is configured for simulating double-gyre circulation, whose setting follows Moore et al. (2004). The model domain is configured as a flat-bottom 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 mid-latitude β-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:

    (3)

where τ0= 0.05 N/m2, y is meridional coordinate and Ly 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:

    (4)
    (5)

where a=0.95℃, T0=5℃, y and z are meridional coordinate and vertical coordinate.

In addition, third-order upstream horizontal advection, second-order horizontal diffusion along z-surfaces, linear bottom drag, and second-order vertical viscosity and diffusion are used. Some other parameters are shown in Table 2.

Table 2 Some parameter values
4.2 Equilibrium states acquisition

We focus on the double-gyre circulation transitions in the multiple-equilibria regime. Viscosity coefficient is one of the important model parameters of ROMS, which is denoted as Vc. With Vc=580 m2/s, we get nonperiod oscillations of double-gyre 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.

Fig.2 The kinetic energy time sequence diagram of doublegyre circulation for 10–30 years with Vc=580 m2/s Wave peaks for t=4 790 days (jet down state) and t=5 040 days (jet up state) are selected, which are denoted by the red crosses

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).

Fig.4 SSH for the multi-equilibrium states a. jet-up equilibrium state; b. jet-down equilibrium state. The white box indicates area Ʌ.

Taking the jet down state for t=4 790 days as the initial state, integrating ROMS for 30 years with different Vc, we find that only under the condition 860 m2/s≤Vc≤1 140 m2/s, double-gyre circulation keeps at the jet down state; 1 140 m2/s≤Vc≤1 280 m2/s, double-gyre 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 Vc, integrating ROMS for 30 years, we also find that only under the condition 860 m2/s≤Vc≤1 140 m2/s, double-gyre circulation keeps at the jet up state; 1 140 m2/s≤Vc≤1 280 m2/s, double-gyre circulation changes from the jet up state to the symmetric state. So for 860 m2/s≤Vc≤1 140 m2/s, multi-equilibrium states can be acquired, whereas for 1 140 m2/s≤Vc≤1 280 m2/s, only one single equilibrium state is obtained.

To characterize the multi-equilibrium states, on the basis of Simonnet et al. (2003), we define a modified transport difference (TD) variable, which is written as follows:

    (6)

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 counter-rotating gyres of equal strength and spatial extent. TD>0 corresponds to the state in which the sub-polar 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 sub-polar gyre. Figure 3 gives the bifurcation diagram of these data, which shows the TDs of the equilibrium states with respect to Vc. The bifurcation diagram shows a perfect pitchfork structure, which agrees with the previous studies of Speich et al. (1995) and Terwisscha et al.(2008).

Fig.3 Bifurcation diagram of double-gyre circulation

According to the shown bifurcation diagram, Vc =860 m2/s is used to evaluate the oscillation between multi-equilibrium states of the double-gyre circulation. The contour plots of sea surface height (SSH) for the multi-equilibrium states in this circumstance are illustrated in Fig. 4. Figure 4a depicts the jet-up equilibrium state, and Fig. 4b depicts the jetdown state, and these two states can co-exist with Vc =860 m2/s. Based on the multi-equilibrium 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 double-gyre 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 setting

In 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, x0 and xt in Eq.7 can be written as follows:

    (7)

where ut, vt and ζt denote the evolutions of initial perturbations u0, v0 and ζ0 at time t, respectively.

Previous research shows that energy norm can distinguish among different conditions of double-gyre 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(x0) is defined as the energy norm of xt over area Ʌ. Thus, J(x0) is written as follows:

    (8)

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:

    (9)

In the following experiments, δ is set to 4.0×1011 m5/ s2 for the jet-up state as the reference state to study the transition from the jet-up state to the jet-down state, 4.5×1011 m5/s2 for the jet-down state as the reference state to study the transition from the jet-down state to the jet-up state. The value of constraint δ is nearly 10%–12% of the energy norm of the corresponding reference states. Other constraints are also experimented. The CNOP-type perturbations for smaller constraints cannot cause the transition between double-gyre states at the specified prediction time, and that for larger constraints can cause the transition between double-gyre states. When constraints are greater than or equal to the above given constraints, patterns of the CNOP-type perturbations are similar for different constraints. The larger the constraint is, the faster double-gyre regime transitions occurs. So, the above given constraint is the minimum for which CNOP-type perturbations can cause the transition between double-gyre states.

In the experiment, the optimized time t is set to 1 000 days. We have also made experiments on CNOP-type perturbations obtained under other optimized times. It is found that the spatial patterns of the CNOP-type perturbations and their evolution are similar to those of 1 000 days when the optimization time is more than 1 000 days, and the double-gyre 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 CNOP-type 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 double-gyre circulations.

Consequently, only the CNOP-type perturbation for δ=4.0×1011 m5/s2 (when taking the jet-up state as the reference state) and δ=4.0×1011 m5/s2 (when taking the jet-down state as the reference state) with t=1 000 days are investigated in this paper. According to the definition of CNOP, the CNOP-type perturbation has the largest nonlinear evolution at prediction time t under the given constraint, and lead to the largest prediction error. So, the CNOP-type perturbation can be regarded as the optimal precursor of double-gyre regime transition.

4.4 PCASA setting

The model domain is configured as a flat-bottom 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 double-gyre 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 CNOP-type 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 double-gyre circulation for 30 years with two different viscosity coefficients, 640 and 320 (units of m2/s). The reason that we use these two coefficients is that the double-gyre circulation is period doubling with Vc=640 m2/s and chaotic with Vc=320 m2/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, x0 is the initial perturbations with 54 776 dimensions, the aim for the step of dimension reduction with PCA is to represent x0 with linear combinations of feature vectors of the data samples, i.e. representing x0 with L·ω. Of course, there is the same constraint, i.e. ||L·ω||≤δ.

L is the feature matrix from PCA, ω is the projection of x0 in feature space. The dimension of ω is far less than that of x0.

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 Rl which is the percentage of energy conserved by the feature space, written as follows:

    (10)

where λi is the ith 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.

Fig.5 Accumulative eigenvalue ratio and gradient The solid line represents the accumulative eigenvalue ratio, and the dashed line denotes its gradient.

As discussed above, x0 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, T0 is 10% of order of magnitude of energy norm of reference states.

Table 3 The parameters setting of SA
5 EXPERIMENT RESULTS AND ANALYSES

We study both the transition from the jet-up state to the jet-down state and the converse transition. PCASA is used to calculate the CNOP-type perturbations and the perturbations located on the boundary of the initial constraints, i.e., δ=4.0×1011 m5/s2 for the jet-up state as the reference state and δ=4.5×1011 m5/s2 for the jetdown state as the reference state, which agrees with the conclusion of Liu (2008).

5.1 From the jet-up state to the jet-down state

For the transition from the jet-up state to the jet-down state, the SSH patterns of the CNOP-type perturbation and its evolution are shown in Fig. 6. The perturbation first shows a two-cell structure located in the area (0 km≤X≤600 km, 750 km≤Y≤1250 km). An anti-cyclonic eddy appears below this structure, and therefore the perturbation evolves into a three-cell structure.

Fig.6 SSH evolution of the CNOP-type perturbation for the transition from the jet-up state to the jet-down state

To investigate whether this type of perturbation can induce the transition between double-gyre states, the perturbation is superimposed on the jet-up reference state, and the nonlinear model is integrated for 1 000 days. Figure 7 shows the SSH evolution of double-gyre circulation after superimposing the CNOP-type perturbation. It is noted that the doublegyre circulation presents a jet-down state at t=1 000 days, which means that the CNOP-type perturbation definitely induces the transition between double-gyre states.

Fig.7 SSH evolution of the double-gyre circulation after superimposing the CNOP-type perturbation for the transition from the jet-up state to the jet-down state

Of further interest, we investigate the evolution of the CNOP-type perturbation and the transition process between double-gyre states to figure out the dynamic mechanism of the variability of double-gyre circulation. After superimposing the CNOP-type perturbation, the SSH pattern is clearly still in a jet-up state at the initial time, and the perturbation mostly shows a two-cell structure with one anti-cyclonic eddy on the top of a larger cyclonic eddy (Fig. 6a), and this type of structure makes the jet-up extension elongated towards the north (Fig. 7b). Subsequently, the cyclonic eddy of the perturbation gradually extends towards the north (Fig. 6cd), which causes the jet-up extension to weaken and shed (Fig. 7cd). At t=50 days, a strong anti-cyclonic eddy occurs below the cyclonic eddy, while the strength of the upper anti-cyclonic eddy gradually decreases (Fig. 6ef). This pattern is similar to the so-called 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 double-gyre regime transitions. This type of spatial structure makes the double-gyre circulation reach a symmetry state (Fig. 7f). The cyclonic eddy at the top gradually becomes stronger (Fig. 6gi), and the double-gyre circulation finally reaches a jet-down 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 two-cell structure. With the extension of the cyclonic eddy towards the north, the jet-up extension is therefore elongated and finally shed, which causes the double-gyre circulation to reach a symmetry state. The second stage occurs from 50 days to 1 000 days, during which an anti-cyclonic eddy appears below the cyclonic eddy, and the perturbation evolves into a three-cell structure. The growth of the upper anti-cyclonic eddy drives the double-gyre circulation to finally reach a jet-down state.

5.2 From the jet-down state to the jet-up state

Similar results can be found for the transition from the jet-down state to the jet-up state. Figure 8 shows the CNOP-type perturbation and its evolution, and Fig. 9 demonstrates the SSH evolution of double-gyre circulation after superimposing the CNOP-type perturbation. It can be observed that the variation occurs due to the CNOP-type perturbation and that the process is symmetrical to the converse transition discussed above.

Fig.8 As in Fig. 6, but for the transition from the jet-down state to the jet-up state
Fig.9 As in Fig. 7, but for the transition from the jet-down state to the jet-up state
5.3 The dynamic mechanism analysis from the point of energetics

We can conclude that the transition between double-gyre 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:

    (11)

where u and v are the velocities in the zonal and meridional directions, respectively; g is the gravitational constant; ρ is the density; and N2 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 right-hand 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 jet-up 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 jet-up state to the jet-down 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.

Fig.10 Potential vorticity of the jet-up reference state at the first layer (a) and the third layer (b)
Fig.11 BTs at t=50 days (a) t=1 000 days (c) and BCs at t=50 days (b) t=1 000 days (d) for the transition from the jet-up state to the jet-down state

Figures 12 and 13 show how the converse case transit from the jet-down state to the jet-up 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.

Fig.12 Potential vorticity of the jet-down reference state at the first layer (a) and the third layer (b)
Fig.13 As in Fig. 11, but for the transition from the jet-down state to the jet-up state
5.4 Random perturbations

To verify whether CNOP-type perturbations are the optimal precursors, we generate 30 random initial perturbations for the jet-up and jet-down reference state respectively, then superimpose these perturbations on the corresponding reference states, and the transitions don't occur after integrating ROMS. This means that CNOP-type perturbations are the optimal precursors of double-gyre regime transitions.

6 CONCLUSION

In this paper, we propose PCASA algorithm to solve CNOP of ROMS aiming at studying the transitions between double-gyre states induced by CNOP-type 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 double-gyre circulations, the jet-up and jet-down 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 CNOP-type perturbations can trigger the transitions between double-gyre states, we superimpose the perturbations on the corresponding reference states, and the transitions occur after integrating ROMS. Taking the transition from the jet-up state to the jet-down state as an example, the evolution process of CNOP can be divided into two stages. In the first stage, the perturbation primarily exhibits a two-cell structure with one anti-cyclonic eddy on top of a larger cyclonic eddy. The northern extension of the cyclonic eddy makes the jet-up extension elongate and shed, which causes the double-gyre circulation reach a symmetry state. In the second stage, an anti-cyclonic eddy appears below the cyclonic eddy, and the perturbation evolves into a three-cell structure. The growth of the upper anti-cyclonic eddy drives the double-gyre circulation to finally reach a jet-down 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 jet-down state to the jet-up state.

Random initial perturbations cannot cause the transitions of double-gyre circulations. CNOP-type 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 double-gyre 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 double-gyre circulation.

7 DATA AVAILABILITY STATEMENT

The data supporting the conclusions of this study are included in this article.

References
Budgell W P. 2005. Numerical simulation of ice-ocean variability in the Barents Sea region. Ocean Dynamics, 55(3-4): 370-387. DOI:10.1007/s10236-005-0008-3
Chang Y L, Oey L Y. 2014. Analysis of STCC eddies using the Okubo-Weiss parameter on model and satellite data. Ocean Dynamics, 64(2): 259-271. DOI:10.1007/s10236-013-0680-7
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(14-16): 2371-2388. DOI:10.1016/S0967-0645(03)00125-5
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. Cross-shelf exchange in a model of the Ross Sea circulation and biogeochemistry. Deep-Sea Research Part Ⅱ, 50(22): 3103-3120.
Duan W S, Mu M, Wang B. 2004. Conditional nonlinear optimal perturbations as the optimal precursors for El Nino-Southern Oscillation events. Journal of Geophysical Research:Atmospheres, 109(D23): D23105.
Feng M, Meyers G. 2003. Interannual variability in the tropical Indian Ocean:a two-year time-scale of Indian Ocean Dipole. Deep Sea Research Part Ⅱ:Topical Studies in Oceanography, 50(12-13): 2263-2284. DOI:10.1016/S0967-0645(03)00056-0
Haidvogel D B, Arango H G, Hedstrom K, Beckmann A, Malanotte-Rizzoli P, Shchepetkin A F. 2000. Model evaluation experiments in the North Atlantic Basin:simulations in nonlinear terrain-following coordinates. Dynamics of Atmospheres and Oceans, 32(3-4): 239-281. DOI:10.1016/S0377-0265(00)00049-X
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): 18-28.
Jiang S, Jin F F, Ghil M. 1995. Multiple equilibria, periodic, and aperiodic solutions in a wind-driven, double-gyre, shallow-water model. Journal of Physical Oceanography, 25(5): 764-786. DOI:10.1175/1520-0485(1995)025<0764:MEPAAS>2.0.CO;2
Liu Y M. 2008. Maximum principle of conditional optimal nonlinear perturbation. Journal-East China Normal University (Natural Science), (2): 131-134. (in Chinese with English abstract)
Mahadevan A, Lu J, Meacham S, Malanotte-Rizzoli P. 2001. The predictability of large-scale wind-driven flows. Nonlinear Processes in Geophysics, 8(6): 449-465. DOI:10.5194/npg-8-449-2001
Marchesiello P, Mcwilliams J C, Shchepetkin A. 2003. Equilibrium structure and dynamics of the California current system. Journal of Physical Oceanography, 33(4): 753-783. DOI:10.1175/1520-0485(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(1-2): 227-258. DOI:10.1016/j.ocemod.2003.11.001
Moore A M. 1999. Wind-induced variability of ocean gyres. Dynamics of Atmospheres and Oceans, 29(2-4): 335-364. DOI:10.1016/S0377-0265(99)00010-X
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: 65-71.
Mu M, Duan W S, Wang B. 2003. Conditional nonlinear optimal perturbation and its applications. Nonlinear Processes in Geophysics, 10(6): 493-501. DOI:10.5194/npg-10-493-2003
Mu M, Duan W S, Wang J C. 2002. The predictability problems in numerical weather and climate prediction. Advances in Atmospheric Sciences, 19(2): 191-204. DOI:10.1007/s00376-002-0016-x
Mu M, Sun L, Dijkstra H A. 2004. The sensitivity and stability of the ocean's thermohaline circulation to finite-amplitude perturbations. Journal of Physical Oceanography, 34(10): 2305. DOI:10.1175/1520-0485(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(3-4): 461-469. DOI:10.1007/s00704-013-0909-x
Nauw J J, Dijkstra H A, Chassignet E P. 2004. Frictionally induced asymmetries in wind-driven flows. Journal of Physical Oceanography, 34(9): 2057-2072. DOI:10.1175/1520-0485(2004)034<2057:FIAIWF>2.0.CO;2
Nauw J J, Dijkstra H A. 2001. The origin of low-frequency variability of double-gyre wind-driven flows. Journal of Marine Research, 59(4): 567-597. DOI:10.1357/002224001762842190
Oey L Y. 2008. Loop Current and deep eddies. Journal of Physical Oceanography, 38(7): 1426-1449. DOI:10.1175/2007JPO3818.1
Peliz A, Dubert J, Haidvogel D B, et al. 2003. Generation and unstable evolution of a density-driven 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 self-sustained oscillations. Journal of Physical Oceanography, 36(8): 1605-1625. DOI:10.1175/JPO2931.1
Pierini S. 2010. Coherence resonance in a double-gyre model of the Kuroshio Extension. Journal of Physical Oceanography, 40(1): 238-248. 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): 1469-1485. DOI:10.1175/2007JPO3658.1
Primeau F. 2002. Multiple equilibria and low-frequency variability of the wind-driven ocean circulation. Journal of Physical Oceanography, 32(8): 2236-2256. DOI:10.1175/1520-0485(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): 1544-1554. 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): 2218-2232. DOI:10.1175/2010MWR3327.1
Sapsis T P, Dijkstra H A. 2013. Interaction of additive noise and nonlinear dynamics in the double-gyre wind-driven ocean circulation. Journal of Physical Oceanography, 43(2): 366-381. DOI:10.1175/JPO-D-12-047.1
Shchepetkin A F, Mcwilliams J C. 2003. A method for computing horizontal pressure-gradient 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 split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9(4): 347-404. DOI:10.1016/j.ocemod.2004.08.002
Shen J, Medjo T T, Wang S. 1999. On a wind-driven, doublegyre, quasi-geostrophic ocean model:numerical simulations and structural analysis. Journal of Computational Physics, 155(2): 387-409. DOI:10.1006/jcph.1999.6344
Simonnet E, Ghil M, Ide K, Temam R, Wang S H. 2003. Lowfrequency variability in shallow-water models of the wind-driven ocean circulation. Part Ⅱ:time-dependent solutions. Journal of Physical Oceanography, 33(4): 729-752.
Speich S, Dijkstra H, Ghil M. 1995. Successive bifurcations in a shallow-water model applied to the wind-driven ocean circulation. Nonlinear Processes in Geophysics, 2(3-4): 241-268.
Sura P, Penland C. 2002. Sensitivity of a double-gyre ocean model to details of stochastic forcing. Ocean Modelling, 4(3): 327-345.
Van Scheltinga A D T, Dijkstra H A. 2008. Conditional nonlinear optimal perturbations of the double-gyre ocean circulation. Nonlinear Processes in Geophysics, 15(5): 727-734. DOI:10.5194/npg-15-727-2008
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): 1153-1161. DOI:10.1007/s00343-013-2301-1
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): 118-134. DOI:10.1007/s00376-011-0199-0
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): 869-884. DOI:10.1002/jgrc.20084
Wang Q, Tang Y M, Pierini S, Mu M. 2017. Effects of singularvector-type initial errors on the short-range prediction of kuroshio extension transition processes. Journal of Climate, 30(15): 5961-5983. DOI:10.1175/JCLI-D-16-0305.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): 81-113.
Wen S C, Yuan S J, Mu B, Li H Y, Ren J H. 2015. PCGD: Principal components-based 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 Long-term 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/p-704860833471.html. Accessed on 2012-12-30. (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: 220-235. DOI:10.1016/j.dsr.2016.08.008
Zhang L L, Yuan S J, Mu B, Zhou F F. 2017b. CNOP-based sensitive areas identification for tropical cyclone adaptive observations with PCAGA method. Asia-Pacific Journal of Atmospheric Sciences, 53(1): 63-73. DOI:10.1007/s13143-017-0005-8
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): 685-699. DOI:10.1007/s00376-017-6263-7