Inertial instability and phase error in Euler forward predictor‑ corrector time integration schemes: Improvement of modeling Great Lakes thermal structure and circulation using FVCOM

This study investigates the inertial stability properties and phase error of numerical time integration schemes in several widely-used ocean and atmospheric models. These schemes include the most widely used centered differencing (i.e., leapfrog scheme or the 3-time step scheme at n-1, n, n+1 ) and 2-time step ( n, n+1 ) 1 st -order Euler forward schemes, as well as 2 nd -stage and 3 rd - and 4 th -stage Euler predictor-corrector (PC) schemes. Previous work has proved that the leapfrog scheme is neutrally stable with respect to the Coriolis force, with perfect inertial motion preservation, an amplification factor (AF) equal to unity, and a minor overestimation of the phase speed. The 1 st -order Euler forward scheme, on the other hand, is known to be unconditionally inertially unstable since its AF is always greater than unity. In this study, it is shown that 3 rd - and 4 th -order predictor-corrector schemes 1) are inertially stable with weak damping if the Coriolis terms are equally split to n +1 (new value) and n (old value); and 2) introduce an artificial computational mode. The inevitable phase error associated with the Coriolis parameter is analyzed in depth for all numerical schemes. Some schemes (leapfrog and 2 nd -stage PC schemes) overestimate the phase speed, while the others (1 st -order Euler forward, 3 rd - and 4 th -stage PC schemes) underestimate it. To preserve phase speed as best as possible in a numerical model, alternating a scheme that overestimates the phase speed with a scheme that underestimates the phase speed is recommended. Considering all properties investigated, the leapfrog scheme is still highly recommended for a time integration scheme. As an example, a comparison between a leapfrog scheme and a 1 st -order Euler forward scheme is presented to show that the leapfrog scheme reproduces much better vertical thermal stratification and circulation in the weakly-stratified Great Lakes.


Introduction
There are three types of numerical or computational stability issues associated with time integration and spatial-discretization schemes in ocean and atmospheric general circulation models. The first stability type is the well-known fast gravity wave stability, which usually applies to the external mode if a model is split into a fast gravity wave mode and a slow wave mode. The fast gravity wave constraint is the Courant-Friedrichs-Lewy (CFL) criterion, which is used to ensure that a linear model is stable (Beckers and Deleersnijder 1993;Wang 1996). The second stability type is linked to non-linear terms, often referred to as non-linear instability associated with spatial discretization. The non-linear instability is often controlled by the energy-conserving spatial numerical schemes (Wang 1996). The third stability type, which is least known and therefore often overlooked, is the inertial instability associated with Coriolis parameter (f) as a result of Earth's rotation (O'Brien 1986;Wang and Ikeda 1997a). The inertial instability is associated with time integration schemes linking to the Coriolis terms, often in the external mode of the mode-split models (Beckers 1999). Although the CFL criterion (usually ranging from several to tens of seconds) is much stricter (smaller) than the time step constraint for inertial stability (i.e., |F|=|f|Δt≤1; Δt<10 4 seconds), the inertial instability can still occur if the Coriolis terms are incorrectly discretized in time in a numerical scheme, leading to inconsistencies with the original physical system (i.e., differential equations). This inertial stability is not explicitly discussed in user's manuals in most hydrodynamic models, as if to imply that the models automatically or naturally meet the inertial stability criterion as long as the CFL criterion is met. Unfortunately, some widely-used models have issues with inertial instability that go unnoticed by novice users, creating substantial, propagating errors in f-related physical processes. Wang (1996) proposed that global stability should be used to gauge the stability of a model, which includes the three types of stabilities mentioned above as well as any other physical constraints in a model (e.g., h min +ζ max >0, where h min is water depth and ζ max is surface water elevation along the coast). Wang and Ikeda (1997a) further defined that the inertial stability is constrained by both a numerical constraint (|F|=|f|Δt≤1, which is automatically met as long as the CFL criterion is met) and a physical constraint, such that the amplification factor (AF) or eigenvalue ( ) of the discretized scheme must be equal to or less than unity, i.e., | |≤1, where | |=1 is defined neutrally stable and | |<1 is stable with damping. Using the inertial stability criterion, Wang and Ikeda (1997a) found that the Euler forward scheme used in ECOM_si (Estuarine and Coastal Ocean Model with semi-implicit scheme, Blumberg 1991) is inertially unstable, which was replaced with a stable 2 nd -stage predictor-corrector scheme (Wang and Ikeda 1997b).
Nowadays, ocean and atmospheric models are made open to broad users. Users often take it for granted that all the released models are ready to be applied to their research and applications, without understanding the stability properties of the numerical schemes used in the models. It should be kept in mind that the substantial differences of geophysical fluid dynamics (GFD) from non-rotating ordinary fluid mechanics (OFM) are 1) rotation and 2) stratification. Particularly in the oceans, the vertical stratification is much stronger and more permanent, than in the atmosphere and freshwater lakes, where frequent deep convection can homogenize density stratification. Rotation in the GFD not only leads to the challenge of inertial instability of a numerical scheme (O'Brien 1986;Wang and Ikeda 1997a, b), but also inevitably introduces phase error in any f-related waves such as inertial waves, coastal trapped waves, tides, Rossby waves, and large-scale planetary waves. Therefore, correct and accurate discretization of the Coriolis terms is essential for any models to reproduce realistic physical processes related to f. It is noted that spatial discretization also affects the wave properties, leading to error in phase and amplitude of a wave (Wang 1996). These processes range from turbulence and small-scale mixing to inertial motion (Austin 2013) to large-scale planetary waves, because they all contain the Coriolis parameter. Therefore, a numerical scheme that is computationally stable in the non-rotating system is not necessarily a stable scheme in the rotating system. In other words, stability, truncation error, and consistency of a stable numerical scheme in the non-rotating system (Lemarie et al. 2015) should be re-examined in the rotating system.
Since the 1990s, several ocean and atmospheric models were widely used, such as Princeton Ocean Model (POM, Blumberg and Mellor 1987), Regional Ocean Model System (ROMS, Song and Haidvogel 1994), MITgcm (Adcroft et al. 2009), Finite Volume Community Ocean Model (FVCOM, Chen et al. 2013), Modular Ocean Model (MOM, Griffies et al. 2009) and Weather Research and Forecast Model (WRF, Skamarock et al. 2008). However, some models introduce several new time integration schemes that were not discussed by Wang and Ikeda (1997a). For example, FVCOM uses 2-time step, 1 st -order Euler forward scheme in the internal (slow wave) mode, and the 4 th -stage Runge-Kutta scheme in the external (fast wave) mode. The ROMS and MITgcm also have several options for time integration schemes (Shchepetkin and McWilliams 2005), which include the leapfrog scheme (default), Euler forward scheme, and 3 rd -stage Adams-Bashforth scheme.
In this paper, we focus only on the inertial instability and its associated phase error caused by the discretized differencing schemes, compared to their analytical (original continuum) solutions. The major purpose of this study is to investigate inertial instability (i.e., computational instability of the inertial mode) of widely used time integration schemes and the related phase error associated with the Coriolis parameter. Though the phase error exists in all time-integration schemes, it can be minimized if we understand the sources of the error. Differences in integration scheme performance (and stability) are illustrated using FVCOM, with model results generated using the default 1 st -order Euler forward scheme compared to those generated using a newly implemented 2 nd -order leapfrog scheme in the Laurentian Great Lakes.

Inertial instability and diffusive behavior observed in the default schemes in FVCOM
Instability issues in FVCOM were first identified when implementing the model in the Laurentian Great Lakes. The model was run using default time integration schemes, such that the external mode is solved using a 4 th -stage Runge-Kutta scheme and the internal mode is solved using a 1 st -order Euler forward scheme. The realistic 5-lake model simulation was forced using 3-hourly NARR (North American Regional Reanalysis) atmospheric data with a cold start on January 1, 1995 where u=v=0 and T=2 o C everywhere. Although the model was able to reproduce reasonable, albeit overly diffuse, surface temperatures in each lake, physically unrealistic oscillations were observed in near-surface currents throughout the year. Instability was especially pronounced in the stratified period ( Fig. 1; day 182, 241), when simulations showed high-speed current "fronts" propagating across the deeper waters of Lakes Superior, Michigan, Huron, and Ontario. Importantly, modeled current speeds and directions were inconsistent with weak atmospheric forcing, and oscillations did not match well-established lake seiche periods or structures. Simulated oscillations were only suppressed in especially shallow waters, like Lake Erie (mean depth: 19m), where bottom friction was able to dampen inertial instabilities.
Model results highlighted significant flaws in FVCOM parameterizations, motivating additional work to identify and correct issues with the default model framework.
Model instability is commonly addressed by increasing the vertical and horizontal viscosity, which effectively stabilizes the model by dampening spurious current and Fig. 1 Simulated lake surface temperature superimposed with lake currents at 10m. Examples of artificially amplified inertial currents (i.e., unstable inertial motion) are shown for day 182 (June 2; upper panel) and day 241 (August 30; lower) in 1995 1 3 temperature oscillations. While this method is effective, it leads to unrealistically high mixing rates, and the viscosity required for stabilization may even be larger than physically permitted. Viscosity stabilization also leads to "smearing" of physical processes, reducing (or eliminating) physically meaningful temperature and current gradients. This issue is especially problematic in permanently stratified oceans, where strong density gradients control many physical and biogeochemical processes. Instead of using excessive viscosity to treat the symptoms of model instability, the underlying cause of the instabilities should be addressed directly. Here, we hypothesize that instabilities observed in FVCOM are driven by unstable time integration schemes, which manifest as computational instability in the inertial mode.
This hypothesis was tested using an idealized FVCOM simulation in Lake Erie. The model was initialized with homogeneous temperature (4 °C) and current (u=v=0) fields and external forcing was limited to heat flux (~80 Wm -2 ) in a circular area in the central basin (Fig. 2a). As such, buoyancy with inertial motions was the driving force in the development of horizontal velocity fields. Two model runs were conducted for 40 days under the same thermal forcing with both the default (Fig. 2b) and leapfrog (Fig. 2c) scheme. The difference between the default and leapfrog time integration schemes (Fig. 2d) shows obvious warming inside the circle, while along the edge of the circle warming and cooling are alternative. The result indicate that default integration schemes produce overly diffusive behavior in temperature (Fig. 2d), which may be directly linked to the truncation error, as will be discussed shortly in "Application to the Great Lakes: Improving vertical thermal stratification" and Appendix C. Below, we provide a theoretical analysis of the inertial stability of various time integration schemes "Theoretical analysis of inertial stability", culminating in a recommendation to use the leapfrog scheme to improve inertial stability. The leapfrog scheme is then applied to FVCOM simulations in the Great Lakes "Application to the Great Lakes: Improving vertical thermal stratification", followed by a proposed phase error correction "Proposed phase error correction" and study conclusions "Conclusions and future efforts".

Theoretical analysis of inertial stability
To reveal the inertial instability in time integration schemes, we introduce a well-known, simplified, pure inertial oscillation equations (system), which are imbedded in both 2D and 3D shallow equations, with no viscosity (or friction) as follows:

3
where w=u+iv, a complex velocity with i 2 =-1 and f is the Coriolis frequency (=2Ωsinφ~10 -4 s -1 , where φ is the latitude and Ω=7.292x10 -5 s -1 ). The reasons we chose this pure inertial model are 1) this inertial model is embedded in the 2-D shallow water equations and the 3-D Reynolds stressaveraged geophysical fluid dynamic system, 2) the analytical solution is known and can be objectively compared to the differencing equations, and 3) unlike the advection-diffusion equation, the pure inertial system has no dissipation (friction or viscosity) so it can be used to examine time integration schemes in an accurate manner. The exact solution to Eq. 1 is The amplitude of the exact solution is always unity with no viscosity, i.e. |w(t)|=√[u(t) 2 +v(t) 2 ]≡1. This is the true criterion to gauge any time integration scheme in the differencing equations used by ocean and atmospheric models with Coriolis force (i.e., in a rotating system) (Wang and Ikeda 1997a). Therefore, if a numerical scheme discretizing (partial differential equations, PDEs) dynamic system (1) preserves the inertial motion with an amplitude of unity (i.e., AF=1), it is defined as neutral inertial stability. If the AF for a numerical scheme is less than unity, it is defined as a stable scheme with numerical damping and if AF is larger than unity, it is defined an unstable scheme for the inertial oscillatory (rotational) system.
It must be pointed out that the inertial wave frequency is f, irrespective of the wave numbers; thus, it is a nondispersive wave. Therefore, theoretically, under no viscosity condition, the clockwise (anticyclonic) inertial wave/ motion [Eq. (1)], once it is excited, will continue forever with unity amplitude. Therefore, the discretized numerical scheme of Eq. (1), under no viscosity condition, should preserve both the inertial motion without decaying and the non-dispersive property. Wang and Ikeda (1997a) have investigated the inertial instability and phase error of a pure inertial motion with (2) (u, v) = cos(ft), − sin(ft) , or w = e −ift a linear friction (dissipation) using a centered differencing scheme, first-order Euler schemes, and second-stage predictor-corrector schemes. The results can be applied to a pure inertial system without friction by simply removing the friction terms. Since these three schemes have been investigated in depth in Wang and Ikeda (1997a) and Durran (2010), they will not be repeated here. However, for the sake of completeness and comparison, we list these schemes in Table 1, 2, 3, 4 and 5 and use the results for comparison throughout this study.

Third-stage Euler forward predictor-corrector scheme
In the default scheme of FVCOM, Euler forward, two-time stepping (n+1, n), 3 rd -or 4 th -stage Runge-Kutta scheme is used in the external mode with Coriolis terms placed at time step n. Since Rung-Kutta scheme (along with Adams-Bashforth scheme) belongs to the same family of the predictorcorrector scheme, we use the simplified predictor-corrector scheme applied to pure inertial system (1) to reveal the stability property of FVCOM. The derivation of stability analysis is shown in detail in Appendix A. If β=0, 1/2, and 1, the amplification factor (AF, i.e., eigenvalue, λ) becomes The 3 rd -stage Euler forward scheme (β=0) has an unstable error growth rate of 1+F 2 /2, equivalent to the 1 st -order Euler forward scheme described in Wang and Ikeda 1997a, see their Fig. 8) and Table 2, since it uses the same old value (at n) for its corrector scheme. When β=0.5, unlike the 2 nd -stage Euler centered scheme (Table 2), the 3 rd -stage Euler centered scheme is stable with a weak dampening and AF decreasing to ~1-F 4 /8, which classifies the 3 rd -stage Euler centered (β=0.5) scheme as a  Chen et al. (2013) Internal mode: 2-time step Euler forward: Δt + fu n = 0, External mode: 3 rd or 4 th -stage Runge-Kutta 1 3 stable scheme. In other words, in a rotating system, the 3 rd -stage predictor-corrector (PC) scheme is better than the 2 nd -stage Euler PC scheme in terms of accuracy, although inertial waves are slightly dampened with this scheme. When β=1, similar to the 2 nd -stage scheme (Table 2), the AF of the 3 rd -stage Euler backward scheme is smaller than unity; the scheme is stable but with significant dampening of the inertial waves.
The phase error due to this scheme can be expressed as, by inserting (20) to (16), Note that, similar to the 2 nd -order centered/leapfrog schemes (Table 2), the 3 rd -stage predictor-corrector scheme also introduces a computational mode, as shown in Eq. (20) and Eq. (4). This computational mode adds a weakness to the 3 rd -stage predictor-corrector scheme, which should be controlled numerically using a frequency filter similar to the Robert-Asselin-William (RAW) filter (Asselin 1972;Williams 2009). Figure 3 shows the normalized inertial wave frequency (ω 1 , physical solution; Fig. 3a) versus F=f∆t . Both forward (β=0) and backward (β=1) schemes produce slower frequency, and the centered scheme (β=0.5) has a much better fit to the exact solution, although still slightly slower than the exact solution. The unphysical computational mode is just symmetric to the physical solution, but around the line of -1 (Fig. 3b), indicating that the artificial waves rotate cyclonically (opposite to the physical solution). Figure 4 shows the inertial wave frequencies for both physical (ω 1 , physical solution) and computational (ω 2 , non-physical solution) modes for the centered scheme (β=0.5). Therefore, a filter, similar to the RAW filter discussed above, should be applied to control non-physical waves and avoid non-physical modes ruining the physical solution, further discussion of which is beyond the scope of this paper.

Fourth-stage Euler predictor-corrector scheme
The derivation of the stability analysis is shown in Appendix B. When β=0, 0.5, and 1, the AF becomes Uncond. unstable, Weakly unstable, Uncond. unstable, Weakly damping, Damping Yes Uncond. unstable, Weakly damping, Damping Yes It is noted that the 4 th -stage Euler forward scheme (β=0) has an unstable error growth rate of 1+F 2 /2, identical to the 1 st -order Euler forward scheme shown in Table 2, which uses the same old value (at n) for its predictor and corrector steps. When β=0.5, the 4 th -stage Euler centered scheme is stable, similar to the 3 rd -stage scheme, but with weaker dampening than the 3 rd -stage scheme since the AF decreases to ~1-F 6 /32. When β=1, the AF of the 4 th -stage Euler backward scheme is smaller than unity (i.e., stable) but always with significant dampening of the inertial waves.
In a similar way, we obtain two frequency roots by inserting (31) to (16): with the normalized phase error defined as Note that the second mode is a computational mode, similar to that produced by the leapfrog scheme and the 3 rd -stage Euler scheme (Table 2). Figure 5 shows the normalized inertial wave frequency (ω 1 , physical solution; Fig. 5a) versus F=fΔt. Both forward (β=0) and backward (β=1) schemes produce slower phase speeds, particularly for the Euler backward scheme. The centered scheme (β=0.5) has a much better fit to the exact solution, although still slightly slower than the exact solution. The unphysical mode is symmetric to the physical solution, but around the line of -1 (Fig. 5b). Figure 6 shows the inertial wave frequencies for both physical (ω 1 ) and computational (ω 2 , non-physical) solutions for the centered scheme only (β=0.5). Therefore, a filter, similar to the RAW filter, should be applied to remove non-physical waves and prevent the non-physical solution from ruining the physical solution, which is beyond the scope of this paper. Figure 7 summarizes the physical frequencies from all five schemes with β=0.5. Both the leapfrog and 2 nd -stage Euler PC schemes (Wang and Ikeda 1997a) overestimate the phase frequency. While the 1 st -order 3 rd -stage predictor-corrector Euler ꞷ 3 rd -stage predictor-corrector Euler ꞷ Euler centered (Wang and Ikeda 1997a) and 4 th -stage Euler centered schemes underestimate the phase frequency, the 3 rd -stage centered Euler scheme produces a relatively better phase frequency than the other two, although with a slight underestimation.

Application to the Great Lakes: Improving vertical thermal stratification and circulation
The hydrodynamic model used in this study is based on the unstructured grid FVCOM (Chen et al. 2013), version 3.1.6. FVCOM solves the primitive equations using the model splitting method: The external mode solves vertically integrated transport equations in which water elevation is solved explicitly using a shorter time step, while the internal mode solves 3-D governing equations using a longer time step. Both internal and external modes include the Coriolis terms, differing from Princeton Ocean Model (POM; Blumberg and Mellor 1987) in which the Coriolis terms are applied only to external mode. FVCOM solves the flux form of the governing equations in unstructured triangular volumes using a discrete-flux scheme following the finite volume approach. This model was applied to simulation of Great Lakes circulation and ecosystem using excessive vertical mixing (Luo et al. 2012;Bai et al. 2013;Rowe et al. 2015;Anderson et al. 2018) The default time integration schemes in FVCOM are the 4 th -stage (modified) Runge-Kutta (RK4) method for the external mode, and the 1 st -order Euler forward scheme with the Coriolis terms discretized at time n for the internal mode. As mentioned above, the Euler forward scheme (Tables 1 and 2; Wang and Ikeda 1997a) is inertially unconditionally unstable, because the AF (or the eigenvalue of the differencing equations) is always greater than one: | |=1+ f 2 Δt 2 >1. The Euler forward scheme can be made stable by using excessively high viscosity, which, nevertheless, can smooth all the realistic physical phenomena such as mesoscale eddies (Wang and Ikeda 1997b). Similarly, very large diffusivity coefficients can diffuse horizontal and vertical thermohaline structures such as vertical stratification and ocean fronts (Luo et al. 2012;Bai et al. 2013). Therefore, in addition to keeping the original default schemes untouched, we implemented the centered differencing (leapfrog) scheme with the RAW filter (Beletsky and Schwab 2001;Wang et al. 2010;Fujisaki et al. 2012Fujisaki et al. , 2013Dupont et al. 2012) to both modes of FVCOM (Fujisaki-Manome and Wang 2016) because both modes in FVCOM suffer inertial instability. Therefore, both modes in the modified version have neutral inertial stability, and the improvement affects both internal (baroclinic flow) and external (barotropic flow) mode. Note that both FVCOM and POM use mode-split approach. However, POM includes the Coriolis terms only in the external mode, while FVCOM applies the Coriolis terms to both external and internal modes. If numerical schemes are stable, the physical processes can be reproduced with both POM and FVCOM, because the adjustment between 2D and 3D systems is conducted for each time step. However, the default schemes are inertially unstable in both modes in FVCOM, smaller internal and external time steps must be used, along with excessive viscosity, to avoid model blowup. This is why the common ratio of the external to internal time step is ~10 in the default schemes, while ~20 in the leapfrog scheme.
Thus, in this work, FVCOM was modified such that the model could be run using one of two different time integration schemes: (1) the default schemes or 2) the leapfrog scheme proposed above.
The modified FVCOM with centered differencing has been successfully applied to the Great Lakes and coastal ocean with validation using measurements Bai et al. 2020;Cannon et al. 2023, this issue). Cannon et al. (2023) ran this revised FVCOM with the leapfrog scheme, coupled with an ice model, from 1979 to 2021 without any data nudging and assimilation, unlike Xue et al. (2017) and many others who nudged observed water temperature to the model temperature. Figure 8 shows the simulated vertical temperature profiles over a seasonal cycle of 1998 in south central Lake Michigan under the 3-hourly forcing of North America regional reanalysis (NARR). The observed temperature shows a strong thermocline with a mixed layer depth around 20m (Fig. 8a). Using the default schemes of FVCOM (RK4 scheme for the external mode and the 1 st -order Euler forward scheme for the internal mode), the temperature profile drifts significantly away from the measurement (Fig. 8b) with a single year of simulation. This can be attributed excessively large vertical and horizontal viscosity produced in the turbulence closure models to dampen the inertial instability (Wang and Ikeda 1997b) and first-order truncation error induced by the first-order Euler forward scheme in the internal mode. The truncation error is the physically meaningful bi-harmonic viscosity (mixing), as shown in Appendix C. The excessive viscosity then destroys the intensified (physical) thermocline. This is why the temperature structure disperses over both the time and with depth (i.e., the entire water column). Indeed, the thermocline is nearly destroyed, and the upper mixed layer depth is very small or shallow--less than 5m. Figure 8c shows the same simulated temperature, but with the centered differencing scheme in both internal and external modes, similar to POM and ROMS. It is clear that the centered differencing scheme reproduces much better temperature structure and stronger thermoclines than the original default schemes, even though the modelled thermocline steepness is still slightly weaker than observations. Figure 8c compares very well with the measurement (Fig. 8a), while default schemes destroy the stratification with diffusive vertical temperature structure, i.e., weak vertical temperature stratification.
We ran the five-lake model  with both default and centered differencing schemes from 2012-2018 with no ice. The setting and forcing (NARR) remain identical in both cases, except for the different time integration schemes. Vertical viscosity was calculated using the 2.5 turbulence closure model, while the horizontal viscosity was calculated using the Smagorinsky parameterization. Figure 9 shows seasonal average lake surface temperatures (LST) for summer and autumn in 2015 (Fig. 9a) modeled using both schemes. Both simulated LST distributions look similar. However, the difference between the default schemes and the leapfrog scheme for all four seasons (Fig. 9b) shows that the Euler forward scheme produces significant larger warming in spring and summer across all five lakes. Maximum warming occurs in summer, with differences as high as 2 -2.5 o C, in deep waters. In autumn, warming produced by the default schemes occurs in deep waters and cooling occurs in shallow waters. The default schemes produce excessive cooling up to 1.5 o C in winter, compared to the leapfrog scheme, indicating that the default schemes reproduce excessive ice in winter due to the cold bias.
To validate the difference between the default schemes and leapfrog scheme with no surface wind-wave mixing, we chose a thermistor chain mooring deployed in western Lake Superior (see the western mooring in Fig. 1 Fig. 7 Enlarged comparison (between 10 -1 <F<10 0 =1) of phase frequencies of all physical solutions in the numerical schemes discussed located between 25 and 40m (Fig. 10a). The model with the leapfrog scheme reproduces a similar mixed layer (depth: ~20m) and thermocline thickness (20 -30m) with no wave mixing (Fig. 8b). Temperature differences between each scheme and the observations are shown in Figure 10(c, d). The default schemes produce excessive warming not only in the upper layer, but also in the lower layer (Fig. 10c), while the leapfrog scheme produces warming only in the upper layer (Fig. 10d). The upper layer warming is likely due to heat flux calculations, while the lower layer warming observed with the default schemes is linked to scheme truncation error that is equivalent to bi-harmonic mixing, as shown in Appendix C. This indicates that the mass is not conserved in the default schemes, as shown in Appendix C. Figure 11 shows the simulated 10-m current distribution using the leapfrog scheme on day 182 and 241 of 1995, which can be vividly compared to the simulation using the default schemes shown in Fig. 1. The inertial instability disappears in Fig. 11, because the leapfrog scheme is of neutrally inertial stability and of second order accuracy in time that does not produce numerical viscosity. The systematic validation of lake circulation, thermal structure using both satellite and in-situ measurements can be found in Bai et al. (2013), Li et al. (2021), and Cannon et al. (2023, this issue).

Proposed phase error correction
When considered in conjunction with Wang and Ikeda (1997a), the above investigation suggests that the Coriolis parameter used in the shallow water equations (i.e., in a rotating system) inevitably introduces phase errors during model integration, regardless of integration schemes. The differencing equations produce phase errors due to discretization at every time step of model integration (O'Brien 1986;Zhang et al. 1987;Wang and Ikeda 1997a;Durran 2010), with errors accumulated (and compounded) over the entire model run duration. The true frequency, f, of a nondispersive inertial wave is strictly held in the corresponding differential equations. Thus, phase error resulting from any numerical schemes should be minimized.
To make a quantitative comparison, Table 3 shows the relative phase error associated with Coriolis frequency (f) in the following form: [ω 1,2 /(-f)-1]x100 (%). The negative (positive) values indicate that the phase of a numerical scheme is slower (faster) than the true inertial frequency. It is clear that the phase derived from the Euler forward scheme, PC3, and PC4 schemes is slower than the physical solution, while leapfrog and 2 nd -stage PC schemes produce faster phase speeds than the physical solution.
Generally speaking, a coastal ocean model time step ranges from a few to tens of seconds in the external model (fast wave), and ~100 seconds in the internal mode. Thus, we take these two time steps, Δt=10s and 100s, as illustrative examples. When Δt=10s, phase error produced by the leapfrog and PC2 schemes is +2x10 -5 %, which is twice that derived from the Euler PC3 and PC4 schemes (-1x10 -5 %). Similarly, when Δt=100s, error produced by both the leapfrog and PC2 schemes is again double that derived from the Euler forward, PC3, and PC4 schemes (+1.67x10 -3 % vs -0.83x10 -3 %). This relative error is two orders of magnitude larger than that when Δt=10s. The same feature can be seen when Δt=1000s, and so on and so forth.
The phase error can significantly affect simulation results, especially for simulations of physical processes associated with the inertial frequency, including currents, tides, storm surges, and Rossby waves. When a model uses the leapfrog and PC2 schemes (see Fig. 5, and Wang and Ikeda 1997a), the phases (exemplified by peaks and troughs) of these waves are always faster than the observed waves. Contrarily, when a model uses the Euler forward, PC3, or PC4 schemes, the phases of these waves are always slower than the observed waves (see Fig. 8 of Rowe et al. 2015). Based on Table 3, phase error in PC3 is almost identical to PC4 when Δt<400 s. As Δt increases, the phase error in PC3 is smaller than that in PC4. As to AF, although PC4 is more accurate than PC3 (F 2 : 1-F 6 /32<1 vs. 1-F 4 /8<1), with each performing as weakly damping schemes, the accuracy of PC3 is sufficient for most model applications. For example, if Δt=100s, F=10 -2 ; then PC3's AF=1-F 4 /8≈1-10 -7 , while PC4's AF=1-F 6 /32≈1-10 -13 . Therefore, considering both AF and phase error, PC3 schemes are advantageous compared to PC4. Another benefit of the PC3 scheme is that the PC3 schemes save one-step calculation and storage compared to PC4.
In order to cancel the (faster) phase error introduced by the leapfrog scheme, Ford et al. (1990) and Wang (1998) alternatively used the Euler backward scheme once after every 20 steps of integration using the leapfrog scheme to slow down the phase of tidal waves and current. The idea Fig. 9 FVCOM-simulated lake surface temperature (LST) using default (Euler forward) scheme (upper panel of Fig. 8a) and leapfrog scheme (lower panel of Fig. 8a) in summer (JAS) and autumn (OND). The LST difference between the default schemes and the leapfrog scheme in four seasons (Fig. 9b) is to correct the wave phase that is accelerated using the leapfrog scheme by the Euler scheme that slows down the wave phase. Although data assimilation methods are widely used to correct the model behavior toward the observation, it is highly recommended that before using data assimilation to constrain a model, the model should be physically correct and numerically accurate. The physical correctness includes model parameters, mixing schemes, and parameterizations of some sub-scale physical processes, etc. The numerical accuracy should include a very basic principle, that is, the model should be neutrally stable or least damping. If a model cannot reproduce basic ocean circulation, temperature and salinity structure, including stratification and fronts, without data assimilation (i.e., the model has strong drift away from the observation), then the model must be either extremely dissipative (AF<1) or unstable (AF>1). In other words, if a model is unable to accurately reproduce basic circulation, temperature and salinity structure without the help of a strong data assimilation (i.e., nudging), then the model should not be used for a forecast system before both physical correctness and numerical accuracy in the model are guaranteed.
A possible application of Table 3 is to program an existing model to reduce the phase error. For example, if a model uses a PC3 scheme for the time integration at Δt=100s for 20 time steps, which slows down the wave phase by -0.00083%x20=-0.0166%, one can alternatively use a PC2 scheme for the time integration in a single time step with Fig. 10 a) Observed water temperature during May 12-December 12, 2012 in Lake Superior; b) FVCOM-simulated lake temperature using leapfrog scheme; Lake temperature difference c) between the default (Euler forward) schemes and observation, and d) between the leapfrog scheme and observation Δt=300s to correct the phase error because the wave phase is accelerated in the PC2 scheme by 0.015%. Alternatively, if a model uses the second-order leapfrog scheme for the time integration at Δt=100s for 10 time steps, which accelerates the wave phase by 0.00167%x10=0.0167%, then one uses the first-order Euler centered scheme for the time integration once at Δt=400s to correct the phase error because the wave phase is slowed down in the Euler centered scheme by -0.0133%. This latter method was previously employed by both Ford et al. (1990) and Wang (1998).
If an ocean tidal model is run for one day continuously without any correction (or data assimilation) using the Euler forward (PC3 or PC4) scheme with ∆t=100 seconds, the simulated wave phase would have a 51 second (~1 minute) lag to the observed tidal wave. In a similar way, if an ocean model with an Euler forward (PC3 or PC4) scheme is used to simulate ocean circulation with mesoscale eddies (i.e., Rossby waves), the propagation of the eddies would have a 51 second lag to the observed eddies over a one-day simulation. That means that one-month simulations would cause 25 minutes of lag, and one-year simulations would result in about 5 hours of lag. The phase lag would become more serious by accumulation if the model is spun up for 10 years (50 hours of lag) and applied to predict or project multi-year and decadal conditions of ocean circulation with mesoscale eddies.
To provide more possible options for using Euler forward and backward schemes, Table 4 and Table 5 provide the phase errors when β=0, and β=1, respectively vs leapfrog scheme. One can use these alternative schemes to minimize the phase error in the model simulations.
In summary, a model with inertial instability cannot simulate correct dynamic phenomena, since the model always needs excessive numerical viscosity to stabilize the simulation, which would drift away from the observation. Similarly, to run a model without correcting the phase, the Table 3 A summary of the dependence of phase errors of different time integration schemes on time step for inertial oscillation for the 1 st -order Euler centered (β=0.5), 2 ndorder centered (leapfrog) differencing schemes, PC2, PC3, and PC4: F=fΔt, where f is the Coriolis parameter (taken 10 -4 s -1 ), and Δt (in seconds) is the time step for each model integration. The error is given by the normalized error, ω/(-f), relative to the true solution, 1, i.e., [ω/(-f) -1]x100 %. Positive (negative) values indicate the faster (slower) phase speed than the exact solution modeled phase speed of the f-related waves can also drift away (i.e. accelerate or decelerate) from the observed wave phase speed. A small instability and a small phase error at each time step can be accumulated to create large errors, with increased drift from observations over long-term integrations. This is particularly critical to seasonal prediction, decadal projection, and long-term climate projection (Cannon et al. 2023, this issue). Nevertheless, it can also seriously affect the prediction skills for a short-term (5-to 7-day) forecast of storm surges (Anderson et al. 2018) if a model drifts away from observations too quickly.

Conclusions and future efforts
The inertial stability is defined by both time-step constraint (F=|f|Δt≤1) and physical constraint, in which the eigenvalue or AF must not be greater than unity, with neutral stability ideally defined as AF≡1. Since the CFL condition for surface gravity waves (or internal wave, Lemarie et al. 2015) is much stricter than F=|f|Δt≤1, the time constraint is not an issue for inertial stability. Based on the above investigations of numerical schemes on the discretization of Coriolis terms, AF, and phase errors, the following conclusions may be drawn: 1) The two time-step Euler forward scheme (β=0) is always inertially unconditionally unstable, since the AF is always 1+F 2 /2>1 (see Table 3 of Wang and Ikeda 1997a and Table 2.2 of Durran 2010). To cure this problem, the velocity in the Coriolis terms must be split into both the old value (at n) and the new or predicted value (at n+1) (i.e., β=0.5 with an equal weight), which is the Euler centered scheme (see chap. 2 of Durran 2010). 2) Although the 1 st -order and 2 nd -stage Euler schemes produce a sole physical solution, the 3 rd -and 4 th -stage Euler forward PC schemes produce both a physical solution and a computational mode. The computational mode should be effectively controlled using a frequency filter, similar to the RAW filter for the centered differencing scheme, rather than simply by increasing viscosity, which would dampen other physical processes. 3) No matter what numerical schemes are used in the rotating system (i.e., f≠0), phase error is inevitably introduced in the discretized differencing equation. The leapfrog scheme and the 2 nd -stage predictor-corrector (PC2) scheme overestimate (accelerate) the phase speed in any waves associated with Coriolis parameter. The 1 st -order Euler forward, PC3 and PC4 schemes always underestimate (decelerate) the phase speed. To overcome this intrinsic shortcoming, the only solution is to alternatively use one scheme of overestimating phase speed and the other scheme of underestimating it. The alternative ratio depends on the magnitude of phase error shown in Tables 3-5. Note that data assimilation can pull the model back to the observed data, but it cannot substantially cure the inherent phase error in the numerical scheme. A model's deficiency should be cured from the root cause, rather than from the symptom. 4) Although the PC4 schemes are more accurate than the PC3 schemes, PC3 schemes are sufficiently accurate and the phase error of PC3 is smaller than that of PC4. 5) The Great Lakes FVCOM with default integration schemes (external mode: RK4; internal mode: 1 st -order Euler forward) produces an overly diffusive thermocline and excessively shallow upper mixed layer, while the leapfrog scheme produces a steeper, more strongly stratified thermocline. LST is significantly warm biased in the summer and spring under the default integration schemes, with cool biases in the winter. Compared to measurements in western Lake Superior, model results from the default integration schemes produced excessively warm water temperatures in the deeper layers.
Biases are linked to a1 st -order truncation error related to the bi-harmonic viscosity generated from by the default Euler forward schemes compared to the non-physical, 2 nd -order truncation error produced by the leapfrog scheme (Appendix C). The lake circulation simulated by the leapfrog scheme (Fig. 11) is free of inertial insta-bility, compared to that simulated by the default schemes ( Fig. 1) In this study, we apply the simplified predictor-corrector schemes to the pure inertial system [Eq. (1)] to prove that FVCOM default schemes, which treat the Coriolis terms at time (n) explicitly in both modes, are inertially unstable. This led to our implementation of the leapfrog scheme in FVCOM, which led to dramatic improvements in model stability, especially over long model runs. Several other time integration schemes may also produce stable, or quasi-stable simulations in FVCOM, but additional analysis would be required before implementation. For example, we did not make an effort to examine the overall stability criteria of 3 rd -and 4 th -stage Runge-Kutta schemes in the rotating system (Ketcheson 2010). However, stability areas and criteria for Runge-Kutta (as well as Adams-Bashforth) schemes were largely derived in a non-rotating system, thus a systematic analysis would be required to re-examine the stability areas and criteria of these two-time stepping, Euler forward schemes before implementation in FVCOM.
In addition to stability, the energy-conservation of time-space discretization schemes (Wang 1996) should also be considered for conducting long-term simulations using FVCOM (Cannon et al. 2023). It is worth noting that for a pure advection scenario in a two-time stepping, Euler forward scheme, energy cannot be conserved (see Appendix D). The discretized solution for the pure advection is unstable without the help of viscosity/diffusivity, since the eigenvalue (or the AF) of the 2-time stepping Euler forward (upwind) scheme is always greater than unity (see the derivation of Appendix D). Therefore, a joint examination of both time integration schemes and spatial discretization schemes in a f-plane should be conducted to guarantee both stability and energy conservation for the 2-time stepping Euler forward Runge-Kutta (and Adams-Bashforth) schemes. It is noted that the flux form of the advection terms is energy conserving in the continuum medium, which is not guaranteed in a discretization medium, which depends on both time and spatial differencing schemes used (see Appendix D).
We have discussed the phase errors caused by different time integration schemes to represent the pure Coriolis terms in "Proposed phase error correction". However, the spatial discretization scheme of a model may also cause phase error of inertial motion. For example, A model with Arakawa A-and B-grids preserves the nondispersive property of the inertial frequency (i.e., true inertial frequency, f, is independent of wavenumber or grid size), while the C-grid scheme degrades the nondispersive frequency in continuum medium into dispersive frequency in the differencing equations, because f c = f cos 2 ( x ∆x/2) cos 2 ( y ∆y/2), which is dependent on wavenumber and grid sizes (see Wang 1996;Beckers 1999), where f c is the inertial frequency in the C-grid model, x and y are wavenumbers, and ∆x and ∆y are the grid sizes, respectively in the x and y direction. As long as ∆x and ∆y are not zero, the inertial frequency in a C-grid model is always dispersive, depending on wave number and grid size, leading to decreasing propagation speeds and dampening the amplitude of inertial waves. This is due to the fact that in the C-grid, an average of velocity u and v to the common grid is needed to calculate the Coriolis terms. As such, in a C-grid model, the high-resolution grids are strongly required in order to minimize dispersion of inertial-related wave frequency.
Phase errors in the f-plane are also affected by bottom friction and vertical and horizontal viscosity (Wang and Ikeda 1997a), and by bottom topography (such as sloping bottom; Wang and Ikeda 1997b) and geometry. Therefore, high-resolution models with accurate topography (including seamounts) and geometry (including small islands etc.) will produce relatively accurate phases for propagating waves in the f-plane.

Appendix A: Stability and phase error analysis of third-stage Euler predictor-corrector scheme
The 1 st -order Euler scheme for the pure inertial system (1) is as follows: and the 2 nd -stage predictor scheme is in a similar way:  Fig. 1, but with the centered differencing (leapfrog) scheme: Simulated lake current at 10m. No artificially amplified inertial current (i.e., unstable inertial motion) was observed on day 182 (June 2;upper panel) and day 241 (August 30; lower) in 1995 Then, the 3 rd -stage corrector scheme can be obtained by Then, we have from (8) and from (9) we have where F=fΔt. So, the corrector scheme of (10) becomes So, we have the following two equations To investigate the inertial stability of Eq. (14), following Wang (1996) and Wang and Ikeda (1997a), the von Neumann or Fourier method is used. The solution to the differencing Eq. (14) has the following Fourier wave form where V 0 =u 0 +iv 0 , =exp(i Δt), i, , Δt, and n are constant initial amplitude, eigenvalue, imaginary number, wave frequency, integration time step, and step number, respectively. Note that ω may be a complex number, i.e. ω = ω r + i ω i , |λ| = 1 when the imaginary part of wave phase is zero, i.e. ω i = 0 and λ > 1 with ω i < 0 leads to an instability with (10) u n+1 −u n Δt − f βv n+1 * * + (1 − )v n = 0, v n+1 −v n Δt + f βu n+1 * * + (1 − )u n = 0.
In summary, a pure advection scheme with 2 nd -order accuracy both in time (leapfrog) and space (centered) conserves energy, because it does not need viscosity to be stable. However, the scheme with leapfrog scheme in time and 1 st -order Euler forward (upwind) scheme in space can guarantee a stable physical solution, but not the computational mode. By contrary, the pure advection scheme with the Euler forward scheme in time is always unstable with either the Euler forward (upwind) or the centered scheme in space. In other words, the Euler forward scheme in time needs extra viscosity (diffusion) to be conditionally stable, which is the typical advection diffusion equations (system). Although the centered differencing scheme produces a computational mode in time and 2∆x short waves in space, which can be removed with certain filters (Zhang et al. 1987), its 2 nd -order accuracy and energy-conserving property (Wang 1996) possess major advantages, compared to the two-time stepping schemes.