Research on Reference State Deduction Methods of Different Dimensions

1.   Introduction

The essence of numerical weather prediction is to solve the equations that govern the movement of the atmosphere. These consist of a set of nonlinear partial differential equations containing multiple prediction variables such as temperature, pressure, wind speed, humidity, etc. Given their complexity, these equations cannot be solved analytically in a direct manner and are primarily resolved through numerical simulation using discretization methods. This study focuses on the spatial discretization method of the atmospheric motion equations.

The atmospheric motion is inherently nonlinear. The high-impact weather events that people concern are generally determined by small- and medium-scale systems overlaid on the large-scale circulation. The accumulation of seemingly minor computational errors can significantly impact the model’s predictive capabilities. One of the most typical problems is “small differences between large numbers” (Zeng, 1963). Calculating the pressure gradient force (PGF) needs the computation of the horizontal partial derivative of pressure using the difference method. Given two biggish and smooth pressures, direct calculations of the much smaller difference between them will lose significant bits and obviously reduce the accuracy of the PGF. The same is true for extended-range and seasonal weather predictions, where differences in the simulation results are determined by the disturbances superimposed on the average climate state. The accumulation of calculation errors can lead to a drift in the long-term integration results (Wu et al., 2008).

Generally, the above problems are solved by deducting a reference state, which is commonly constructed by using a temperature profile or large-scale flow field. By utilizing this reference state, the flow field is separated into a basic flow and a perturbation flow. This approach transforms the issue of dealing with “small differences between large numbers” into the direct calculation of small quantities, significantly improving the discretization accuracy of the equations and enhancing the computational stability of the model. Similarly, the issue of “drift” in climate models can be avoided or alleviated through the deduction of the reference state (Wu et al., 2008).

To use large time steps, the semi-implicit time integration method (Staniforth and Côté, 1991) is widely utilized in numerical models. Explicit integration is used for slowly changing synoptic scale waves, whereas implicit integration is used for rapidly changing gravitational waves, sound waves, and the like. By incorporating a three-dimensional (3D) reference state into the semi-implicit time integration (Su et al., 2018, 2020), the magnitude of the nonlinear term is reduced, thereby enhancing the accuracy of time integration. Additionally, in the semi-Lagrangian scheme (Staniforth and Côté, 1991), addressing numerical noise in steep terrain is crucial for improving the accuracy and stability of the model. Ritchie and Tanguay (1996) emphasized that terrain forcing remains constant during time integration and can be accounted for through the reference state. Qian et al. (2021) summarized the distinctions between full forecasting and perturbation forecasting, emphasizing that the latter can extract key information more clearly and has diverse applications in model forecasting, synoptic analysis, and diagnostic evaluation.

In summary, subtracting the reference state is crucial for enhancing the accuracy of numerical models. When the reference state is appropriately set during model integration, the flow field can be effectively segregated into basic flow and perturbation flow. This approach avoids issues such as the discretely problematic “small differences between large numbers” and thereby improves the accuracy of the dynamic core.

Over the years, model developers have conducted extensive research on the selection of reference states and the linearization of equations. Nowadays, most weather and climate models widely adopt methods such as the isothermal atmosphere or one-dimensional (1D) and two-dimensional (2D) reference state deductions. In this paper, the temperature reference state is categorized based on its dimension under geometric height coordinates. For the zero-dimensional (0D), 1D, and 2D reference state, some basic work and part of the current operational models are introduced, for the 3D and four-dimensional (4D), developments in recent years are introduced in detail.

The 0D reference state, also known as 1D reference state in some literature, refers to the isothermal atmosphere. Zeng (1963) suggested that the accuracy of the dynamic core could be improved by subtracting the reference profiles of temperature or altitude from the static equations. Early research work (Zeng, 1979; Zeng et al., 1985; Chen et al., 1987; Roeckner et al., 1996) linearized the prediction equations by subtracting the isothermal atmosphere, which is essentially a constant value of the temperature profile that does not vary with altitude. This approach is straightforward for implementation and it addresses the issue of “small differences between large numbers” to a certain extent, making it suitable for cases with low model tops. In the initial versions of the Global Regional Assimilation and Prediction System (GRAPES) model (Xue and Chen, 2008; Shen et al., 2017), the concept of the isothermal atmosphere was employed. To maintain accuracy at lower altitudes, the isothermal atmosphere was typically set to approximately 288 K.

For the 1D reference state, the temperature profile only varies with altitude. Zeng et al. (1989) and Zhang et al. (1990) constructed a 1D reference profile of temperature and pressure based on the standard atmospheric temperature stratification, using the hydrostatic balance relation. Chen and Shu (1994) subtracted the 1D reference profile in the ECMWF IFS (Integrated Forecast System), effectively reducing the discretization error, following the approach outlined by Zeng et al. (1989). Zhang (1990) and Liang (1996) adopted a 1D temperature reference profile from standard stratification in Institute of Atmospheric Physics (IAP) atmospheric general circulation model. Both Wang et al. (2004) and Wu et al. (2008) utilized a similar approach in IAP GAMIL (Grid Atmospheric Model IAP/LASG) model and in CAM3 (Community Atmosphere Model version 3), respectively. All of the fifth-generation Mesoscale Model (MM5) (Dudhia, 1993), the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019), the Prediction Across Scales (MPAS) model (Skamarock et al., 2012), and the ICOsahedral Non-hydrostatic (ICON) model (Prill et al., 2019) used a similar methodology to solve the dynamic equations by subtracting a constant profile as a function of height. Lafore et al. (1998) deducted the 1D reference profile of potential temperature and dimensionless pressure by averaging the initial field in the Meso Non-Hydrostatic (Meso-NH) model. The profile varies depending on the starting time, partially accounting for tempo-ral changes in the profile and effectively reducing the discretization error in the model. Zhang et al. (2017) also employed a similar approach in the GRAPES_MESO model to enhance its typhoon prediction capabilities.

For the 2D reference state, the temperature state is a 2D plane varying with altitude and latitude. Based on the method proposed by Zhang et al. (1990), Sheng et al. (1992) successfully integrated a 2D reference state into the global spectral model of the Australian Bureau of Meteorology Research Centre, effectively improving the prediction skill at the upper level. Wu et al. (2008) further developed the CAM3 model by deducting the topography-related surface pressure component from the surface pressure prediction equation, thus effectively improving the spatial discretization accuracy. In GRAPES_GFS, a 3D reference state (to be introduced later) is firstly implemented (Su et al., 2018, 2020); and when the horizontal resolution was increased from 0.25 to 0.125 degree, for the sake of integration efficiency and stability, Zhang et al. (2023) simplified the reference state from 3D to 2D.

For the 3D reference state, the temperature state is a 3D field that varies with altitude, longitude, and latitude. The IFS model of ECMWF utilizes surface pressure and temperature as prediction variables. Temperton et al. (2001) deducted the topography-related surface pressure component from the mass equation, and constructed a 3D temperature state based on this component at vertical levels under the pressure-based terrain following coordinate, then subtracted it from the energy equation, which effectively improved the accuracy of the PGF and Semi-Lagrangian (SL) advection at steep terrain. Note that the 3D reference state deducted by this method is not close to the real atmosphere, in terms of pressure intensity; it is less than half of the real atmosphere. The novel Finite-Volume Module of ECMWF, namely the IFS-FVM (Kühnlein et al., 2019) integrates the fully compressible non-hydrostatic equations using semi-implicit time stepping, and the perturbations of potential temperature and Exner pressure correspond to deviations from an ambient state, which satisfies a general balanced subset of control equations; for example, a stationary state in thermal wind balance, or simple stationary vertical profiles in hydrostatic balance. Based on the GRAPES_MESO, the Tropical Regional Assimilation Model for the South China Sea (TRAMS) model developed by the Guangzhou Institute of Tropical Marine Meteorology of CMA incorporates a unique 3D reference state scheme introduced by Chen et al. (2016). This 3D reference state is derived from the initial field’s smoothness, satisfying hydrostatic balance and remaining constant throughout the integration. Su et al. (2018) further improved the reference state upon Chen et al. (2016) in GRAPES_GFS by shifting the horizontal partial derivative of the reference state from the nonlinear term to the linear term. This ensures that the magnitude of the linear term is significantly larger than the nonlinear term, and the SISL (Semi-implicit SL) solver is re-derived and implemented accordingly. The 3D reference state now adopts a climate field (Su et al., 2020), maintaining hydrostatic balance and consistency throughout the integration. This improved 3D reference state scheme is more rational, significantly enhancing the accuracy and stability of the dynamic core (Shen et al., 2023). Similarly, Xu et al. (2020) made improvements to the scheme, following Su et al. (2018) and transferring the horizontal partial derivative of the reference state from the nonlinear term to the linear term in the TRAMS model. Chen et al. (2021) further refined the scheme by separating certain components from the perturbation flow, based on the hydrostatic balance at each time step, and incorporating them into the SISL solver. This new scheme effectively reduces biases in the near-surface strata. The dynamic core built on the multi-moment constrained finite-volume (MCV) method, independently developed by CMA Earth System Modeling and Prediction Centre (CEMC) (Li et al., 2013; Qin et al., 2019), utilizes a cubed-sphere grid, a local high-order reconstruction algorithm, and an implicit–explicit (IMEX) Runge–Kutta time integration scheme. A key aspect of its design involved deduction of the 3D reference state based on climate fields in the discretization process, effectively improving the accuracy and stability of the model (Chen et al., 2023).

For the 4D reference state, the temperature state is a 4D field that varies with altitude, longitude, latitude, and time. The reference state of Met Office’s Even Newer Dynamics for General atmospheric modelling environment (ENDGame) model can be considered as a quasi-4D one (Wood et al., 2014). At each time step, the potential temperature from the previous step is adjusted to achieve hydrostatic balance, serving as the reference state, and the perturbation state is the component of non-hydrostatic balance. This approach is equivalent to using a 3D reference state that evolves with time, ensuring that the reference state remains close with the current model atmosphere throughout the integration process. In the algorithm, the vertical variation of the 3D reference state is deducted from both sides of the complete equation, while the horizontal variation term is ignored. For resolutions above the kilometer scale, the vertical change in the reference state is much greater than the horizontal change. This algorithm not only guarantees precision but also simplifies the solver’s complexity, reduces the implementation difficulty, and ensures the computational stability of the model. For Finite Volume 3 (FV3) of NCEP (Chen et al., 2013), both the vertical Lagrangian coordinate and vertical Eulerian coordinate versions utilized a hydrostatic pressure component as the background when solving the vertical momentum equation, with only the perturbation of pressure being prognostic. This method improves the stability and accuracy of the system. When the hydrostatic pressure state is a 3D field and evolves with time, it can also be classified as a 4D reference state. Smolarkiewicz et al. (2019) proposed a generalized perturbation form for atmospheric dynamics across all scales; the proposed semi-implicit solver showed the potential to accurately separate the background state from the global atmospheric circulation; and they verified the usability through ideal experiments, providing a theoretical basis for further work. Qian et al. (2021) and Qian and Du (2022) separated the atmospheric equations into climatic equations and perturbation equations by incorporating time-varying reference states of all variables. This method provides more options for model design, but it also faces two fundamental challenges: the need for high-resolution temporal and spatial reanalysis data, as well as the requirement for physical process tendency data in climatic states.

Based on the SISL dynamic core of GRAPES_GFS (Shen et al., 2023; Zhang et al., 2023), this article first introduces the implementation process of 0D, 1D, 2D, 3D, and 4D reference states. Then, it compares the effects of each dimension of reference state through several benchmark tests, followed by introducing the problems encountered with the use of the 4D algorithm in the real atmospheric simulation of GRAPES_GFS. Finally, conclusions are drawn and suggestions are offered.

2.   Implementation of reference states of different dimensions

Based on the GRAPES_GFS dynamic core, the following shows the process of solving the equations by deducting the 0D, 1D, 2D, 3D, and 4D reference state, respectively. The dynamic core is based on a spherical latitude and longitude grid, shallow atmospheric approximation, and non-hydrostatic balance equations:

where u, $v$, and $w$ are the 3D wind components, $\theta\mathrm{_v}$ is virtual potential temperature, $\pi$ is Exner pressure [$\pi =$$C_p \left( {p}/p_0 \right)^{\gamma}$], f is the Coriolis force constant, $g$ is gravitational acceleration, $C\mathrm{_{\mathit{P}}}$ is specific heat capacity at constant pressure, $\gamma=C_{\mathrm{\mathit{P}}}R^{-1}$, R is gas constant, $\delta\mathrm{_{NH}}$ is the switch for hydrostatic balance, ${D_3}$ is 3D divergence, and $F_\theta ^*$ is source and sink of heat.

According to the SISL algorithm, a reference state is introduced to linearize the original equations [Eqs. (1)–(6)] into a set of linearized equations. These linearized equations are then converted into a prediction equation based on the time integration scheme. The following section takes the 3D reference state as an example to show the detail of the solver, and then describe the differences between 0D, 1D, 2D, and 4D reference state.

2.1   3D reference state

Under geometric height coordinates, the 3D reference state is:

where $\tilde \pi$ and $\tilde \theta$ represent the reference state of $\pi$ and $\theta$ under the hydrostatic balance, while ${\pi '}$ and ${\theta '}$ are perturbation components.

For $u$, convert Eq. (1) into height-based terrain following coordinate system, and then use Eq. (7) to expand the PGF term,

where the expression for coordinate transformation terms are $\hat{z}=Z_{\mathrm{T}}\dfrac{z-Z\mathrm{_s}}{Z\mathrm{_T}-Z_{\mathrm{s}}}$, $\varphi_{\mathrm{s}x}=\dfrac{1}{a\cos\phi}\dfrac{\partial Z\mathrm{_s}}{\partial\lambda}$, $\varphi_{\mathrm{s}y}=\dfrac{1}{a}\dfrac{\partial Z_{\mathrm{s}}}{\partial\phi}$, $\Delta Z\mathrm{_s}=$$Z\mathrm{_T}-Z\mathrm{_s}$, $\Delta Z{_z}=Z\mathrm{_T}-z$, and $\Delta Z_{\hat{z}}=Z_{\mathrm{T}}-\hat{z}$, $Z_{\mathrm{s}t}=\dfrac{Z\mathrm{_T}}{\Delta Z\mathrm{_s}}$, $Z\mathrm{_T}$ is model top, $Z_{\mathrm{s}}$ is terrain height. This term is further organized into linear terms (L) and nonlinear terms (N):

For $v$, similar to the $u$ equation, we can obtain

For $w$, substitute Eq. (7) of reference state into Eq. (3):

Utilize the hydrostatic balance of reference state,

and then convert Eq. (11) into height-based terrain following coordinate system,

By sorting out Eq. (13), we can obtain:

For $\pi$, substituting Eq. (7) of reference state into $\dfrac{\mathrm{d}\pi}{\mathrm{d}t}$:

The reference state does not change with time, so the second term on the right hand side of Eq. (15), i.e., ${\partial\tilde{\pi}}/{\partial t},$ equals 0. We may convert the last three terms on the right hand side of Eq. (15) to the terrain following coordinates as follows:

where $\hat w$ is the vertical velocity of the terrain following coordinate. It can be proved by the relation of the coordinate transformation that:

By substituting Eqs. (16)–(18) into Eq. (15), and simplifing them with Eqs. (19)–(20), we can obtain

For the first term on the right hand side of Eq. (4):

where ${D_3}|_{\hat{z}}=\left[\dfrac{1}{a\cos\phi}\dfrac{\partial u}{\partial\lambda}+\dfrac{1}{a\cos\phi}\dfrac{\partial(\cos\phi v)}{\partial\phi}+\dfrac{\partial\hat{w}}{\partial\hat{z}}\right]_{ \hat{z}}$.

Substitute Eqs. (21)–(22) into Eq. (4) and organize them into linear and nonlinear terms:

For $\theta$, Similar to $\pi$, we can get

Substitute Eq. (24) into Eq. (5) and organize them into linear and nonlinear terms:

The above linearized equations [Eqs. (9), (10), (14), (23), (25)] are discretized in time according to the predictor–corrector SISL scheme (Zhang et al., 2022, Shen et al., 2023), and then transformed into prediction equation of $u$, $v$, $w$, and ${\theta '}$ used in ${\pi '}$(Su et al., 2018). For the 3D reference state, the derivation process of the prediction equation becomes complicated, because the horizontal partial derivative of the reference state is classified into linear term. Next, the Preconditioned Classical Stiefel Iteration (PCSI) solver (Zhang et al., 2023) is used to solve the Helmholtz equation and obtain ${\pi '}$. Finally, the other prediction variables are inverted. The advection of water vapor is calculated separately by using the Piece-wise Rational Method (PRM) scheme (Su et al., 2013).

In addition, the method of setting the reference state is very important, and the 3D reference state can be obtained from the initial field by filtering or smoothing, or it can be set to climatic states from reanalysis data (Su et al., 2020). The 3D reference state needs to satisfy the hydrostatic balance between $\tilde \pi$ and $\tilde \theta$, and ensure the monotonicity of $\tilde \theta$ in vertical direction to maintain the stability of the model (Wood et al., 2014).

2.2   0D reference state (isothermal atmosphere)

For the separation of the prediction variables $\pi$ and $\theta$, the 0D reference state is:

The setting of the reference state is given according to the definition of $\pi$ and $\theta$ combined with the hydrostatic balance relationship:

${T_0}$ is generally selected as 288 or 300°C. For the 0D reference state, $\tilde \pi$ and $\tilde \theta$ are both functions of altitude $z$.

For the linearization of governing equation at terrain following coordinates, Eq. (9) is deformed in the 0D reference state:

In contrast to Eq. (9), because the reference state is categorized based on its dimension under geometric height coordinates, the zonal partial derivative of $\tilde \pi$ at the geometric height coordinates is equal to 0.

For linearization of $\pi$, the $\dfrac{{\partial \tilde {\pi}}}{{a\cos \phi \partial \lambda }}$ and $\dfrac{{\partial \tilde {\pi}}}{{a\partial \phi }}$ in ${L_\pi }$ of Eq. (23) are the horizontal derivatives under the terrain following plane, so its form does not change in 0D, 1D, 2D, and 4D. For example, in the case of isothermal atmosphere, the horizontal derivative of $\tilde \pi$ is 0 under geometric height plane, and not 0 under terrain following plane at the location with terrain.

2.3   1D reference state

The separation of 1D reference state from the prediction variables is the same as that of 0D reference state, but $\tilde \pi$ and $\tilde \theta$ are a function of altitude [Eq. (26)], and the linearization equation is also the same [Eq. (28)]. The vertical profile of the 1D reference state can be obtained by averaging the initial field at the geometric height plane (Lafore et al., 1998; Zhang et al., 2017). Then, $\tilde \pi$ and $\tilde \theta$ on the terrain following plane are derived according to the hydrostatic balance, and they can also be derived from the hydrostatic balance based on the average temperature decline rate of atmospheric stratification (Dudhia, 1993; Chen and Shu, 1994; Prill et al., 2019; Skamarock et al., 2019).

2.4   2D reference state

For the 2D reference state, $\tilde \pi$ and $\tilde \theta$ are a function of both latitude and altitude [see Eq. (26)]. The linearization equation for $u$ is the same as Eq. (28), while that for $v$ retains the meridional partial derivative of the reference state on the geometric height plane, as Eq. (10). The 2D reference state can be obtained by averaging the initial field in the zonal direction on the geometric height plane, or by averaging the climatic states from reanalysis data in the zonal direction (Zhang et al., 2023). Then, $\tilde \pi$ and $\tilde \theta$ on the terrain following plane are derived according to the hydrostatic balance.

2.5   4D reference state

In the 4D reference state, both $\tilde \pi$ and $\tilde \theta$ are a function of longitude, latitude, altitude, and time. In the algorithm, the reference states deducted at each time step are independent of each other and are not passed into the next time step. Thus, there is no need to incorporate the time variation term of the reference state into the equations. One of the key problems of 4D algorithm is the adjusting of the frequency of the 3D reference state, which can be adjusted at every step or at a certain time interval, and theoretically, adjusting at every step is the best option. In the subsequent baroclinic benchmark test, we tested the difference between adjusting at each step and adjusting at one-day intervals. In the SISL solver, after adjusting the reference state, it is necessary to recalculate the linear term and the A-matrix of the Helmholtz equation, which will increase the computational cost at a certain extent. The other part of the solver remains the same as the 3D reference state. The 4D reference state can be obtained by filtering or smoothing the current state or by employing curved surface fitting methods based on the current state, it is necessary to ensure the hydrostatic balance relationship of $\tilde \pi$ and $\tilde \theta$, as well as monotonicity of $\tilde \theta$ in vertical direction.

3.   Benchmark test

3.1   Rossby–Haurwitz wave

This test mainly simulates the propagation of weather-scale waves and can demonstrate the dissipation, mass conservation, and phase error of the dynamic core. The specific test setup is referred to Jablonowski et al. (2008), the horizontal resolution is 0.5 degrees, 60 evenly spaced vertical layers are used, the layer height is 400 m, the time step is 900 s, and the integration period is 15 days.

In the test, the 0D reference state is the isothermal atmosphere of 300°C, which is given as $\tilde \pi$ (Fig. 1a) and $\tilde \theta$ according to Eq. (12); the 1D reference state is the horizontal average profile of the initial field (Fig. 1b); the 2D reference state is the 2D distribution obtained after averaging the initial field in the zonal direction (Fig. 1c); the 3D reference state directly sets the initial field of $\pi$ and $\theta$ to $\tilde \pi$ (Fig. 1d) and $\tilde \theta$, with the initial perturbation ${\pi '}$ and ${\theta '}$ being zero; for the 4D reference state, at the beginning of each time step during the model integration, a 5-point second-order precision Shapiro filter (Shapiro, 1975) is alternately applied 10 times for $\pi$ and $\theta$ in both the zonal and meridional directions. The filtered variables are then taken as $\tilde \pi$ and $\tilde \theta$, then $\tilde \pi$ on the ground and $\tilde \theta$ is used to recalculate $\tilde \pi$ (Fig. 1e) above the ground according to hydrostatic balance.

Download

Larger image

PowerPoint slide

Fig  1.  Exner pressure ($\tilde \pi$) at 200-m model layer after 15 days integration for (a) 0D, (b) 1D, (c) 2D, (d) 3D, and (e) 4D reference states.

Ideally, the Rossby waves should not decay as they orbit the earth from west to east (Jablonowski et al., 2008, Wan, 2009). Compared with the surface pressure at the initial time (Fig. 2a), after 15 days of integration, the dissipation is more obvious in 0D case, with the Rossby wave around the equator significantly deformed; the dissipation is reduced in 1D and 2D; and is further weakened in 3D and 4D. Among them, the accuracy of 4D (Fig. 2f) is the highest and its wave pattern is most similar to that in the initial field (Fig. 2a), which indicates that gradually adjusting the reference state closer to the current state can reduce the perturbation during the model integration, improving the spatial discretization accuracy and reducing the numerical dissipation.

Download

Larger image

PowerPoint slide

Fig  2.  (a) Initial surface pressure and surface pressures after 15 days of integration for (b) 0D, (c) 1D, (d) 2D, (e) 3D, and (f) 4D reference states.

3.2   The accuracy of PGF over steep terrain

The accuracy of PGF over steep terrain in the stationary atmosphere was investigated by using the benchmark test No.2_0 in the Dynamical Core Model Intercomparison Project (DCMIP) (Ullrich et al., 2012). A circular, Schar-like mountain was constructed in the tropics, centered at 0°N, 150°W, with a radius of 135° and a maximum height of 2000 m (Fig. 3). The simulation utilized a horizontal resolution of 1.0° and 30 vertical layers, each 400 m high. It ran for 6 days with a 600-s time step, initializing from a steady-state atmosphere. Detailed parameters are provided in Ullrich et al. (2012). Although the test simulated a stationary atmosphere, the accumulation of discrete errors also altered the distribution of temperature and pressure, so we also conducted a 4D reference state test. The specific 0–4D reference state setting was identical to that used in the Rossby–Haurwitz wave test.

Ideally, the model atmosphere should remain stationary as its initial state, but for a discretized numerical model, the coordinate transformation at steep terrain and the deviation of the numerical discretization will alter the pressure gradient force, which in turn will trigger a certain movement. The smaller this motion is, the smaller the discretization error will be. The performance of the models in this test varied significantly (Ullrich et al., 2012). Compared the results, the discretization error is largest for the 0D reference state (Fig. 3b), and was reduced in 1D (Fig. 3c) and 2D (Fig. 3d), and further in 3D (Fig. 3e). Compared to 3D, the magnitude of error in 4D (Fig. 3f) does not change significantly, but its distribution becomes more random. The correlation between the error and the terrain has weakened, indicating the effect of deducting the basic flow. This test demonstrates that the spatial discretization error of the model is significant in steep terrain, and the accuracy can be effectively improved by constructing a high-precision reference state and deducting the influence of complex terrain.

Download

Larger image

PowerPoint slide

Fig  3.  (a) Vertical profile of the terrain height at the equator. Vertical profile of zonal wind at the equator after 6 days of integration for (b) 0D, (c) 1D, (d) 2D, (e) 3D, and (f) 4D reference states.

3.3   Mountain-induced Rossby wave

This test aims to simulate the generation, propagation, and extinction of the Rossby waves influenced by large-scale topography. A bell-shaped mountain was placed at the mid-latitudes of the Northern hemisphere, centered at 30°N, 90°E. The half-width and maximum height of the mountain are 150 km and 2000 m. The background atmosphere was set to a constant westerly wind speed of 20 m s^–1. The simulation used a horizontal resolution of 0.5 degrees, with 60 even vertical layers, each 400 m high, ran for 15 days with a 900-s time step. Detailed parameters are provided in Jablonowski et al. (2008). The 0–4D reference state setting is identical to that used in the Rossby–Haurwitz wave test.

The mountain range prompts a specific vortex motion within the zonal wind field, which gradually transforms into a Rossby wave train propagating downstream from the mountain. Although no analytic solution exists for this test, our simulations are able to replicate the Rossby-wave pattern compared to Jablonowski et al. (2008) and Wan (2009). For different dimensions of the reference state, the intensity of the trough and ridge for 1D (Fig. 4c) and 2D (Fig. 4e) shows no significant change compared to 0D (Fig. 4a). For 3D (Fig. 4g) and 4D (Fig. 4i), the intensity of the system is increased, and the gradient of the height is obviously strengthened. We further assess accuracy by examining the relative vorticity, calculated as $\dfrac{{\partial v}}{{\partial x}} - \dfrac{{\partial u}}{{\partial y}}$. Similar to the conclusion regarding the height field, 1D (Fig. 4d) and 2D (Fig. 4f) did not change much, while 3D (Fig. 4h) and 4D (Fig. 4j) showed system enhancement, with 4D being slightly stronger than 3D. This indicates that spatial discretization accuracy can be effectively enhanced by subtracting the reference state of higher dimensions, which is beneficial for simulating the emergence and evolution of synoptic-scale systems.

Download

Larger image

PowerPoint slide

Fig  4.  The 700-hPa geopotential height field (left; gpm, denoted by “m” on the color bar) and relative vorticity (right; 10^–5 s^–1) after 15 days of integration. (a, b) The 0D reference state, (c, d) 1D minus 0D, (e, f) 2D minus 0D, (g, h) 3D minus 0D, and (i, j) 4D minus 0D.

3.4   Baroclinic wave

The baroclinic wave is triggered when using the steady-state initial conditions with the overlaid zonal wind perturbation, the model can capture the wave’s occurrence, development, and extinction (Jablonowski et al., 2008). The simulation, using a 0.5-degree horizontal resolution, with 60 even vertical layers, each 400 m high, ran for 15 days with a 900-s time step. Detailed parameters are provided in Jablonowski et al. (2008). Because the initial field of $\pi$ and $\theta$ are uniform in meridional direction, the 2D and 3D reference states are the same, and thus only 0D, 1D, 3D, and 4D are compared. The specific setting of the reference state is identical to that used in the Rossby–Haurwitz wave test. Since this test is most sensitive to the 4D reference state, we tested different reference state adjustment intervals: adjusting each day (4D_day) and adjusting every step (4D_step).

No analytic solution is available for this test. When comparing the numerical result with Jablonowski et al. (2008), our tests accurately reproduce the occurrence and propagation of the baroclinic wave. For different dimensions of the reference state, the wave trains are of the same phase, the intensity of baroclinic wave shows a similar but not significant enhancement for 1D (Fig. 5c) and 3D (Fig. 5e) compared to 0D (Fig. 5a), the enhancement for 4D_day (Fig. 5g) is more obvious and further increased for 4D_step (Fig. 5i). In terms of the change in relative vorticity, it also shows similar characteristics as the temperature field, the relative vorticity is streng-thened in 1D (Fig. 5d) and 3D (Fig. 5f), and further strengthened in 4D_day (Fig. 5h), with significantly strengthening in 4D_step (Fig. 5j). The 3D and 1D results change little, but the 4D results change significantly, the reason for this is that the temperature and pressure of initial field in this test are uniform in the zonal direction, and the 3D reference state of $\pi$ and $\theta$ does not contain perturbations, so its results are similar to 1D. However, the 4D reference state can always be close to the real wave during the development process of baroclinic waves. This test demonstrates that when the reference state is close to the real state, it can effectively improve the simulation accuracy of the model.

Download

Larger image

PowerPoint slide

Fig  5.  The 850-hPa temperature field (left; °C) and relative vorticity (right; 10^–5 s^–1) after 9 days of integration. (a, b) The 0D reference state, (c, d) 1D minus 0D, (e, f) 3D minus 0D, (g, h) 4D_day minus 0D, and (i, j) 4D_step minus 0D.

3.5   Comparison of mass conservation of the benchmark test

The above tests are based on the GRAPES_GFS dynamic core, which adopts the SISL scheme. Because the point mapping algorithm from the upstream point to the arrival point cannot ensure the mass conservation, a correction algorithm (Su et al., 2016) to ensure the conservation of dry air mass is used. Here, the correction operator of mass conservation is turned off, and the effects of different dimensional reference state on the mass conservation in three sets of benchmark tests (Rossby–Haurwitz wave, Baroclinic wave, and Mountain-induced Rossby wave) are compared.

The total mass of the atmosphere is represented by the global mean surface air pressure, with initial values slightly varying in three tests, for ease of comparison, these values are scaled to a unified 1000 hPa. As shown in Fig. 6, the 0D reference state exhibits the most significant mass change, with a 0.8-hPa increase in the Rossby–Haurwitz wave test, a 0.3-hPa increase in the mountain-induced Rossby wave test, and a decrease of 0.3 hPa in the Baroclinic wave test. Adoption of the high-dimensional reference state enhances spatial discretization accuracy and significantly improves mass conservation, resulting in global mean sea level pressure fluctuations within 0.2 hPa across multiple tests. In terms of dimensionality, after 15 days of integration, the average absolute change in global mean sea level pressure across the three tests was 0.46 hPa for 0D, 0.13 hPa for 1D, 0.14 hPa for 2D, 0.10 hPa for 3D, and 0.11 hPa for 4D. Overall, 3D and 4D exhibit superior mass conservation compared to 1D and 2D, and both significantly outperform the 0D reference state.

Download

Larger image

PowerPoint slide

Fig  6.  Change of global mean surface air pressure during the integration. Baroclinic wave (abbreviated as B; solid line), Rossby–Haurwitz wave (abbreviated as H; long dashed line), and mountain-induced Rossby wave (abbreviated as M; dotted line). 0D (black), 1D (red), 2D (pink), 3D (blue), and 4D (green). The subplot is the mass change at Day 15.

4.   Real atmospheric simulation

In the real atmospheric simulation, Zhang et al. (2017) achieved 1D reference state in the GRAPES_MESO model, while Su et al. (2018, 2020) achieved 3D reference state in the GRAPES_GFS model, Zhang et al. (2023) simplified the reference state from 3D to 2D in GRAPES_GFS for efficiency and stability, and these efforts yielded the desired results. However, when attempting to implement the 4D algorithm in the GRAPES_GFS model, the model frequently exhibited integral overflow. Analyses revealed that the reference state changes over time closely mirror the actual atmosphere, encompassing information across various scales. The horizontal partial derivative of the reference state exhibits significant numerical jitters in the vicinity of steep terrain, especially at high resolution. Although authentic, these jitters behave as numerical noise within the dynamic solver. In addition, GRAPES_GFS uses latitude and longitude grids, this kinds of numerical noise is more pronounced at higher latitudes. In the process of solving the dynamical equations, the horizontal derivative of the reference state is very important, it is classified into the linear term and participates in the discrete process of SISL. This kind of numerical noise is introduced into the A-matrix, making the Helmholtz equation difficult to converge and often leading to integration overflow. The UKMO encoun-tered similar difficulties with the implementation of the 4D reference state for the ENDGame model (Wood et al., 2014). To circumvent the issue, they deducted only the vertical change of the 3D reference state on both sides of the full equation, ignoring the horizontal change. This approach simplifies the scheme’s implementation and ensures computational stability. Due to the different discrete methods used in the equations for the GRA-PES_GFS and ENDGame, we cannot adopt a similar approach. We also tested reducing the update frequency of the 3D reference, not updating it every step but once a day, which did not solve the stability problem. Because it is a daily variation of the reference state, the 3D reference state is still constructed based on the forecast field at a certain time step, which dose not address the problem of smooth derivatives, even if we add a certain number of filters.

5.   Conclusions and suggestions

Through the review and benchmark test, the paper draws the following conclusions: (1) The atmospheric motion is inherently nonlinear, when solving the control equations, it can be linearized by deducting the reference state. This separation of the basic flow and perturbation flow effectively addresses the issue of “small differences between large numbers” encountered when solving terms such as PGF, thereby enhancing the spatial discretization accuracy of the dynamic core. (2) The high-dimensional reference state closely approximates the real atmosphere across different longitudes, latitudes, altitudes, and time, leading to a more comprehensive and effective improvement in spatial discretization accuracy.

Based on the test results and various problems encountered in the research, the following suggestions are given for the deduction method of the reference state:

(1) After deducting the reference state, the model prediction variables transition from full variables to perturbation quantities, the reference state is calculated beforehand, and any calculation error therein becomes a constant error source during integration, necessitating cautious handling. Typically, the reference state must satisfy hydrostatic balance in the vertical direction. If the predicted variable is potential temperature, then its reference state must satisfy monotonicity in the vertical direction, a statically stable profile will ensure the model’s stability. If the height is derived from temperature, the vertical derivative of temperature is determined by using the difference method, while the vertical derivative of height is provided analytically, and vice versa. For the 3D and 4D reference states, the horizontal derivative is commonly computed by using a second-order horizontal difference method. However, the use of high-precision difference algorithms may introduce jitters near steep terrain, creating a persistent source of computational noise in the integration process, potentially leading to integration overflow. It is therefore essential to strike a balance between calculation accuracy and integration stability.

(2) Regarding the 3D reference state, some designers aspire to eliminate the terrain’s influence from it. However, the horizontal derivative of the reference state included in the equation may result in significant jitter, potentially compromising the model’s stability. Alternatively, some designers opt to subtract large-scale stationary motion by adjusting the reference state analogously to the climate field, enabling the model to solely predict small-scale transient motion.

(3) When utilizing explicit time integration schemes, such as Runge–Kutta, the reference state can be straightforwardly subtracted. Different dimensions have minimal impact on the equation solver, resulting in a relatively simple implementation. Apart from subtracting the reference state of pressure and temperature, you might also explore subtracting the reference state of horizontal wind. Conversely, for models employing implicit time integration schemes, the 2D, 3D, and 4D reference states necessitate the consideration of horizontal variations in the linear term. Consequently, the equation solver must be re-derived whenever the linear term changes, leading to an increase in implementation complexity.

(4) The 4D reference state changes with time, and the reference state can be adjusted to remain closely aligned with the real atmosphere during the integration process. Theoretically, it is best to update the reference state at each step, decreasing the update frequency of the reference state will gradually lead to a simulation that is closer to 3D reference state, and we also verified this through benchmark test. The 4D algorithm appears to be the most optimal choice, but the question of how to achieve this in simulating the real atmosphere is quite specific. Currently, the ENDGame of Met Office has successfully implemented a quasi-4D reference state in its operational model, while the FV3 of NCEP deducts the hydrostatic pressure component when solving the vertical momentum equation, which is also equivalent to deducting the 4D reference state.

(5) In climate models or global medium-range weather models, it is advisable to avoid using a simple 0D reference state (i.e., an isothermal atmosphere). Instead, the use of 2D, 3D, or even 4D reference states is recommended to effectively separate the basic flow from the perturbation flow. This approach ensures high calculation accuracy across different altitudes and latitudes. In the context of high-resolution limited area models, a vertically varying 1D reference profile can significantly enhance spatial discretization accuracy compared to an isothermal atmosphere. Given the smaller computational area, the adoption of 2D and 3D reference states does not yield a noticeable improvement in prediction accuracy over the 1D reference profile.