HTML
--> --> -->Previous studies have shown that the bimodal feature of the Kuroshio south of Japan has considerable effect on storm tracks, extratropical cyclone activity, and local precipitation (e.g., Xu et al., 2010; Nakamura et al, 2012; Hayasaki et al., 2013). Furthermore, it also has important influence on local fishery activity and ship navigation. Hence, it is very meaningful to be able to predict the Kuroshio path.
Recently, with greater understanding of the Kuroshio path variations and with the development of high-performance computers, predictability studies have been conducted based on the dynamics of Kuroshio path variations. Such work can be classified into two groups: estimation of the prediction time limit of LM events (e.g., Komori et al, 2003; Kamachi et al., 2004; Miyazawa et al, 2005; Usui et al., 2006), and examination of the effects of initial errors on predictions of the LM path. For example, (Fujii et al., 2008) used singular vector analysis to explore the effects of initial errors on the prediction of the LM. Their findings indicated that baroclinic instabilities resulted in further development of the LM path. (Wang et al., 2013) adopted the conditional nonlinear optimal perturbation (CNOP) method to overcome the linear approximation of the singular vector method in studying the effects of nonlinear physical processes on error growth in LM prediction. They found that the advection process has an important impact on error growth; however, they used a 1.5-layer reduced gravity model in which baroclinic instabilities were not considered. The above reflects the lack of consensus regarding the mechanism of nonlinear error growth and of its impact on LM prediction. Thus, it is interesting to investigate the evolution of nonlinear optimal initial errors during LM prediction using a multilayer model.
This study adopted the Regional Ocean Modeling System (ROMS), in conjunction with the CNOP approach proposed by (Mu et al, 2003), to explore the effects of optimal initial errors on the prediction of the LM path and to clarify the growth mechanism of the initial errors. The remainder of the paper is organized as follows: The model and methods adopted, including the CNOP calculation procedure, are described in section 2. Section 3 presents the derived results, including the patterns of the optimal initial errors, their evolution and effects on the prediction of the LM path, and the error growth mechanism. Finally, a summary and a discussion are presented in section 4.
2.1 Model simulation
ROMS, which is a state-of-the-art hydrostatic, free-surface, three-dimensional primitive equations ocean model, was adopted in this study to simulate the Kuroshio path variations. Considering the computational resource constraints, a one-way nested numerical simulation was used, where a 3/2 grid ratio was used to decrease the grid size from (1/8)°× (1/8)°cosφ (～14 km, where φ is latitude) in Nest 1 to (1/12)°× (1/12)°cosφ (～8 km) in Nest 2. As shown in Fig. 1, the domain of Nest 1 extended from 20°S to 60°N and from 100°E to 70°W, i.e., covering the entire North Pacific Ocean, while the domain of Nest 2 covered (23°-46°N, 122°-162°E), i.e., including the region within which the Kuroshio path variations occur. It should be noted that both the harmonic horizontal viscosity coefficient VIS2 and the diffusion coefficient TNU2 were set to 50.0 m2 s-1 in Nest 1 but to 40.0 m2 s-1 in Nest 2. Nest 1 was forced by monthly climatology data from the Comprehensive Ocean-Atmosphere Data Set (Diaz et al., 2002). The initial velocity and the free-surface elevation were set to zero, and the initial temperature and salinity values were prescribed by the monthly climatology of the World Ocean Atlas 2009 (available online at https://www.nodc.noaa.gov/OC5/WOA09/netcdf_data.html). All lateral boundaries were open and, while the boundary temperature and salinity values were also obtained from the World Ocean Atlas 2009, all other variables were derived by means of the geostrophic relation. The monthly averaged sea surface height (SSH), temperature, salinity, and zonal and meridional velocities were stored as output. Nest 1 was integrated for 50 years. After a spin-up of 30 years, the final 20 years' output of Nest 1 was used to drive Nest 2, which was integrated for 20 years under the same climatology forcing as Nest 1.
Figure1. Model domain and topography (units: m). The outer frame is Nest 1, while the inset black frame is Nest 2.It was necessary to examine the capability of ROMS in simulating the Kuroshio Current. The averaged SSHs obtained from the final 20 years' output and from satellite altimetry (AVISO, available online at https://www.aviso.altimetry.fr/en/data.html) are presented in Figs. 2a and b, respectively. ROMS appears successful in reproducing both the spatial pattern of the Kuroshio Current south of Japan and the Kuroshio Extension system, with a reasonable recirculation gyre located in the Shikoku Basin and two quasi-stationary meanders east of 140°E. Furthermore, we also validate the skill of the model to reproduce the seasonal variation of the Kuroshio path. Although some deviations from the observation existed in summer and winter, the model basically also captures the seasonal variation (figure not shown). To quantify the bimodal feature, the Kuroshio path index is introduced, which is defined as the latitude of the southernmost point of the Kuroshio axis within the longitudinal band from 135°E to 140°E. In this study, the Kuroshio axis refers to the 16-cm isoline of SSH, which corresponds approximately to the location of the maximal velocity of the Kuroshio Current. Following (Yoshida et al., 2006), we also regard 32°N as the threshold for identification of the Kuroshio LM path. Thus, when the index is smaller than 32°N, the Kuroshio south of Japan is considered to take the LM path; otherwise, it is considered to follow the NLM path. The time series of the Kuroshio path index, based on the monthly outputs of Nest 2, is shown in Fig. 3. It is found that the Kuroshio bimodal feature shows strong interannual variability and that four LM events occur during the 20 years. Although the simulated LM is slightly shorter and weaker than observed, which might be attributable to its high eddy activity (Akitomo, 2008), ROMS is able to capture the essential characteristics of the Kuroshio path variations under climatological forcing.
Figure2. Averaged SSH (units: m): (a) ROMS-simulated result for 20 years in Nest 1; (b) AVISO result averaged from 1993 to 2012.
Figure3. Kuroshio path index based on monthly outputs when horizontal coefficients VIS2= TNU2=40.0 (gray bands denote the LM path).2
2.2 CNOP approach and calculation procedure
In this study, the CNOP approach was used to search for the optimal initial errors affecting prediction of the LM path. Next, we briefly describe the method. The nonlinear model can be formally written as: \begin{equation} {X}_t={M}_t({X}_0) , \ \ (1) \end{equation} where X0 is the initial state vector, Xt is the state vector at time t, and Mt denotes the nonlinear propagator from the initial time to the future time t. When we superimpose an error x0 that satisfies the constraint $\|$x0$\|$A ≤δ (δ is the constraint radius) on the initial state X0, a nonlinear constraint optimization problem is defined as: \begin{equation} \label{eq2} J({x}_{0\delta})=\mathop{\max}\limits_{\|x_0\|_A\le\delta}J({x}_0)=\mathop{\max}\limits_{\|x_0\|_A\le\delta}\|{M}_t({X}_0+{x}_0)-{M}_t({X}_0)\|_B^2 , \ \ (2)\end{equation}where J(x0) represents the objective function that measures the nonlinear evolution of the initial error x0 at time t, $\| \cdot \|$A is the constraint norm, $\| \cdot \|$B is the objection function norm, and x0δ is the initial error that brings about the largest nonlinear evolution at time t, which is called the CNOP or the optimal initial error.
Here, we present the CNOP calculation procedure using the nonlinear optimization system based on the ROMS and the Spectral Projected Gradient the Spectral Projected Gradient algorithm (Birgin et al., 2000; Zhang et al., 2016):
(1) Selection of the initial constraint. In this study, we chose the sum of the quadratic errors nondimensionalized and normalized by the standard deviation as the initial constraint, similar to (Li et al, 2014). It contains all possible factors affecting the LM, which can be expressed as follows: \begin{equation} \|\delta {x}_0\|=\sqrt{\left(\frac{u'}{\bar{u}_{\rm std}}\right)^2+\left(\frac{v'}{\bar{v}_{\rm std}}\right)^2+ \left(\frac{T'}{\bar{T}_{\rm std}}\right)^2+\left(\frac{s'}{\bar{s}_{\rm std}}\right)^2+\left(\frac{h'}{\bar{h}_{\rm std}}\right)^2} , \ \ (3)\end{equation} where u'(v',T',s',h') and $\bar{u}_{std}(\bar{v}_{std},\bar{T}_{std},\bar{s}_{std},\bar{h}_{std})$ represent the error and the standard deviation of the SSH (zonal velocity, meridional velocity, temperature, and salinity), respectively, within the model domain (23°-46°N, 122°-162°E).
(2) Selection of the objective function. In this study, the objective function was defined as the kinetic energy over the area (30°-35°N, 132°-140°E), in which the LM path occurs. Greater kinetic energy means the amplitude of the LM is greater. The specific expression can be written as: \begin{eqnarray} \label{eq4} &&\|{M}_t({X}_0+{x}_0)-{M}_t({X}_0)\|_{{\rm KE}}=\|{x}(t)\|_{{\rm KE}}\nonumber\\ &&=\frac{\rho_{\rm ref}}{2}\int_{z=0}^{z=1000}\int_{y=25}^{y=35}\int_{x=135}^{x=140}[(u'_t)^2+(v'_t)^2]dxdydz ,\quad \ \ (4) \end{eqnarray} where u't and v't are the zonal and meridional velocity perturbations at time t, and ρ ref denotes the reference density.
(3) Selection of appropriate LM events in searching for CNOP-type errors. There are four LM events reproduced by the model. The first and second LM path starts to form in winter, while the third and fourth LM path starts to form in autumn. To avoid a large amount of calculation, we eventually chose the second and fourth LM events (denoted Case 1 and Case 2 in Fig. 3) as the background states to conduct the analysis in this paper.
(4) Choosing the appropriate optimization time. With consideration of the period of transition between the two types of path and the simulated results, we selected 70 days as the optimization time. Thus, Case 1 encompassed 1 December, Model Year 7 to 10 February, Model Year 8; and Case 2 encompassed 1 October, Model Year 17 to 10 December, Model Year 17. Taking Case 1 as an example, the formation process of the simulated LM path based on the SSH is shown in Fig. 4a. We can see that a small meander is generated to the east of Kyushu, and that it propagates downstream prior to the formation of the LM, which is similar to the LM formation based on the AVISO data (Fig. 4b). This also validates the ability of ROMS to simulate the mesoscale processes during the LM formation from another perspective.
Figure4. Formation process of the Kuroshio LM path: (a) Case 1; (b) AVISO. Shading represents the SSH field (units: m).(5) Selection of the constraint radius. Based on practical observation and the findings of previous studies, the constraint radius was set to δ=1× 107. Actually, when the change range of the radius is small, the patterns of the calculated optimal initial errors are similar to each other; therefore, only the results for δ =1× 107 are described here.
3.1. Structures and effects of CNOP-type initial errors
Following the optimization calculation, two CNOP-type initial errors (denoted as CNOP1 and CNOP2) were obtained for each case. To investigate their spatial patterns, Fig. 5 shows the distribution of the computed errors integrated vertically from the surface to the maximum depth in accordance with the initial constraint in Eq. (5). It can be seen that the areas of large amplitude of CNOP1 and CNOP2 for both cases are located mainly in the region upstream of the LM, which indicates the initial errors southeast of Kyushu around (28°-32°N, 130°-134°E) are most important in the prediction of the LM path. In addition, we also examined the vertical structures of the CNOP-type initial errors. It was found that the large amplitudes are restricted to areas in the upper 2500 m (especially 500-2000 m), indicating that the initial errors at these depths have the greatest effect on LM prediction.
Figure5. Spatial patterns of calculated CNOPs integrated vertically from the surface to the maximum depth according to the initial constraint in Eq. (5): CNOP1 (left); CNOP2 (right).To explore the relationship between CNOP1 and CNOP2, we introduce the similarity index S=$\langle$e_1,e_2$\rangle$/$\|$e_1$\|\|$e_2$\|$, where $\|$e_i$\|$=$\sqrt{(e_i,e_i)}$(i=1,2) is the Euclidean inner product over the area with the large amplitude CNOPs, and e1 and e2 represent CNOP1 and CNOP2, respectively. By calculation, the similarity coefficient between CNOP1 and CNOP2 was determined as -0.5943 (Case 1) and -0.5349 (Case 2), prompting consideration of the effect of the high negative correlation between CNOP1 and CNOP2 on LM prediction.
We superimposed CNOP1 and CNOP2 on the initial fields of the background states and integrated the nonlinear model for 70 days to explore their effects on LM prediction. Figure 6 presents the Kuroshio axes in different states for the two cases, where the black, red, and blue lines represents the background state, the field with CNOP1, and the field with CNOP2, respectively. It can be seen for both cases that CNOP1 causes the LM path to become strengthened and to move southwestward on Day 70, whereas CNOP2 causes the LM path to become weakened and to move northeastward. To quantify the effects of the CNOP-type errors, we considered the prediction error denoted by the objective function value as a measure, as in (Wang et al., 2012). Table 1 lists the objective function values for the CNOP-type initial errors in both cases, which indicate the prediction error caused by CNOP1 is larger than that caused by CNOP2. Furthermore, the objective function value for CNOP1 (CNOP2) in Case 1 is not identical to that for CNOP1 (CNOP2) in Case 2, although the difference is small. Overall, the two types of CNOP impose effects on the prediction of the LM path that act in two different directions.
2
3.2. Nonlinear evolution of CNOPs
To explore how the initial errors influence the LM path, we examined the nonlinear growth processes of the CNOPs to help us understand the physical mechanism of LM formation. The nonlinear evolutions of the SSH anomaly components for the final 40 days in the error fields of Cases 1 are shown in Fig. 7, where the scales of shading are different for different dates. It can be seen that the evolution pattern of the SSH anomaly component for CNOP1 is similar to that for CNOP2, except the sign is opposite. Specifically, for CNOP1 in Case 1 (left-hand panels of Fig. 7), the negative anomaly that exists around (32°N, 135°E) on Day 40 gradually moves downstream along the Kuroshio axis, while both its amplitude and its intensity increase. On Day 70, the negative anomaly is advected to the region where the Kuroshio path variations occur, and its further development tends to extend the LM path southwestward in the background field. In addition, it should be noted that while the negative anomaly moves and develops, the positive anomaly to its east part also develops rapidly. The spatial pattern of the negative and positive anomalies slows the downstream propagation of the LM, contributing to its further enlargement. For CNOP2 in Case 1 (right-hand panels of Fig. 7), almost the opposite evolution pattern of the SSH anomaly component can be seen, although there are some minor differences. The positive anomaly positioned to the east of the Kii Peninsula on Day 40 moves downstream along the Kuroshio axis and its amplitude increases slowly (Day 50 and Day 60). After Day 60, the positive anomaly grows rapidly, as presented on Day 70. This strong positive anomaly directly reduces the intensity of the negative anomaly corresponding to the LM path in the background field, causing the predicted LM path to become weakened. Furthermore, the negative anomaly downstream of the positive anomaly on Day 40 moves closer to the Japanese coast. To a certain extent, it offsets the primary positive anomaly and it hinders downstream propagation of the meander, causing the LM path to become weakened. Similar results were obtained for the calculated CNOPs of Case 2.
Figure6. Kuroshio axis on Day 70, with reference to the 16-cm isoline of SSH: (a) Case 1; (b) Case 2. The black, red, and blue lines represent the Kuroshio axis for the background state, the Kuroshio axis with CNOP1, and the Kuroshio axis with CNOP2, respectively.
Figure7. Nonlinear evolution of the SSH anomaly components in the error field for Case 1 (units: m; left panels: CNOP1; right panels: CNOP2).2
3.3. Growth mechanism of optimal initial errors
To explore the error growth mechanism, we employed the eddy-energetics analysis method (e.g., Tsujino et al., 2006; Fujii et al., 2008; Oey, 2008; Chang and Oey, 2014). The calculation of eddy energetics usually requires decomposition of an instantaneous field into the mean and the eddy field. In this study, the background state is exactly treated as the mean field and the error field is treated as the eddy field. Then, by calculating the time-tendency of the eddy kinetic energy (EKE) in the quasi-geostrophic framework, the barotropic conversion rate (BT) of the mean kinetic energy to the EKE, which is associated with barotropic instabilities, can be expressed as: \begin{equation} \frac{\partial({\rm EKE})}{\partial t}\!=\!\overbrace{-\rho_0\!\left(\!u'u'\frac{\partial\bar{u}}{\partial x}\!+\! u'v'\!\left(\!\frac{\partial\bar{u}}{\partial y}\!+\!\frac{\partial\bar{v}}{\partial x}\!\right)\!+\! v'v'\frac{\partial\bar{v}}{\partial y}\!\right)}^{{\rm BT}}+{\rm other\;terms\,}, \ \ (5)\end{equation} where $\bar{u}$ and $\bar{v}$ are the zonal and meridional velocities, respectively, in the mean field, while u' and v' are the corresponding velocity errors in the eddy field. It should be noted that the EKE is transferred to the mean kinetic energy when BT is negative. Similarly, by calculating the time-tendency of the eddy potential energy (EPE) in the quasi-geostrophic framework, the baroclinic conversion rate (BC) of the mean potential energy to the EPE, which is relevant to baroclinic instabilities, is defined as: \begin{equation} \frac{\partial({\rm EPE})}{\partial t}=\overbrace{\frac{g}{\rho_0\rho_z}\left(u'\rho'\frac{\partial\rho}{\partial x} +v'\rho'\frac{\partial\rho}{\partial y}\right)}^{{\rm BC}}+{\rm other\ terms }, \ \ (6)\end{equation} where g is gravitational acceleration and ρ0 is the reference density. Here, ρ (ρ') denotes the density in the mean (eddy) field. Similarly, the EPE is absorbed into the mean field when BC is negative.It needs to be mentioned that although the meanings of BC and BT are not exactly the same as those in the eddy-energetics analysis because of the different definitions of the background and eddy fields, they are still considered as the indicators of baroclinic instability and barotropic instability, as (Fujii et al., 2008) stated. Then, based on the above definitions, BT and BC for CNOP1 and CNOP2 in the two cases were computed, as shown in Figs. 8 and 9. The contour lines denote the SSH fields for the CNOP-type errors and the shading represents the spatial patterns of BT and BC at 400 m in Case 1. Because the propagation speed and the growth rate of the meander at the beginning of the development of the LM are both slow, the figures show only the results for the final 30 days. Obviously, BT and BC both work during the predicted LM formation and they are almost of the same order of magnitude. Furthermore, their spatial patterns are also similar, indicating that both barotropic and baroclinic instabilities are related to the growth of the optimal initial errors.
Figure8. SSH (contours) and energy conversion rates (shading) at 400 m during the predicted LM formation caused by CNOP1 for Case 1 (units: 10-7 m2 s-3).
Figure9. As in Fig. 8 but for CNOP2 in Case 1.For CNOP1 (Fig. 8), BT and BC are mostly positive on both flanks of the growing meander, where the errors acquire energy from the background field and demonstrate substantial evolution, implying flow instabilities have important effects on the LM formation process. This explains why CNOP1 causes the LM to move westward, while the LM path in the background field is strengthened, as shown in Fig. 6a. In fact, the positive-negative-positive patterns of the energy conversions correspond to the error evolution shown in the left-hand panels of Fig. 7. The positive values of BT and BC on both sides of the meander indicate that the background field provides the energy for the error evolution through barotropic and baroclinic instabilities. It seems that the error growth to the east slows the downstream propagation of the growing error to the west, thus extending the LM path southwestward. For CNOP2 in Case 1, although the evolution leads to a weakened LM path, similar conclusions regarding BT and BC can be derived. As shown in Fig. 9, BT and BC are both positive to the west of the meander but negative to the east, which means the error in the west develops more but the error evolution is suppressed (Fig. 7). Therefore, growth of the positive SSH anomaly in the west offsets the negative anomaly of the background LM path, causing the LM path to become weakened. For Case 2, the spatial patterns of BT and BC are similar to those of Case 1, i.e., the meander tends to move toward the region in which BT and BC are positive.
Irrespective of whether considering CNOP1 or CNOP2, both BT and BC show large positive values in similar regions, and barotropic and baroclinic instabilities play important roles in the prediction of the LM path. They are complementary to each other and none should be ignored.
The above analysis is dependent on our simulation; therefore, for further verification, we analyzed the error growth mechanism of the 2004/05 LM path based on JCOPE2 reanalysis data (Miyazawa et al, 2009). Similar to (Usui et al., 2008), the daily outputs were decomposed into values for the previous five days (mean field) and the corresponding deviation (eddy field). Then, BT and BC at 400 m during LM formation were computed, as shown in Fig. 10. It can be seen that not only are the distributions of BC and BT similar to each other, but also their magnitudes are the same. Especially on 20 July 2004, the positive values of BC and BT at the tip of the meander suggest the direction of meander growth in the following moment, as shown on 20 August. At that time, the positive values of energy conversion are located in the western region of the meander, indirectly explaining why the meander moves westward during further development. The analysis based on JCOPE2 data also suggests the importance of barotropic and baroclinic instabilities in the LM formation process, which is consistent with our results based on ROMS simulation.
Figure10. SSH (contours) and energy conversion rates (shading) at 400 m during the 2004/05 LM formation process based on JCOPE2 reanalysis data (units: 10-7 m2 s-3).Using selected LM events (Case 1 and Case 2) as background fields, two CNOP-type initial errors (CNOP1 and CNOP2) were obtained. The results showed that the large amplitudes of CNOPs were concentrated mainly in areas in the upper 2500 m in the upstream of the LM, i.e., to the southeast of Kyushu. The rapid nonlinear growth of CNOP1 was found to strengthen the LM through extension of the LM path southwestward, while the rapid nonlinear growth of CNOP2, whose evolution pattern is opposite to that of CNOP1, tended to weaken the LM through contraction of the LM path northeastward. Energy analysis of the error field indicated that BT and BC have similar magnitudes and both show positive values that indicate the direction of movement of the LM path. Thus, barotropic and baroclinic instabilities both play important roles in the error growth, although the specific effects of baroclinic instabilities were not analyzed in detail because they lie beyond the scope of this study. In addition, it needs to be mentioned that we also calculated the CNOP-type initial errors of the first and third LM events and obtained similar results (figure not shown).
This study constitutes the first attempt to explore the three-dimensional optimal initial errors of LM formation and to reveal the growth mechanism of the initial errors using the ROMS and CNOP method. Because of computational demands, the study focused on the results for just two cases. In future work, additional cases will be examined. In addition, the location of the obtained initial errors could contribute to the identification of sensitive areas for targeted observation for improved prediction of the LM path, which would also provide useful guidance for the deployment of observational instrumentation.
