Preferred equilibrium solutions of the barotropic vorticity equation

The long-term behaviour of the barotropic vorticity equation (BVE), which is beneficial for understanding long-term climate variations, has been investigated by applying an explicit quadratic conserved finite difference integration scheme. The results suggest that an initial pattern predominated at the beginning of time evolves to an equilibrium state after a long transition period in which the conserved prerequisites, such as the integral kinetic energy, square vorticity, and vorticity, are broken by accumulated errors. This could cause the BVE to be indeterministic. Similar phenomena could reappear in integration cases starting at the same initial pattern with different small perturbations, initial patterns, time and grid intervals, and difference schemes. The eventual equilibrium states for the above cases are different from each other but spatially similar. They all possess a single-wave structure in the horizontal plane but with random phases. The results further imply that despite the indeterministic characteristics due to the nonlinearity, the long-term behaviour of the BVE is deterministic in a spatial structure, shedding light on long-term climate change studies.


Introduction
The linearized barotropic vorticity equation (BVE) led the discovery of atmospheric Rossby-Haurwitz waves (Rossby 1939;Haurwitz 1940), which greatly motivated atmospheric and oceanic dynamic research. The nonlinear BVE can be simplified through the use of a double Fourier series and with the omission of all but the largest scales of motion (Lorenz 1960). The socalled low-order model (LOM) reduces the nonlinear BVE to an autonomous dynamical system with three ordinary nonlinear differential equations whose analytical solutions are elliptic functions of time, and then, the equilibrium solutions, stability, and trajectories can be further analysed (Lakshmivarahan et al. 2006). Due to the convenience of reducing the nonlinearity, LOMs have become a standard method and have been extensively applied in studies of complex atmospheric phenomena (Saltzman 1962;Lorenz 1963;Charney and DeVore 1979;Wiin-Nielsen 1992;Gluhovsky et al. 2002).
Numerical modelling has been greatly developed since the first successful numerical solution of the BVE over a limited area of the Earth's surface in 1950 (Charney et al. 1950). Many finite difference schemes with integral constraints on quadratic quantities have been designed to integrate the BVE and avoid integration instability (Lilly 1965;Arakawa 1966). These analogues have been widely adopted in many theoretical studies (Galewsky et al. 2004;Thuburn and Li 2000) and numerical general circulation models (Bryan 1963;Grimmer and Shaw 1967;Kurihara and Holloway 1967;Smagorinsky et al. 1965). Computational instability for long-term integration may occur if no artificial dissipation term is added (Grammeltvedt 1969). This is because the quadratic conservation constraints are broken when the leap-frog time difference scheme is used (Zeng and Ji 1981). The computational stability for long-term integration can be improved by designing schemes with complete quadratic conservation constraints (Ji and Zeng 1982;Ji and Wang 1991). These schemes offer powerful tools for exploring the long-term behaviour of the nonlinear BVE, which, to our knowledge, has not yet been studied, unlike its well-known short-term evolution. Therefore, the long-term integration behaviour of the BVE under periodic boundary conditions has been explored in this paper by using a finite difference scheme. The different schemes used for calculations are described in Section 2. The results are presented in Section 3, and the conclusions are summarized in Section 4.

Finite difference integration schemes
The dimensionless BVE for two-dimensional incompressible flow (Lorenz 1960) can be written as where t is time; u and v are horizontal velocity vector components in the x direction and y direction, respectively; and ς ¼ ∂v ∂x − ∂u ∂y is the vertical vorticity. By introducing the stream function ψ, the horizontal velocity and the vertical vorticity can be expressed as where ∇ 2 is the horizontal Laplacian operator. Following Lorenz (1960), we neglect the latitudinal variation of the rotating effect of the Earth, which has vast applications in atmospheric and oceanic motions. The case for the Coriolis parameter is discussed separately. The quadratic conservation features of Eq. (1) can be derived by setting double periodic condition boundaries Þis the kinetic energy (KE) and S is the integration domain. As mentioned above, a finite difference scheme should satisfy such prerequisites to directly solve the long-term numerical solution of Eq. (1). The quadratic conservations in a finite difference scheme are written as where superscript n and n + 1 denote the time indexes; the notation A k k 2 ¼ d 2 ∑  Applying the incompressible condition, Eq. (1) can be rewritten as The direct discretization of Eq. (5) can produce the implicit finite difference scheme that is adopted in this paper are notations for conveniently expressing the scheme; the subscripts i, j are the grid indexes in the x and y directions; and dt, d are the time and grid intervals, respectively. Equation (6) is a stable quadratic conservation scheme (Ji and Zeng 1982); however, the computation time is large, especially for long-term integration. Therefore, explicit quadratic conservation difference schemes are designed to speed up the computing time. An explicit finite difference scheme can be constructed by modifying Eq. (6) as follows: where ε n dtBς n is a term that is introduced to ensure that the scheme satisfies quadratic conservation, ε n is an undetermined coefficient, and Bis a dissipation operator. Equation (7) can be shortened as where A denotes the difference operator in Eq. (7). Ji and Wang (1991) indicated that Eq. (8) would be a quadratic conservation scheme by taking where K 1 = ‖Aς n ‖ 2 /(Bς n , ς n ), K 2 = (Bς n , Aς n ) ⋅ d/(Bς n , ς n ),K 3 = ‖Aς n ‖ ⋅ ‖Bς n ‖ ⋅ d/(Bς n , ς n ), and e ς nþ1 ¼ ς n −dtAς n . The notation (F, G) denotes the inner product of the two variables F and G in the domain S. Its expression is defined as F⋅GΔS. Refer to Ji and Zeng (1982) and Ji and Wang (1991) for detailed information. , and vorticity (c) curves for the whole integration period. To glorify the evolution features, the whole integration period is divided into three parts that are bound by the vertical lines in each figure Lorenz (1960) expanded ψ by a double Fourier series and retained three terms to conduct a theoretical analysis: where A, F, G are intensity coefficients and k, l are the wavenumbers in the x and y directions, respectively. It is necessary to retain at least three terms to theoretically analyse the nonlinear interaction among different terms (Fjörtoft 1953). Therefore, the integration process is started by specifying A = 0.12, F = 0.24, G = 0, k = 2, and l = 1 (see Fig. 1a) and integrating for a long time span (t = 10 5 , approximately 34 years for the unit time scale is 3 h) with grid and time intervals of 2π/40 (approximately 157 km for the unit length scale is 1000 km) and 0.005, respectively. In contrast, numerical solutions do not always reproduce the variations in the intensity of the average strength of the mid-latitude westerly winds, a phenomenon known as the index cycle (Namias 1950). At the very beginning, the stream function pattern evolves into an index cycle with circulation cells appearing and displaced alternately ( Fig. 1b-g). However, the index cycle does not last an indefinitely long time as the analytic solution predicted. After the beginning period, the stream function pattern gradually evolves towards a single-wave structure ( Fig. 1h-o), which tends to be steady with continuously accumulated integration time. Compared with the initial pattern (Fig. 1a), the equilibrium state (Fig. 1o) differs not only in spatial pattern but also in intensity. This introduces an interesting question of how the transformation process from the initial pattern to the equilibrium state occurs. We analyse variables with conservation properties, such as the integral KE, square vorticity, and vorticity, to explore the variations (Fig. 2). As shown in Fig. 2, these variables do not meet the conservation requirements for the whole integration period. The integral KE curve (Fig. 2a) exhibits good conservation properties at the very beginning (left panel of Fig. 2) in which the influence of the initial pattern is more significant, and the stream function pattern evolves as an Fig. 3 The six equilibrium stream function patterns at t = 100, 000 analogical index cycle as mentioned above. The stream function no longer meets the conservation constraint with a significant climbing trend (middle panel of Fig. 2), implying that the integration results encounter a predictability problem due to errors and are not reliable. Stopping at this point, the calculation seems to be perfect. Fortunately, we continue to integrate for a longer period of time. A new conservation property (right panel of Fig. 2), significantly different from the previous one, eventually replaces the climbing trend, accompanying the stream function pattern approaching an equilibrium state. Both the integral square vorticity and vorticity curves experience similar variations (Fig. 2b-c). This demonstrates that the developed errors in the second period will not develop forever and disturb the computation but will gradually be suppressed by the strong conservation property of the flow system. As Fjörtoft (1953) once noted, the nature of the turbulence development in Eq. (1) is born from the conservation requirements it has to fulfil. Although the flow system eventually reaches an equilibrium state, the system does not regain its predictability after experiencing a long unpredictable period because slightly differing initial states could evolve into considerably different states (Lorenz 1963). We further artificially add six random normal distributed perturbations with a mean value of zero and a standard deviation of 0.001 to the initial pattern to simulate a sort of 'noise' superposed on the 'accurate' initial state. All perturbations are insignificant compared with the initial pattern. The six 'new' initial patterns are then integrated for the same time span. The stream function patterns initially evolve with little discrepancy (figure omitted), gradually diverging from each other and then eventually approaching different equilibrium states (Fig. 3). Although the six equilibrium states (Fig. 3) are different, they possess similar spatial patterns. Each of them exhibits a single-wave structure similar to the equilibrium state, and it seems that any one of them could be reproduced by the translation of one of the others by certain distances in the x and y directions. All the equilibrium states have a spatially similar single-wave structure, which demonstrates the largest scale of motion (wavelength) in the bounded region. This supports the arguments derived analytically, Fig. 4 The six equilibrium stream function patterns at t = 100, 000 which hypothesises that the larger scale (wavelength) is favoured over the smaller scale (wavelength) due to the strong conservation constraints in two-dimensional incompressible flow (Merilees and Warn 1975;Van Delden 1984). Therefore, we could infer that the stream function patterns with a single-wave structure might be one kind of eventual equilibrium state of the nonlinear BVE.
To further confirm the above analysis, we use larger wavenumbers (k = 10, l = 5) to construct a new initial pattern for smaller motion scales. We also use a random normal distributed initial pattern with a mean of zero and a standard deviation of 0.001 to simulate how the nonlinear BVE evolves with a complete random perturbation. We also apply the Lorenz (1960) initial pattern but set a fine grid resolution (d = 2π/80) and a coarse time resolution (dt = 0.01) to explore the influence of different integration scheme parameters. Finally, we apply the implicit finite difference integration scheme described in Eq. (6) to integrate with the Lorenz (1960) type initial value and the complete random normal distributed initial value, respectively, to test whether the phenomenon is scheme dependent. The eventual stream function equilibrium patterns for the above six cases (Fig. 4) are spatially similar to the previous equilibrium patterns (Fig. 3), further implying that the single-wave structure equilibrium states are ubiquitous for the numerical integration of the nonlinear BVE.
The 13 stream function patterns in equilibrium (Figs. 1o, 3, and 4) could inspire an intuitive conjecture that the eventual equilibrium patterns should have a similar expression as Eq. (11), namely, where A * and F * are intensity parameters; k * and l * are the wavenumbers in the x and y directions, respectively, and should be equal to one according to the above analysis; and θ 2 and θ 1 are the phases in the x and y directions, respectively. The specific values of the six parameters can be identified by applying Eq.
(12) to fit the equilibrium patterns through the least squares method. According to the results, the values of both k * and l * are very Fig. 5 Stream function patterns from t = 10, 000 to t = 100, 000 with a time interval of 10, 000. Streamlines are drawn at intervals of 0.1 close to one, further validating the single-wave structure speculation; both θ 2 and θ 1 seem to be random numbers, implying that the deterministic single-wave structure is indeterministic in phases. Since Eq. (12) can represent the equilibrium states, there should be no variation if the initial values are taken to continuously integrate the nonlinear BVE, regardless of how long the integration time lasts. Therefore, one of the 13 equilibrium states is selected (Fig. 1o) as the initial value to continuously integrate for the same time span. As expected, the stream function maintains its initial pattern and intensity during the whole integration period (Fig. 5). The integral KE, square vorticity, and vorticity curves slightly oscillate against the conserved values, showing good conservation properties (Fig. 6), and the period in which the conservation constraint is not satisfied (seen in Fig. 2) does not occur, further demonstrating that the flow system is independent of time after approaching the equilibrium state. The equilibrium states in Figs. 3 and 4 are also used as the initial values to integrate the BVE for the same time span. The results are similar (figure omitted).

Conclusion
In the present study, the nonlinear BVE has been directly solved by applying finite difference schemes with complete quadratic conservation features. The longterm numerical solutions show new characteristics that have not been revealed before. The stream function initially evolves under the influence of the initial value, portraying an analogical index cycle similar to the analytical prediction made by the LOM method. The integral KE, square vorticity, and vorticity conserve their values well during this period. After the period, the influence of the initial value is gradually weakened, and the conserved variables are not conserved again. This result is significantly different from the theoretical solution that predicts an endless index cycle. Computational instability will not be developed due to the accumulated error. In contrast, it will be gradually suppressed by the strong conservation constraint that the Fig. 6 The integral KE (a), square vorticity (b), and vorticity (c) curves for the whole integration period. To glorify the evolution features, the whole integration period is divided into three parts are bound by the vertical lines in each figure flow system obeys. Therefore, the stream function evolution is gradually independent of time and approaches an equilibrium state with a single-wave structure. The integral KE, square vorticity, and vorticity regain their conservation attributes, accompanied by the gradually steady stream function pattern. Six random patterns are added to the initial value to simulate types of errors that cannot be eliminated in the real world. The results suggest that the stream function patterns for the six cases initially evolve concurrently with the undisturbed case. However, the stream function patterns gradually diverge, and each eventually approaches an equilibrium state when the integration process continues. The eventual equilibrium states do have certain similarities despite the obvious differences. Each of them shows a single-wave structure in the domain with different phases in the x and y directions. Six cases are further designed to test the performances of the different initial patterns, spatial resolutions, and time resolutions, as well as an implicit difference scheme. The results are similar. The 13 equilibrium stream function patterns together prove the argument that the equilibrium states of the BVE are a single-wave structure with different phases.
Classical theory predicts an indeterministic characteristic of the deterministic system due to the nonlinearity (Lorenz 1963). The results of this study confirm this long-held view and further demonstrate that the indeterministic characteristics contain a deterministic component. For the nonlinear BVE, the eventual equilibrium pattern is indeterministic, but its spatial pattern is deterministic, i.e., a simple single-wave structure with random phases. This implies that we do not know the accurate future state of the flow system, but we do know what its spatial pattern is. Climate systems also obey many conservation laws. Therefore, this result offers the possibility of exploring long-term climate variability, which is currently a concern.
Funding information This study was jointly funded by the National Natural Science Foundation of China (Grants 41505042, 41805041); the National Program on Global Change and Air-Sea Interaction (GASI-IPOVAI-03); the National Basis Research Program of China (2015CB953601, 2014CB953903); and the Fundamental Research Funds for the Central Universities. The authors acknowledge the support from the National Marine Environmental Forecasting Center (KJHX2015190).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.