the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Towards the closure of momentum budget analyses in the WRF (v3.8.1) model
TingChen Chen
ManKong Yau
Daniel J. Kirshbaum
Budget analysis of a tendency equation is widely utilized in numerical studies to quantify different physical processes in a simulated system. While such analysis is often postprocessed when the output is made available, it is well acknowledged that the closure of a budget is difficult to achieve without temporal and/or spatial averaging. Nevertheless, the development of errors in such calculations has not been systematically investigated. In this study, an inline budget retrieval method is first developed in the WRF v3.8.1 model and tested on a 2D idealized slantwise convection case with a focus on the momentum equations. This method extracts all the budget terms following the model solver, which gives a high accuracy, with a residual term always less than 0.1 % of the tendency term. Then, taking the inline values as truth, several offline budget analyses with different commonly used simplifications are performed to investigate how they may affect the accuracy of the estimation of individual terms and the resultant residual. These assumptions include using a lowerorder advection operator than the one used in the model, neglecting grid staggering, or following a mathematically equivalent but transformed format of the governing equations. Errors in these postprocessed analyses are found mostly over the area where the dynamics are the most active, thus impairing the subsequent physical interpretation. A maximum 99th percentile residual can reach >50 % of the concurrent tendency term, indicating the danger of neglecting the residual term as done in many budget studies. This work provides general guidance not only for budget diagnoses with the WRF model but also for minimizing the errors in postprocessed budget calculations.
The atmosphere is a complex system with different scales of motion. Its dynamics are governed by a set of fluid equations based on the fundamental laws of physics. Although the equation set cannot be solved analytically, numerical models can be used to simulate the observed weather and climate systems to improve our understanding of the atmosphere. Due to the complexity and nonlinearity of the numerical models, budget analysis is often employed to interpret the results by quantifying the contribution of each term (i.e., physical process) in a tendency equation that governs the evolution of a certain quantity in the simulated system. The accuracy of a given budget analysis can be estimated from the residual term, defined as the difference between the tendency term on the lefthand side (lhs) of the equation and the summation of all the forcing terms on its righthand side (rhs). Budget analysis has been performed on diverse properties (e.g., momentum, temperature, water vapor, vorticity) of many systems on various scales, including the Madden–Julian oscillation (MJO; e.g., Kiranmayi and Maloney, 2011; Andersen and Kuang, 2012), tropical cyclones (e.g., Zhang et al., 2000; RiosBerrios et al., 2016; Huang et al., 2018), squall lines (e.g., Sanders and Emanuel, 1977; Gallus and Johnson, 1992; Trier et al., 1998), supercell thunderstorms (e.g., Lilly and Jewett, 1990), and so on.
Despite the popularity of the budget analysis, it is generally acknowledged that, in model postprocessing analysis, obtaining a closed budget with a negligible residual is difficult (e.g., Kanamitsu and Saha, 1996) and has been accomplished mostly in time or domainaveraged budget calculations (e.g., Lilly and Jewett, 1990; Balasubramanian and Yau, 1994; Arnault et al., 2016; Kirshbaum et al., 2018; Duran and Molinari, 2019). Even in the case of averaged budgets, the residual term that contains nonexplicitly diagnosed physics can be larger than the tendency term (e.g., Liu et al., 2016), and many studies simply do not display the residual, making the proper interpretation of the budget analysis difficult.
The “residual analysis method” is sometimes utilized to obtain an indirect estimation of the physical processes that are hard to diagnose or are unresolved in a set of analysis or observational data. In such cases, a nonnegligible residual is sometimes used to gain insight into such processes. However, as just discussed, the residual term also contains the inaccuracies associated with the calculations within the budget analysis (e.g., Kornegay and Vincent, 1976; Abarca and Montgomery, 2013). It is thus unclear whether the unresolved physics in such data sets do indeed comprise the main component of the residual without considering the contributions of other sources of errors in the budget calculation (Kuo and Anthes, 1984). Whereas it is almost impossible to separate the subgridscale, unresolved processes from other errors in reanalysis or observational data (e.g., Hodur and Fein, 1977; Lee, 1984), the focus of this study is on numerical model data where the local tendency and all the associated resolved and parameterized physics can be obtained from the model. Thus, the residual term in this study specifically refers to errors in the budget calculation.
To reduce the residual, an inline budget analysis that extracts all the terms of a prognostic equation directly from the model during its integration is generally the most accurate. However, the procedure has been reported only in a few studies (e.g., Zhang et al., 2000; Lehner, 2012; Moisseeva, 2014; Moisseeva and Steyn, 2014; Potter et al., 2018; see Appendix A for a summary and comparison among these works). Most other studies still conduct the offline or postprocessing budget analysis when the output is made available after the model integration. Some specific suggestions have been given in the past regarding how to reduce the error of postprocessed budget analysis. For example, Lilly and Jewett (1990) emphasized the importance of evaluating terms using the same differencing scheme, grid stretching, and grid staggering as that used in the simulation model. However, it is uncertain whether these rules have been widely followed, and how much of a reduction in residual can be obtained with this approach.
In some postprocessed budget analyses, transformed equations with different assumptions from those in the model are used and naturally lead to errors in the budget results. On the other hand, even when the same form of the equations is followed, errors can still arise from multiple sources during the postprocessing. Some errors are inherent in the time discretization scheme of the model, some are traced to the numerical methods in solving the temporal or spatial derivatives with finite differencing (e.g., Kuo and Anthes, 1984), and others might emerge during the interpolation or extrapolation from model grids to analysis grids (e.g., Lilly and Jewett, 1990). While the tendency term is often the result of a few cancelations among competing forcing terms, the seemingly nondominant terms may be as important as the large forcing terms in determining the sign and the value of the tendency. Thus, an incorrect estimation of even a small term may result in a residual with magnitude comparable to the tendency term, hindering the subsequent physical interpretation.
A few models, such as the Cloud Model 1 (CM1; Bryan and Fritsch, 2002) and the High Resolution Limited Area Model (HIRLAM; Undén et al., 2002), include inline budget diagnoses that users can choose to include in the model output. However, many other commonly used models (e.g., FifthGeneration NCAR/Penn State Mesoscale Model (MM5; Grell et al., 1994), Weather Research and Forecasting Model (WRF; Skamarock et al., 2008), the Advanced Regional Prediction System (ARPS; Xue et al., 2000, 2001), and the Regional Atmospheric Modeling System (RAMS; Pielke et al., 1992)) do not have this capability. In this study, we develop an inline momentum budget retrieval tool in the Advanced Research WRF model, one of the most widely used numerical weather prediction models. During the period 2011–2015, there were on average 510 peerreviewed journal publications involving WRF per year (Powers et al., 2017). Given the widespread use of WRF for both realcase and idealized modeling, such a budget tool may prove useful in numerous applications. In our budget diagnosis, each contributing term is extracted during the model integration and stored as a standard output. In so doing, we essentially solve the prognostic variables as done in the model so that the two sides of the tendency equation are always in balance regardless of the output time interval. By taking the results from the inline budget analysis as truth, we then perform several different postprocessing budget analyses with commonly made simplifications or a different format of equation. Comparisons between the postprocessed budgets and the inline/true values are made to investigate the potentially large errors in each forcing term and the resultant residuals.
2.1 Model and momentum equations
The WRF configuration used in this study is a twodimensional [(y, z); no variation in the x direction], fully compressible, nonhydrostatic, and idealized version of the Advanced Research WRF model, version 3.8.1 (Skamarock et al., 2008). Here we briefly revisit the parts that are relevant to the momentum budget analysis. The governing equations in the WRF model are cast on a terrainfollowing dryhydrostatic pressure coordinate. This vertical coordinate, η, is defined as
where p_{dh} is the hydrostatic pressure of the dry air and μ_{d} represents the mass of the dry air per unit area in the column; ${\mathit{\mu}}_{\mathrm{d}}={p}_{\mathrm{dh}\mathit{\_}\mathrm{sfc}}{p}_{\mathrm{dh}\mathit{\_}\mathrm{top}}$, where p_{dh_sfc} and p_{dh_top} indicate the values of p_{dh} at the surface and the top of the dry atmosphere, respectively.
To ensure conservation properties, the model equations are formulated in flux form, with the prognostic variables coupled with μ_{d}. The fluxform momentum components are defined as
where u, v, and w are the two horizontal and vertical velocities, respectively. Note that the drymasscoupled velocities (U, V, W) on coordinates (x, y, z) have units of pascal meter per second, and the drymasscoupled vertical velocity on η coordinate, Ω, has a unit of pascal per second. For the idealized 2D case on an f plane as in this study, the momentum equations in the WRF model are written as
where
is the fluxform advection, p is the full pressure with inclusion of vapor, ϕ is the geopotential, f is the Coriolis parameter, r_{e} is the mean earth radius, and α and α_{d} are the full and dryair specific volume, respectively. In our selected microphysics scheme (Thompson et al., 2008), six hydrometeors are included, and thus $\mathit{\alpha}={\mathit{\alpha}}_{\mathrm{d}}{\left(\mathrm{1}+{q}_{\mathrm{v}}+{q}_{\mathrm{c}}+{q}_{\mathrm{r}}+{q}_{\mathrm{i}}+{q}_{\mathrm{s}}+{q}_{\mathrm{g}}\right)}^{\mathrm{1}}$, where q_{v}, q_{c}, q_{r}, q_{i}, q_{s} and q_{g} are the mixing ratios for water vapor, cloud, rain, ice, snow, and graupel, respectively. The rhs forcing terms for the V tendency include the fluxform advection (ADV), horizontal pressure gradient force (PGF), Coriolis force (COR), vertical (earthsurface) curvature (CUV), and the remaining physics (P_{V}). For the W tendency, the rhs forcings contain the fluxform advection (ADV), net force between the vertical pressure gradient and buoyancy (PGFBUOY), curvature effect (CUV), and the remaining physics (P_{W}). The remaining physics may include diffusion, damping processes, and other parameterized physics, depending on the model setup. Note that for closing the budget analysis, all the known physics processes that come into play should be explicitly written in the equation and be diagnosed or directly retrieved from the model. The residual (res) is added on the last rhs term in Eqs. (1) and (2) to represent the imbalance between the two sides of the equation during budget analysis, but it is not part of the original equations solved in the model.
To develop an inline budget retrieval tool, it is important to understand how these prognostic variables are advanced in the WRF model. Governing equations are first recast to perturbation forms with respect to a dry hydrostatically balanced reference state that is a function of height only (defined at initialization) to reduce truncation errors and machine rounding errors. Specifically, variables of p, ϕ, α_{d}, and μ_{d} are separated into reference and perturbation components, e.g., $p(x,y,\mathit{\eta},t)=\stackrel{\mathrm{\u203e}}{p}\left(z\right)+{p}^{\prime}(x,y,\mathit{\eta},t)$. The introduction of these perturbation variables only changes the expressions for the rhs terms PGF and PGFBUOY in Eqs. (1) and (2), which will not be shown here for simplicity. Readers can refer to Skamarock et al. (2008, chap. 2.5) for more details.
Based on Skamarock et al. (2008), Fig. 1 summarizes the WRF integration strategy. The integration is wrapped by a thirdorder Runge–Kutta (RK3) scheme, in which the prognostic variables (generalized as Φ here) are advanced from t to t+Δt given their corresponding partial differential equations, $\frac{\partial \mathrm{\Phi}}{\partial t}=F\left(\mathrm{\Phi}\right)$, following a threestep strategy:
where Δt is the model integration time step and F, the largestep forcing, represents the summation of all the rhs terms of Eqs. (1) and (2) excluding the residual. Although the parameterized forcings stay fixed from step one to three as most of the parameterization schemes are called only once at the first RK3 step, the rest of the nonparameterized forcings and thus the total F are changed with the updated Φ^{*} and ${\mathrm{\Phi}}^{**}$ at the second and third RK3 step. Within each RK3 step, a subset of integration with a relatively smaller time step is embedded to accommodate highfrequency modes for numerical stability (Wicker and Skamarock, 2002; Klemp et al., 2007; Skamarock et al., 2008). A maximum number of small steps in one model integration step can be specified by the user. To improve accuracy in the temporal solver, the variables being advanced in this smallstep integration are the temporal perturbation fields, defined by the deviation from their more recent RK3 predictors: ${\mathrm{\Phi}}^{\prime \prime}=\mathrm{\Phi}{\mathrm{\Phi}}^{t*}$, where ${\mathrm{\Phi}}^{t*}={\mathrm{\Phi}}^{t}$, Φ^{*} and ${\mathrm{\Phi}}^{**}$ for the first, second, and third RK3 step, respectively. Thus, the perturbation momentum equations to be solved are driven by the largestep forcings and the smallstep (sometimes referred as “acousticstep” although it deals with both acoustic and gravity wave modes (e.g., Klemp et al., 2007; Skamarock et al., 2008)) corrections:
where τ indicates the time in the smallstep integration, and C as well as ${c}_{\mathrm{s}}^{\mathrm{2}}$ are soundwaverelated terms (Skamarock et al., 2008, chap. 3.1.2). Here we leave out the details regarding the smallstep terms that are irrelevant to the inline budget retrieval. Note that the overbar in Eq. (6) indicates a forwardintime averaging operator for the smallstep modes to damp instabilities associated with vertically propagating sound waves (see Eq. 3.19 in Skamarock et al., 2008). Equations (5) and (6) are the ones used to integrate the prognostic momentum fields in the WRF model. For each RK3 step, after the total largestep forcing F is determined, V^{′′} and W^{′′} are defined and advanced within the smallstep scheme by a loop that adds F multiplied by a time interval, Δτ (varies with different RK3 steps; see Fig. 1), and the smallstep forcing (ACOUS). After the smallstep integration loop ends, V and W are then recovered from their temporal perturbation fields and moved forward to the next RK3 step. While it is not relevant to the momentum equations discussed here, for some variables directly contributed by the microphysics scheme, the associated contribution should be considered after the RK3 integration loop ends as the microphysics are integrated externally using an additive time splitting (Fig. 1) (Skamarock et al., 2008, chap. 3.1.4).
2.2 Experimental setup
The main discussion of this study will focus on a 2D (y, z) idealized simulation of slantwise convection. This process releases conditional symmetric instability (CSI), which can be idealized by assuming no flow variations along the direction of thermal winds (denoted as the x direction in our setup; Markowski and Richardson, 2010, chap. 3.4). The initial field consists of a thermally balanced uniform westerly wind shear in x. This baroclinic environment contains no conditional (gravitational) instability, no inertial stability, and no dry symmetric instability but does contain some CSI. A twodimensional bubble containing perturbations of potential temperature and zonal wind is added to initiate convection and a slanted secondary circulation (v, w). See Appendix B for more details about the experimental setup. The domain size is 1600 and 16 km in the y and z direction, respectively, with a horizontal grid length of 10 km and 128 vertical layers. The model integration time step is 1 min. For simplicity, the only parameterization used is the Thompson microphysics scheme (Thompson et al., 2008). In addition, the upperlevel implicit Rayleigh vertical velocity damping (damp_opt = 3) is also activated (Skarmarock et al., 2008, chap. 4.4.2). The former does not directly contribute to the momentum fields (although it can affect the momentum field indirectly through density and pressure variations), and the latter, contained in P_{W} in Eqs. (2) and (6), affects only the W momentum budget. No subgrid turbulence scheme is used (diff_opt = 0). The WRF model offers different orders of advection operators, and the default third and fifthorder operators are selected for the vertical and horizontal in this case, respectively. Most of the subsequent analyses and discussion are based on this slantwise convection case with a 10 km grid length unless specified otherwise. Two other simulations, one of which uses the same setup but with an increased horizontal resolution of 2 km, will be discussed in Sect. 4.
Figure 2 shows the 48 h evolution of the 99th percentiles of v and w (hereafter the lowercase indicates that the calculation uses the uncoupled momentum field) and their tendencies. For the 10 km case, the horizontal velocity reaches its peak in about 20 h, a few hours after the vertical velocity reaches its maximum, and then undergoes a weakening. Both v and w tendencies are maximized at around 15 h. To understand the evolution of the associated flow dynamics, a momentum budget analysis serves as a natural choice. However, as a preliminary step prior to carrying out such analysis, we focus only on the technical discussion of the budget analysis methodology. The physical interpretation of the motion is beyond the current scope and will be presented in a subsequent paper.
3.1 Inline momentum budget analysis
For the inline budget analysis, all the terms are retrieved directly from the model for all the integration time steps, and therefore they represent the “instantaneous” terms that act over the specified short integration time window. For the largestep forcing, the WRF model accumulates all forcing terms at the beginning of each RK3 step. To separate them, we simply take the difference before and after WRF calls the subroutine for each largestep forcing, store their values separately, and output only the values at the third RK3 step (the total forcing is $F\left({\mathrm{\Phi}}^{**}\right)\mathrm{\Delta}t$ as shown in Fig. 1). As for the contribution of the smallstep modes, they are obtained by accumulating over all the small steps in the third RK3 step (ACOUS sum shown in Fig. 1). It is worth noting that Eq. (6) is a vertically implicit equation that couples with the geopotential tendency equation (Skarmarock et al., 2008; Klemp et al., 2007). A tridiagonal equation for the vector W (involving three grid points in the vertical direction) is thus solved (Satoh, 2002). This means that W (the scalar at a given grid point) is not advanced by linear additions in the smallstep or acoustic scheme. To ensure the closure of the inline retrieval budget, we simply take the total changes that are contributed by the implicit solver in the acoustic scheme as smallstep modes of W in the third RK3 step. Note that this way does not violate the original W equation in Eq. (6). The contribution from these accumulated smallstep modes in the V and W tendency budgets are combined with their largestep PGF and PGFBUOY, respectively, as they share the same mathematical expressions. Finally, we add the inline calculation for the tendency term outside of the RK3 integration loop, after the microphysics scheme:
where Δt is the model integration time step and Φ represents V or W (coupled momentum; hereafter the momentum tendency with capital V or W refers to the lhs term derived for the budget analysis). The values of Φ at times t and t+Δt, the latter denoted by superscripts, are termed the current and predicted states, respectively. Note that while variables of momentum tendencies (specifically named “ru_tend, rv_tend and rw_tend”) can be directly outputted from the WRF model by modifying the registry file, these variables do not necessarily represent the actual momentum changes that consider all the physical (e.g., microphysics, smallstep modes) and nonphysical processes (e.g., damping) but only the summation of all the largestep forcings.
Figures 3 and 4 present the results of the inline budget analysis for horizontal momentum and vertical momentum, respectively, at three selected times (6, 12, and 16 h). To demonstrate the momentum changes in a common physical unit (velocities; meter per second), every term of the fluxform budget equation shown in this paper is divided by the dryair mass, ${{\mathit{\mu}}_{\mathrm{d}}}^{t+\mathrm{\Delta}t}$ (so that, for example, the V tendency has a unit of meter per second squared). The magnitude of the V tendency intensifies during this period with local maxima on the order of 10^{−4} to 10^{−3} m s^{−2} (Fig. 3). Two forcing terms, PGF and COR, are a few times larger than the ADV term but generally offset each other, making the ADV term of comparable importance in determining the tendency. The CUV term for V tendency is generally small and thus not shown in Fig. 3. The residual, obtained from Eq. (1) with P_{V} equal to 0, is always smaller than ${\mathrm{10}}^{\mathrm{7}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}$ during the entire 48 h simulation (not shown). To understand how the peak error evolves with time and to avoid reaching misleading conclusions based on one or more outlying values, the evolution of the 99th percentile magnitude of the residual term is shown. Figure 5 shows that it reaches a value of about $\mathrm{7}\times {\mathrm{10}}^{\mathrm{9}}$ m s^{−2} at around 15 h. Recall that the 99th percentile magnitude of the simulated v tendency has a peak of $\mathrm{7}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2} (Fig. 2b). Thus, the relative magnitude of the 99th percentile residual is about 0.001 % of the 99th percentile tendency term during the peak intensifying stage. Compared to the V tendency, the W tendency exhibits narrower features in the horizontal direction (Fig. 4) with an overall smaller magnitude in every term. The two largest forcings, PGFBUOY and CUV, usually have opposite signs, so their combined effect is on the same order as the ADV and the W tendency term. While the contribution from the upperlayer vertical velocity damping is not shown in Fig. 4, it is included as part of the rhs (P_{W}) of Eq. (2) when calculating the residual for the inline budget analysis. The residual in the inline W budget is generally 4 orders of magnitude smaller than its tendency term. The 99th percentile residual for W budget is about $\mathrm{2}\times {\mathrm{10}}^{\mathrm{10}}$ m s^{−2}, around 0.0003 % of the 99th percentile w tendency during the peak intensifying stage of the convection (not shown).
3.2 Postprocessed momentum budget analyses
3.2.1 Key features and methodologies
In contrast to extracting terms directly from the model during its integration, most of the studies in which the momentum budget analysis is conducted use the model output files after the completion of the integration. Note that since the suboutput timestep information is not available between successive outputs, only the largestep forcing terms can be estimated in these postprocessed budget analyses. Generally, the neglect of the acoustic or smallstep modes is expected to have little impact on the results as the highfrequency modes are often considered meteorologically insignificant. However, it is mentioned in Klemp et al. (2007) and Skamarock et al. (2008) that the WRF smallstep integration scheme includes not only the acousticwave but also some gravitywave modes, which may not be insignificant. These gravitywave modes form during the smallstep integration due to the designated terms that are required for acousticwave propagation and “Consequently, in this vertical coordinate (i.e., terrainfollowing hydrostatic pressure coordinate), the terms governing the acoustic and gravity wave modes are intermingled to the extent that it does not appear feasible to evaluate any of the gravity wave terms on the large time steps, even if one desired to do so” (Klemp et al., 2007).
Most of the studies did not reveal the complete details about how their analysis was done, so we cannot presume their methodologies and the possible errors. However, a few simplifications commonly made in the postprocessed budget analyses may introduce errors that result in deviations from the simulated results and thus a significant residual. Below we revisit the relevant features of the WRF model that should be considered and discuss how they might affect the postprocessed budget if they are ignored. Then, the results are shown for different postprocessed budget analyses with different simplifications (Table 1). The aim herein is to identify these potential errors hidden in the budget calculation and show how severely they affect the resulting interpretation.
(a) Diagnosed tendency
In a postprocessed budget analysis, the tendency term of a given variable is approximated by the difference between the value of this variable at two successive output times divided by the output time interval. Thus, the accuracy may be sensitive to the output time interval. The value at the predicted state has a form of
If the output interval is longer than the model integration time step, the diagnosed tendency would deviate from the model prediction of the instantaneous tendency. To increase the accuracy, the output time interval Δt_{output} needs to be similar to the integration time step Δt.
(b) Spatial discretization on the C staggered grid
For computational efficiency and accuracy, WRF utilizes a Cgrid staggering system (Arakawa and Lamb, 1977). This staggering system is pertinent to the numerical solution for spatial derivatives. For most of the spatial derivatives other than advection (e.g., the pressure gradient force), the secondorder finite difference operator is used in the WRF model. For example, the y derivative of variable Φ is calculated using the discrete operator:
The index ($i,j,k$) corresponds to a location with $(x,y,\mathit{\eta})=(i\mathrm{\Delta}x,j\mathrm{\Delta}y,k\mathrm{\Delta}\mathit{\eta})$, where Δx, Δy and Δη are the grid lengths in the two horizontal and vertical directions (can be vertically stretched), respectively. The same expression applies for the x or the η derivatives. Grid staggering implies that different variables may be located on different grids, i.e., shifted by a halfgrid point from the others as illustrated in Fig. 6. Depending on what variable the spatial derivatives are intended for, Eq. (9) should be carried out on the corresponding grid, which is not necessarily the same as the Φ grid. For example, for the V tendency, all the associated forcing terms involving the spatial derivatives should be performed on the V grid. More specifically, to calculate the PGF term for the V tendency equation, the term $\frac{\partial p}{\partial y}$ and the term $\frac{\partial p}{\partial \mathit{\eta}}$ in Eq. (1) should be calculated on the V grid but not the pressure grid (p grid). Applying Eq. (9) for $\frac{\partial p}{\partial y}$, the V grid with location indices of $(i,j\frac{\mathrm{1}}{\mathrm{2}},k$) and $(i,j+\frac{\mathrm{1}}{\mathrm{2}},k$) falls exactly on the p grid and hence no interpolation is required (red arrows in Fig. 6a). However, for $\frac{\partial p}{\partial \mathit{\eta}}$, the pressures on the V grid with indices of $(i,j,k\frac{\mathrm{1}}{\mathrm{2}}$) and $(i,j,k+\frac{\mathrm{1}}{\mathrm{2}}$) must be obtained (red arrows in Fig. 6b) through linear interpolation using their surrounding closest four pressure values, e.g.,
which is weighted by the irregular (stretched) vertical gridlengths (Fig. 6b).
If the Cgrid staggering is not considered during the postprocessing analysis, i.e., all the variables have been interpolated on the universal grids before carrying out the budget calculation, in addition to the potential errors brought on by the interpolation method, the term $\frac{\partial p}{\partial y}$, for example, would essentially involve pressure differences over a larger grid interval of 2×Δy instead of Δy, with larger associated truncation errors.
(c) Advection operators
For advection, higherorder operators for finite differencing are provided as the default WRF setup. Taking the y component of the fluxform advection for V momentum in Eq. (3) as an example, with a fifthorder operator as selected in the present simulation, it is written as
where V and v are the masscoupled and massuncoupled velocities, respectively:
and
The oddorder advection operators include a spatially centered evenorder operator and an upwind diffusion term. A detailed discussion on the advection scheme in the WRF model with differentorder operators can be found in Wicker and Skamarock (2002) and Skamarock et al. (2008). Simplifying the advection estimation using an operator with an order that differs from the numerical setup would contribute to errors in the ADV estimation.
(d) Forward or backward Euler method
Conceptually, the WRF model can be considered more of a forward scheme, i.e., using the known variables from the current state to calculate the forcing and then advancing the variables forward until reaching the prediction time. However, there are a few implicit components during the integration. For example, as discussed in Sect. 2.1, the largestep forcings are updated using a predictor–corrector method in the second and third RK3 steps. In addition, the W equation is coupled with the geopotential tendency equation and includes a forwardintime weighting that utilizes predicted states of the geopotential and temperature in solving the W (Eq. 3.11, 3.12, and 3.19 in Skamarock et al., 2008).
In numerical analysis for solving ordinary differential equations, the (explicit) forward Euler method approximates the change of a system from t to t+Δt using the current states (t), while the (implicit) backward Euler method finds the solution using the predicted states (t+Δt):
Consistent with this concept, the rhs forcing terms of a budget equation can be estimated using two different instantaneous states in analogous ways. However, we emphasize that the postprocessed budget analysis does not solve the tendency equation per se but only diagnoses the relationship between the two sides of the equation. Note that for postprocessing analyses, the availability of the data depends on the output time interval (Δt_{output}), which is often much larger than the integration time step (Δt). Thus, for the tendency at a given time t+Δt, when applying the forward Euler method to estimate the associated rhs forcings, the “current states” one can use are the most recent prior output at $t+\mathrm{\Delta}t\mathrm{\Delta}{t}_{\mathrm{output}}$ (see Fig. 7):
If Δt_{output} is the same as Δt, Eq. (14) reverts to Eq. (12). If Δt_{output} is much larger than Δt, the backward Euler method using predicted states at t+Δt may better estimate the true model forcing terms as they are calculated using variables at a closer time to the real integration window in the model (Fig. 7).
The above two diagnostic methods estimate the forcing terms using instantaneous states. However, as mentioned in Sect. 3.2.1(a), the diagnosed lhs tendency depends on two successive model output times. Thus, an average between forcings diagnosed explicitly and implicitly are often considered. For a postprocessed analysis, this translates into estimating the forcings using both predicted states and the most recent prior available current states:
(e) Flux or advective form of equation
While the momentum equations solved in the WRF model are in flux form, their corresponding advective forms can be derived and are often used for postprocessed budget analyses for convenience. To derive the advective form, the fluxform V momentum equation (Eq. 1 excluding residual) is first multiplied by a factor of $\frac{\mathrm{1}}{{\mathit{\mu}}_{\mathrm{d}}}$ and V is rewritten as μ_{d}v:
Then, by adding the mass continuity equation in WRF (multiplied by a factor of $\frac{v}{{\mathit{\mu}}_{\mathrm{d}}}$):
to the rhs of Eq. (16), we obtain
Moving the first term on the rhs of Eq. (17) to the lhs, the second rhs term can be combined with the fluxform advection using the vector identity $\mathrm{\nabla}\cdot \left({\mathit{\mu}}_{\mathrm{d}}\mathit{v}\right)={\mathit{\mu}}_{\mathrm{d}}\left(\mathrm{\nabla}\cdot \mathit{v}\right)+\mathit{v}\cdot \left(\mathrm{\nabla}{\mathit{\mu}}_{\mathrm{d}}\right)$. Then, the advective form of the horizontal momentum equation is obtained as
3.2.2 Results of horizontal momentum budget
Table 1 summarizes all the postprocessed budget analyses tested in this study. In the present section, we first present the results one by one, and then a qualitative intercomparison among them and the inline retrieval method is discussed. The first postprocessed method (POST10minE) for V budget follows all the approaches in the model as closely as possible using the 10 min output data. The fluxform equation, C staggering grids, and the same orders of advection operators as the experimental setup are used. The diagnosis of the largestep forcing is applied directly to the model outputs on η levels using the explicit or forward Euler method as shown in Eq. (14). The diagnosed forcing terms are compared with their corresponding true values from the inline retrieval (Fig. 8). Errors smaller than, but on the same order of 10^{−4} m s^{−2} as the V tendency, are observed in all terms including the diagnosed tendency term. These errors grow in magnitude and areal coverage with the growth of the disturbance. Aside from COR, the absolute errors in the tendency, ADV and PGF can exceed $\mathrm{6}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2}, the former two of which are more than 50 % of the magnitude of their true (instantaneous) values locally.
The second postprocessed analysis (POST1minE) is done following the same approach but applied to the 1 min (same as the integration time step for this simulation) output data, and the results show strongly reduced errors in all terms (Fig. 9). The errors that remain are mostly in the PGF term and likely stem from the fact that the smallstep modes and the RK3 integration scheme are not considered in the postprocessed budget. These inherent errors result in a small residual term with a general order of 10^{−5} m s^{−2}, 1 to 2 order(s) smaller than the maximum V tendency. In terms of local maxima, the 99th percentile magnitude of the residual obtained in POST1minE gives a relative magnitude of about 7 % of the 99th percentile v tendency during the peak intensifying stage of the convection at around 15 h (Figs. 2b and 5). Although reducing the model output interval to be close to the integration time step helps to balance the budget without the need for inline diagnoses, it is computationally expensive especially for large, dataintensive simulations.
Given that computational cost is often a major consideration, we also test whether the implicit or backward Euler method (POST10minI) can improve the estimation of instantaneous forcing terms relative to the explicit method for the same 10 min output data (POST10minE). POST10minI follows the same strategy as POST10minE except that all the rhs terms, following Eq. (13), are diagnosed with the predicted states instead of the previous output states. As depicted in Fig. 10, POST10minI does indeed better capture the true model estimated forcing values as errors in all the rhs forcing terms diminish greatly to an accuracy similar to POST1minE. However, as these forcings are calculated at a given instant, the imbalance of the budget would remain if the diagnosed tendency term is not calculated instantaneously (the second column from the right in Fig. 10). Therefore, if budget analysis at an instant of time is desired, we recommend adding the tendency calculation within the model as a standard output and diagnosing the forcing terms implicitly, which yields a residual term on a similar order to the one obtained in POST1minE (the rightmost column in Fig. 10 and Fig. 5a).
For the more common situation, the postprocessed analyses diagnose rhs terms using two successive outputs over an output time interval, i.e., taking the averages of the explicitly and implicitly calculated forcings using Eq. (15) on the 10 min output (POST10min(E+I)/2). Comparing the averaged rhs forcings with the analogously diagnosed lhs momentum tendency (Eq. 8) gives a small residual to a similar accuracy level as POST1minE and POST10minI (the rightmost column in Figs. 11 and 5b).
We now investigate the impact of other common simplifications on top of the reference experiment, POST10min(E+I)/2. The first such simplification is to approximate the fluxform advection term using the secondorder operator (Eq. 9) for both vertical and horizontal components (POST2oadv(E+I)/2) instead of the third and fifthorder operators as used in the model setup. In our simulation, such inconsistency of advection operators introduced errors in the ADV term with a maximum value $>\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2}, more than 50 % of its true magnitude along the slantwise convective band (Fig. 12). Next, we repeated POST10min(E+I)/2 but the calculation is applied after all the model output variables have been interpolated to the universal or unstaggered grid (pressure grid) (POSTnonstag(E+I)/2). This is a common way to postprocess model output data for plotting purposes. As mentioned earlier, this approach would reduce the accuracy when solving the spatial differential terms, and indeed, the results do indicate significant errors over a large area in both ADV and PGF (Fig. 13). Their combined errors result in widespread residual values $>\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ m s^{−2} even over the area where the tendency term is smaller than $\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2} (error is at least of 30 % magnitude of the tendency term over a wide area and is reaching 100 % over the band head).
Finally, a different format of the V equation, the advective form, is used for postprocessed analysis (POSTadvF(E+I)/2). Mathematically, the fluxfrom momentum equation can be rewritten in the advective form without making any additional approximation, only with the aid of the conservation law of dryair mass in the WRF model as shown in Eqs. (16)–(18). However, during the interchange of the expression for the tendency and advection terms, truncation errors may be introduced. We reiterate that the tendency term in the advective form is not equivalent to the one in the flux form divided by μ_{d}; however, calculation suggests that they are approximately equal, i.e.,
with a maximum error that is on the order of ${\mathrm{10}}^{\mathrm{7}}\sim {\mathrm{10}}^{\mathrm{8}}$ m s^{−2} (3 orders of magnitude smaller than the simulated v tendency) in our study. The summation of the tendency term and advection term in these two forms of the momentum equation should be mathematically identical, so we would expect to see a small difference in the advection term as in the tendency term. However, we find that the advection term in the advective form has a strong positive bias compared to that in the flux form (Fig. 14). The residual term in the POSTadvF(E+I)/2 is thus negatively biased over the entire convective band with a magnitude exceeding $\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2} (reaching 100 % error near the upper half of the convective band). If the residual is neglected or not shown, authors and/or readers may falsely consider the advection process to be the dominant term governing the evolution of the slantwise updraft.
A quantitative comparison of the 99th percentile of the magnitude of the residual term in the domain (excluding the boundaries) among different analysis methods is shown in Fig. 5. The residuals between the instantaneously diagnosed forcings and the true model tendency term (calculated inline) are shown in Fig. 5a while the ones between the averaged forcings of two consecutive outputs and the diagnosed tendency term are shown in Fig. 5b. The evolution of the 99th percentile residual shows generally larger magnitudes when the momentum tendency is larger (Fig. 2b), suggesting that these errors may amplify in stronger convection cases. While the postprocessed budget analysis in POST1minE, POST10minI, and POST10min(E+I)/2 can achieve a relatively small 99th percentile residual (peak at $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ m s^{−2} , or about 7 % of the concurrent 99th percentile v tendency), the inline budget analysis always gives a much smaller magnitude ($<{\mathrm{10}}^{\mathrm{8}}$ m s^{−2}, or 0.001 % of the tendency, during the entire simulation). Figure 5 also shows that any simplification that is inconsistent with the model solver can severely degrade the accuracy of the postprocessed budget analysis. Both POSTnonstag(E+I)/2 and POSTadvF(E+I)/2 can lead to a 99th percentile of the residual magnitude peaking at around $\mathrm{4}\times {\mathrm{10}}^{\mathrm{4}}$ m s^{−2} or more, which corresponds to >50 % of their concurrent 99th percentile simulated v tendency, respectively. Generally, a higher relative magnitude of residual to v tendency is reached if the maximum instead of the 99th percentile is examined (despite larger fluctuation with time). We also examined the 95th percentile of the residual magnitude and obtained qualitatively similar results although the relative magnitudes of such chosen residuals among the three postprocessing methods with simplifications (POST2oadv(E+I)/2, POSTnonstag(E+I)/2, and POSTadvF(E+I)/2) vary due to their different error distributions.
3.2.3 Results of vertical momentum budget
For the W equation, the closure of the postprocessed budget appears not to be practicable even when the output time interval is reduced to the integration time step (Fig. 15). One partial reason is that the spatially noisy smallstep modes, neglected in the offline budget analysis, are surprisingly large with a general order of 10^{−4} m s^{−2} over the growing band, which is 1 order of magnitude larger than the W tendency (see the blue and red contours overlapped on the residual subplots in Fig. 4). These highfrequency modes not only include vertically propagating sound waves but also some gravity wave modes (Klemp et al., 2007). Furthermore, as indicated in Eq. (6) and mentioned in Sect. 3.1, the W equation solved in the WRF model is implicit, coupled with geopotential tendency equation and includes a forwardintime averaging operator that is applied to the smallstep modes:
where β is a userspecified parameter and Δτ indicates the small time step in the acoustic scheme. This means that the smallstep modes at a current small step, ${\stackrel{\mathrm{\u203e}}{\left(\mathrm{ACOUS}\right)}}^{\mathit{\tau}}$, are calculated using information (e.g., geopotential, potential temperature and density) at the forecast time τ+Δτ (see Eq. 3.11 and 3.12 in Skamarock et al., 2008). All these components are not feasible for an offline budget calculation.
The application of POST1minE for the W tendency shows that this method accurately estimates most of the processes, but large errors $>\mathrm{2}\times {\mathrm{10}}^{\mathrm{3}}$ m s^{−2} remain in the PGFBUOY term resulting in a widespread residual that reaches the same magnitude of the peak W tendency term (Fig. 15). The fact that these errors exceed the smallstep modes (contributing to PGFBUOY) suggests that such imbalance does not solely come from the neglect of the smallstep modes. A close comparison of the postprocessed and the inline PGFBUOY shows that our estimation is close to the inline value to an accuracy of at least three significant figures at the first RK3 step before the acoustic contribution is considered (not shown). However, this largestep forcing term adjusts rapidly, sometime even with a sign change, from step to step within the RK3 integration. Although it is feasible to estimate F(Φ^{t}) via postprocessing, it is however impossible to retrieve $F\left({\mathrm{\Phi}}^{**}\right)$ in Eq. (4), leading to the poor estimation of vertical pressure gradient and buoyancy force in the W budget. This result also suggests that the budget closure for vertical velocity is difficult by nature due to its rapid variation on small scales.
The growth of the residual as the convection intensifies (Fig. 5) motivates a test for a different case with stronger momentum tendencies. A WRF idealized 2D squall line test case (em_squall2d_y; Skamarock et al., 2008) is selected with a horizontal resolution of 250 m and 3 s integration time step, and the simulation is integrated for 1 h. A subgrid turbulence scheme based on the prognostic turbulent kinetic energy equation is activated (diff_opt=2 and km_opt=2; Skamarock et al., 2008, chap. 4.2.4). The simulated v tendency in this case is 2 orders of magnitude stronger than the one in the slantwise convection case. The inline retrieval budget tool works well with 99th percentile residuals generally 2 orders smaller than the tendency terms in the domain during this simulation. However, as compared to the slantwise convection case, this case features a larger relative magnitude of 99th percentile residual to the 99th percentile tendency term of about 0.1 %. Furthermore, the postprocessed budget analysis applied to the output data with an output interval the same as the integration time step (analogous to POST1minE but in this case, it is termed POST3secE; Fig. 16), with no simplification made, does not work as well as in the slantwise convection case. POST3secE shows that the largest error appears in the PGF term with a magnitude of 50 % of its true value at a given instant. The error in diffusion only accounts for about 10 % of the error at the same time. One possible reason is that unlike the case of slantwise convection where the PGF exhibits a horizontally rather uniform structure with almost the same sign (Fig. 3), the PGF term in this case has a more complex spatial structure with several sign changes over a horizontal distance of 10–15 km. Thus, large errors appear at the edge of these positive or negative patches where the sign changes. Despite the small spatial scales of these errors, the large error magnitude would render accurate interpretation of the physical process difficult based on such postprocessed budget analysis. This result suggests that the postprocessed budgets, even when done with care, do not always work well, and that the associated residual or errors might be sensitive to the intensity of the simulated system, the spatial or temporal resolution, and the nature of the physical processes governing the different systems.
While an increase in spatial resolution often requires a shorter integration time step for numerical stability and may result in stronger simulated convection, it is almost impossible to separate all these factors. We can, however, conduct the same slantwise convection simulation with a higher resolution of 2 km (and a shorter integration time step of 10 s) to exclude the effect of different physical processes in different systems and discuss the changes in the accuracy of the budget analysis when spatial resolution is increased from 10 km. As shown in Fig. 2b, in the 2 km simulation the maximum of the simulated 99th percentile v tendency is $\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{3}}$ m s^{−2}, almost twice the magnitude in the 10 km run. The magnitude of the residual from the inline budget analysis also becomes larger with the 99th percentile value almost 1 order larger than that in the 10 km simulation (Fig. 5). However, its relative magnitude is still small and amounts to about 0.005 % of the tendency in the 2 km case. For the postprocessed budget analysis applied to the 2 km simulation, the 99th percentile residual with the instantaneous calculation of POST10minI2km appears only slightly larger yet sometimes smaller than those in its 10 km case (Fig. 5a). For the method using two model outputs for both diagnosed tendency and forcing terms, the peak 99th percentile residual in POST10min(E+I)/22km is about 4 times larger than that in its 10 km counterpart (POST10min(E+I)/2). This is likely due to the larger deviation caused by the longer diagnosed window (10 min) with respect to the integration time step (10 s) in the 2 km case. In addition, it appears that the simulated fields adjust more rapidly with more complex structures on smaller scales in the 2 km simulation as compared to the 10 km simulation (not shown). If the same analysis is performed using the 1 min output (POST1min(E+I)/22km) as opposed to the 10 min output, the residual can be greatly reduced to be similar to that obtained in POST10min(E+I)/2 (Fig. 5b).
The results presented above suggest that the relative magnitude of errors in budget analysis vary with different systems or cases. Furthermore, while the absolute errors in the inline momentum budget analyses do indeed increase with increasing horizontal resolution, the relative magnitude with respect to the simulated tendency does not increase substantially. The accuracy of the postprocessed budget analysis using the averages of two consecutive model outputs is highly dependent on the ratio of the output interval and the integration time step. A ratio of 10 as used in the POST10min(E+I)/2 results in an acceptable accuracy (99th percentile residual of about 7 % of the tendency), while a lower value of 6 is required for highresolution simulations (e.g., the 2 km case) to reach a similar accuracy. For cases with a more complex physical process like the squall line test case, the inline budget retrieval appears necessary for adequate budget closure.
Budget analysis is a commonly employed tool in numerical studies to understand the underlying mechanisms for certain simulated features of interest. However, many studies still have difficulties in achieving a balanced or closed budget especially when a fullphysics model is used and when the budget is calculated instantaneously over a local area. Aside from the complexity of various (some implicit) parameterization schemes, the main challenge in closing the budget involves the analysis of postprocessed data using algorithms that are inconsistent with the model solver. In this study, an inline momentum budget retrieval tool is developed for the WRF model, and its advantages for momentum budget analysis are demonstrated. The 99th percentile residual obtained from this inline retrieval is always smaller than or about 0.1 % of the actual tendency term in all the tested cases, which include idealized, 2D simulations of slantwise convection and squall lines. Taking the results from the inline retrieval as “truth”, we investigate the potential errors in each term and the resultant residual for postprocessed budget analyses under different assumptions.
The comparison among different postprocessed diagnoses is focused on the horizontal momentum (V) budget. The reason is that postprocessed vertical momentum (W) budget analysis fails to produce reasonably accurate results due to the noisy vertical pressure gradient and buoyancy forces that are tied closely to the small timestep modes and the implicit scheme used for the vertical momentum integration. Thus, inline retrieval is necessary for an accurate W budget analysis. The errors in the postprocessed V budget arise from both the lefthandside tendency term and the righthandside (rhs) forcing terms. To improve the accuracy of the diagnosed momentum tendency estimation, one can reduce the output interval to the model integration time step, which incurs a large computational cost and consumes a large amount of disk space. An alternative and cheaper solution is to add the tendency calculation within the model as a standard output. Our test case of slantwise convection shows that the diagnosed tendency using two successive model outputs with a 10 min interval to approximate the instantaneous true tendency (with an integration time step of 1 min) could create an error exceeding 50 %, which greatly limits the effectiveness of such a budget for physical interpretation.
For the rhs forcing terms in the V equation, errors can be limited if the postdiagnosis is done with care using the same form of the model equation, the same spatial discretization, and the same order of the advection operators and performing the calculation on the original (e.g., C staggering and vertically stretched) model grids. However, these steps are necessary but not necessarily sufficient for the closure of the budget, as the forcing term diagnosis also largely depends on the selected input states. If the budget at an instant of time is desired, the explicit or forward Euler method using the previous states might result in large and widespread errors in the advection and horizontal pressure gradient terms (local peak errors are about 50 % and 25 % of their true values in our simulation, POST10minE) unless the output interval is reduced to the integration time step. In the latter case (POST1minE), an error <5 % for each individual term and a residual generally 1 to 2 order(s) smaller than the maximum tendency can be achieved (the 99th percentile residual is about 7 % of the 99th percentile v tendency). An alternative way to reach a similar level of accuracy for instantaneous fields without compromising the computational cost is to diagnose the rhs forcings using the implicit or backward Euler method (POST1minI). This method diagnoses the forcings using the predicted states and thus can better capture the true model forcings by using inputs at a closer time to the model integration window.
Instead of performing the calculation using model output at one given instant, a more general postprocessed budget analysis can use two successive model outputs (POST10min(E+I)/2). This method seems to work well with the 99th percentile residual being about 7 % of the 99th percentile v tendency in our 10 km slantwise convection case with 10 min output intervals. However, the accuracy of this method varies among the test cases of different systems and is sensitive to the ratio of the output interval to the integration time step. Among the tests conducted in this study, an upper limit of 10 for this ratio is suggested, and it should be even smaller for highresolution simulations of highamplitude weather systems, as rapid adjustments occur on the small scales.
Three other common assumptions in postprocessing analysis are made on top of the POST10min(E+I)/2 to examine their potential impacts on the accuracy of the horizontal momentum budget analysis. First, utilizing an advection operator with a lower order than the one used in the model setup degrades the accuracy of the advection term with up to 50 % error over the area where the advection is the strongest (POST2oadv(E+I)/2). Second, the neglect of the staggering grids would negatively impact the estimation of all the spatial differential terms, leading to a widespread residual of at least 30 % of the local tendency (POSTnonstag(E+I)/2). Last, when the advective form of the momentum equation is used for postdiagnosis rather than the flux form, although it is mathematically equivalent to the flux form solved in the model solver, a strong negatively biased residual results (POSTadvF(E+I)/2). Both POSTnonstag(E+I)/2 and POSTadvF(E+I)/2 give a peak 99th percentile residual of about >50 % of the concurrent 99th percentile of the v tendency. All the above errors do not just appear randomly; rather, they are spread over the area where the dynamics are the most active, thus undermining the physical interpretation of the dynamics of the simulated system. We thus emphasize the importance of revealing the magnitude of the residual (relative to the tendency term) in publications on budget analysis, to enable readers to gauge the validity of the results.
While the postprocessed V budget analysis can reach an acceptable accuracy in some cases, the resultant residual may vary from case to case even when the same analysis method is adopted. Our test of an idealized squall line case with strong momentum tendencies shows that the application of the postprocessed budget analysis method without any simplification using the 3 s (same as the model integration time step) output data nevertheless results in a large relative error magnitude (∼50 %) in the horizontal pressure gradient force, with very smallscale error structures.
In summary, different assumptions or simplifications made in a postprocessed budget analysis may severely impact the estimation of each forcing term and result in a large imbalance of the budget. Based on our experiments, we conclude that the inline retrieval method like that developed herein is the most reliable one for budget analysis in numerical studies. While the budget analyses shown in this study are only for V and W momentum under the 2D idealized configurations, this newly developed budget tool also retrieves budget terms for U momentum and potential temperature. It can be applied to 3D idealized and real cases as the map projection is also considered, following the original governing equations as shown in Skamarock et al.’s (2008) Eq. (2.23)–(2.25) with map factors, which are equal to 1 for an idealized setup on the Cartesian coordinate. We also stress that in some budget studies where a coordinate transformation is necessary (e.g., from Cartesian to polar), some errors are unavoidable. In such cases, it is best to perform the budget calculation using the inline retrieval method on the model grid and then transform the budget to a new coordinate (e.g., Zhang et al., 2000). Finally, in situations where the inline coding cannot be done, this study also provides general guidance to minimize the error in the budget. Thus, our results are beneficial to budget analyses in numerical studies in general and are not limited to the WRF model.
To our knowledge, there are at least three other similar inline budget retrieval works that have been done in the WRF model:

Lehner (2012) applied to v3.2.1:
Lehner (2012) provides a very detailed instruction of how an inline budget retrieval is done for the WRF model. The method or code was utilized in Lehner and Whiteman (2014) to study the mechanisms of the thermally driven crossbasin circulation. However, the code was never made publicly available. From the document, it appears that Lehner's (2012) general procedure of retrieving the rhs budget terms during the model integration is essentially the same as our approach, which considers both the large timestep and the small or acoustic timestep contributions. Furthermore, the individual contribution from different parameterization schemes that are activated in her study was also separately retrieved. While the general method appears highly similar to our code, the momentum budget retrieval in Lehner (2012) only applies to the horizontal momentum (U and V) whereas our tool includes the budget retrieval for the vertical momentum (W) as well.

Moisseeva (2014) and Moisseeva and Steyn (2014) with v3.4.1:
The code is publicly available. The developed budget retrieval is also for the horizontal momentum equations only. The method is simpler than Lehner (2012) as it does not include the acoustic or smallstep correction terms. Furthermore, while the large timestep, nonparameterized terms (e.g., pressure gradient terms, advection, Coriolis terms) are individually retrieved, their modified registry file only outputs one (summarized) term for all the parameterized physics.

Potter et al. (2018) with v3.8.1:
The code is publicly available. This budget retrieval uses the code adapted from Moisseeva (2014), taking references from Lehner (2012), and is applied to the same version of the WRF model as used in this study (v3.8.1). More components are added from the version used in Moisseeva (2014), including the potential temperature budget, vertical velocity budget, the sixthorder diffusion term, the parameterized physics term decomposed to boundary layer, and radiation schemes. A major difference from our retrieval tools exists in that the smallstep components are neglected in Potter et al. (2018). Comparing the budget analysis results using our retrieval tool with those using theirs for the same idealized test case of the 2D squall line, the largest differences appear in terms that involve the smallstep contributions (e.g., PGF and PGFBUOY), which result in larger residual terms with Potter et al.'s (2018) retrieval method (not shown). While the relative magnitudes of such residuals to the tendency term still appear small for the horizontal momentum budget, they become larger for the vertical momentum budget. This is consistent with our result that the smallstep modes are more important in the W budget equation than in the V budget equation, and thus ignoring them results in larger errors.
Furthermore, calculations of the lhs tendency terms are added as new variables in our tool while the tendency terms used in the above studies are the model variables ru_tend, rv_tend, rw_tend, etc., which only represent the summation of all the largestep forcings to their corresponding fields (can be directly outputted via changing the WRF registry file) instead of their true local changes with time.
To construct an initial condition that contains conditional symmetric instability (CSI) but to avoid dry symmetric instability and dry and conditional (gravitational) instability is a challenging task (Persson and Warner, 1995). Therefore, the initial profile in our test case is decided by a trialanderror method and follows the following steps.

We first prescribe a horizontally uniform Brunt–Väisäilä frequency, ${N}^{\mathrm{2}}=\frac{g}{{\mathit{\theta}}_{\mathrm{v}}}\frac{\partial {\mathit{\theta}}_{\mathrm{v}}}{\partial z}$ with a vertical profile of
$$\begin{array}{}\text{(B1)}& \left(\right)open="\{">{N}^{\mathrm{2}}=\begin{array}{c}\mathrm{1.25}\times {\mathrm{10}}^{\mathrm{4}}\\ \mathrm{9}\times {\mathrm{10}}^{\mathrm{5}}\\ \mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}\end{array}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}\phantom{\rule{0.25em}{0ex}}\mathrm{for}\phantom{\rule{0.33em}{0ex}}\begin{array}{c}z\mathrm{0.5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\\ \mathrm{5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\le z\mathrm{10.5}\phantom{\rule{0.33em}{0ex}}\mathrm{km},\\ z\ge \mathrm{13.5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\end{array}\end{array}$$where z is the height and there is a linear transition for the layers $\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\le z<\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{km}$ and $\mathrm{10.5}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\le z<\mathrm{13.5}\phantom{\rule{0.125em}{0ex}}\mathrm{km}$ using the specified values beneath and above the layers.

A constant geostrophic vertical zonal wind shear is given: $\frac{\partial {U}_{\mathrm{g}}}{\partial \mathrm{z}}=\mathrm{5.8}\times {\mathrm{10}}^{\mathrm{3}}$ s^{−1}. Thermal wind balance gives
$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\partial {U}_{\mathrm{g}}}{\partial \mathrm{z}}}={\displaystyle \frac{g}{f{\mathit{\theta}}_{\mathrm{v}}}}{\displaystyle \frac{\partial {\mathit{\theta}}_{\mathrm{v}}}{\partial y}}.\end{array}$$ 
Based on Eqs. (B1) and (B2), we can specify the value of θ_{v} at any point and then derive the θ_{v} for the entire domain. In this case, $\phantom{\rule{0.33em}{0ex}}{\mathit{\theta}}_{\mathrm{v}}\left({y}_{\mathrm{0}},{z}_{\mathrm{0}}\right)=\mathrm{287.5}$ K, where (y_{0},z_{0}) indicates the grid at the surface on the southern boundary.

The relative humidity (RH) field is constructed by specifying a horizontally uniform background profile (RH_{background}) with some enhancement (RH_{bubble}) over an elliptical area where the initial perturbation will be later added. The enhanced humidity over a limited area hastens the release of CSI and avoids convection developing near the southern boundary.
$${\mathrm{RH}}_{\mathrm{background}}\left(z\right)=\left\{\begin{array}{l}\mathrm{0.81}\\ min\left[\mathrm{0.81},\mathrm{1}\mathrm{0.9}{\left(\frac{z\mathrm{5}}{\mathrm{7.5}}\right)}^{\mathrm{0.8}}\right]\\ \mathrm{0.1}\end{array}\right)$$$$\mathrm{for}\phantom{\rule{0.33em}{0ex}}\begin{array}{c}z\le \mathrm{5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\\ \mathrm{5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}<z<\mathrm{12.5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\\ z\ge \mathrm{12.5}\phantom{\rule{0.33em}{0ex}}\mathrm{km}\end{array},$$$${\mathrm{RH}}_{\mathrm{bubble}}\left(y,z\right)={\mathrm{RH}}_{\mathrm{background}}\left(z\right)\times {f}_{\mathrm{enhancement}}\left(y,z\right),$$where
$$\begin{array}{rl}& {f}_{\mathrm{enhancement}}\left(y,z\right)=\\ & \left\{\begin{array}{l}\mathrm{1.22}\\ \mathrm{1.22}\mathrm{0.11}\left(e\mathrm{1}\right)\\ \mathrm{1}\end{array}\right)\phantom{\rule{0.33em}{0ex}}\mathrm{for}\phantom{\rule{0.33em}{0ex}}\begin{array}{c}e\le \mathrm{1}\\ \mathrm{1}e\le \mathrm{3}\\ e\mathrm{3}\end{array},\end{array}$$where $e={\left(\frac{y\mathrm{410}}{{e}_{\mathrm{b}}}\right)}^{\mathrm{2}}+{\left(\frac{z\mathrm{1}}{{e}_{\mathrm{a}}}\right)}^{\mathrm{2}}$, ${e}_{\mathrm{b}}=\mathrm{100},{e}_{\mathrm{a}}=\mathrm{3}$, and y and z are the horizontal distance from the southern boundary and height, respectively, with units of kilometers. The constructed initial profile has a maximum RH of 98.82 % over an elliptical area centered at y=410 km and z=1 km.

A constant surface pressure is specified: p_{sfc}=1000 hPa.

We then iteratively solve for the hydrostatically balanced pressure, water vapor mixing ratio, potential temperature, dry and full (moist) air density, and geostrophic zonal wind for the entire domain.
The constructed initial environment contains some CSI, which is identified by the presence of negative saturated geostrophic potential vorticity (Chen et al., 2018). In this test case, CSI only exists over the southern half of the domain and never extends higher than 5 km.
To initiate convection, a 2D bubble of potential temperature and zonal wind perturbations is inserted in the area where RH is maximized and where the saturated geostrophic potential vorticity has a value of about $\mathrm{0.2}\sim \mathrm{0.1}$ pvu. The center of the bubble, located at y_{c}=400 km and z_{c}=1.5 km, has a maximum potential temperature perturbation Δθ_{max}=0.5 K and zonal wind perturbation $\mathrm{\Delta}{u}_{\mathrm{max}}=\mathrm{6}$ m s^{−1}. Both perturbation fields decrease to 0 with radius, following Δθ=Δθ_{max}cos^{2}(0.5πr) and Δu=Δu_{max}cos^{2}(0.5πr), for $r<=\mathrm{1}$, where $r=\sqrt{(\frac{y{y}_{\mathrm{c}}}{R}{)}^{\mathrm{2}}+(\frac{z{z}_{\mathrm{c}}}{H}{)}^{\mathrm{2}}}$, R=50 km is the horizontal radius, and H=1.5 km is the vertical radius.
The standard version of WRF v3.8.1 is publicly available at http://www2.mmm.ucar.edu/wrf/users/download/get_sources.html (last access: 6 June 2019; National Center for Atmospheric Research, 2016). The inline budget retrieval tool in the WRF v3.8.1 described in this study can be found at https://doi.org/10.5281/zenodo.3373872 (Chen, 2019). In this repository, all the files that remain unchanged from the defaults are tagged as “Initial commit”. The modified files for the budget retrieval include the Registry.EM_COMMON within the directory registry; module_diag_misc.F, module_diagnostic_driver.F, and module_physics_addtendc.F within the directory phys; module_after_all_rk_steps.F, module_big_step_utilities_em.F, module_em.F, module_first_rk_step_par2.F, module_small_step_em.F, and solve_em.F within the directory dyn_em. The current version includes retrieval for terms of local tendency, advection, horizontal pressure gradient force, net force resulting from vertical pressure gradient and buoyancy, Coriolis force, curvature, upper damping (damp_opt=2 and 3), turbulence or diffusion (diff_opt=2), verticalvelocity damping (w_damping=1) and parameterized physics from the planetary boundary layer scheme (bl_pbl_physics), the radiation scheme (ra_lw_physics and ra_sw_physics), the cumulus scheme (cu_physics), and the shallow cumulus scheme (shcu_physics).
TCC designed and performed the numerical experiments under the supervision of MKY and DJK. MKY proposed the idea of comparing the inline and postprocessed budget analyses. TCC developed the code of the inline budget retrieval tool in the WRF v3.8.1 model and the postprocessed analyses. DJK provided useful suggestions to improve the work. TCC prepared the paper and all coauthors contributed to the writing and editing of the paper.
The authors declare that they have no conflict of interest.
We thank Patrick Hawbecker, Ian Dragaud, and one anonymous reviewer for their valuable comments and suggestion that helped to improve this study.
This research has been supported by the NSERC/HydroQuebec Industrial Research Chair program and the Fonds de recherche du Québec – Nature et technologies (FRQNT) doctoral research scholarship grant.
This paper was edited by Chiel van Heerwaarden and reviewed by Patrick Hawbecker and one anonymous referee.
Abarca, S. F. and Montgomery, M. T.: Essential Dynamics of Secondary Eyewall Formation, J. Atmos. Sci., 70, 3216–3230, https://doi.org/10.1175/JASD120318.1, 2013.
Andersen, J. A. and Kuang, Z.: Moist Static Energy Budget of MJOlike Disturbances in the Atmosphere of a Zonally Symmetric Aquaplanet, J. Climate, 25, 2782–2804, https://doi.org/10.1175/JCLID1100168.1, 2012.
Arakawa, A. and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, Method. Comput. Phys., Adv. Res. Appl., 17, 173–265, https://doi.org/10.1016/B9780124608177.500094, 1977.
Arnault, J., Knoche, R., Wei, J., and Kunstmann, H.: Evaporation tagging and atmospheric water budget analysis with WRF: A regional precipitation recycling study for West Africa, Water Resour. Res., 52, 1544–1567, https://doi.org/10.1002/2015WR017704, 2016.
Balasubramanian, G. and Yau, M. K.: Baroclinic Instability in a TwoLayer Model with Parameterized Slantwise Convection, J. Atmos. Sci., 51, 971–990, https://doi.org/10.1175/15200469(1994)051<0971:BIIATL>2.0.CO;2, 1994
Bryan, G. H. and Fritsch, J. M.: A benchmark simulation for moist nonhydrostatic numerical models, Mon. Weather Rev., 130, 2917–2928, https://doi.org/10.1175/15200493(2002)130<2917:ABSFMN>2.0.CO;2, 2002.
Chen, T.C.: WRFV3.8.1_inline_budget_retrieval: Inline budget retrieval tool for 3D momentum components and potential temperature, Zenodo, https://doi.org/10.5281/zenodo.3373872, 2019.
Chen, T.C., Yau, M. K., and Kirshbaum, D. J.: Assessment of Conditional Symmetric Instability from Global Reanalysis Data, J. Atmos. Sci., 75, 2425–2443, https://doi.org/10.1175/JASD170221.1, 2018.
Duran, P. and Molinari, J.: Tropopause Evolution in a Rapidly Intensifying Tropical Cyclone: A Static Stability Budget Analysis in an Idealized Axisymmetric Framework, J. Atmos. Sci., 76, 209–229, https://doi.org/10.1175/JASD180097.1, 2019.
Gallus, W. A. and Johnson, R. H.: The momentum budget of an intense midlatitude squall line, J. Atmos. Sci., 49, 422–450, https://doi.org/10.1175/15200469(1992)049<0422:TMBOAI>2.0.CO;2, 1992.
Grell, G. A., Dudhia, J., and Stauffer, D.: A description of the fifthgeneration Penn State/NCAR Mesoscale Model (MM5), NCAR Tech. Note NCAR/TN398+STR, 121 pp., https://doi.org/10.5065/D60Z716B, 1994.
Hodur, R. M. and Fein, J. S.: A vorticity budget over the Marshall Islands during the spring and summer months, Mon. Weather Rev., 105, 1521–1526, https://doi.org/10.1175/15200493(1977)105<1521:AVBOTM>2.0.CO;2, 1977.
Huang, Y.H., Wu, C.C., and Montgomery, M. T.: Concentric eyewall formation in Typhoon Sinlaku (2008), Part III: horizontal momentum budget analyses, J. Atmos. Sci., 75, 3541–3563, https://doi.org/10.1175/JASD180037.1, 2018
Kanamitsu, M. and Saha, S.: Systematic tendency error in budget calculations, Mon. Weather Rev., 124, 1145–1160, https://doi.org/10.1175/15200493(1996)124<1145:STEIBC>2.0.CO;2, 1996.
Kiranmayi, L. and Maloney, E. D.: Intraseasonal moist static energy budget in reanalysis data, J. Geophys. Res., 116, D21117, https://doi.org/10.1029/2011JD016031, 2011.
Kirshbaum, D. J., Merlis, T. M., Gyakum, J. R., and McTaggartCowan, R.: Sensitivity of idealized moist baroclinic waves to environmental temperature and moisture content, J. Atmos. Sci., 75, 337–360, https://doi.org/10.1175/JASD170188.1, 2018.
Klemp, J. B., Skamarock, W. C., and Dudhia, J.: Conservative splitexplicit time integration methods for the compressible nonhydrostatic equations, Mon. Weather Rev., 135, 2897–2913, https://doi.org/10.1175/MWR3440.1, 2007.
Kornegay, F. C. and Vincent, D. G.: Kinetic Energy Budget Analysis During Interaction of Tropical Storm Candy (1968) with an Extratropical Frontal System, Mon. Weather Rev., 104, 849–859, https://doi.org/10.1175/15200493(1976)104<0849:KEBADI>2.0.CO;2, 1976.
Kuo, Y. and Anthes, R. A.: Accuracy of diagnostic heat and moisture budgets using SESAME79 field data as revealed by observing system simulation experiments, Mon. Weather Rev., 112, 1465–1481, https://doi.org/10.1175/15200493(1984)112<1465:AODHAM>2.0.CO;2, 1984.
Lee, C.S.: The bulk effects of cumulus momentum transport in tropical cyclones, J. Atmos. Sci., 41, 590–603, https://doi.org/10.1175/15200469(1984)041<0590:TBEOCM>2.0.CO;2, 1984.
Lehner, M.: Observations and largeeddy simulations of the thermally driven crossbasin circulation in a small, closed basin, Ph.D. thesis, University of Utah, available at: https://collections.lib.utah.edu/ark:/87278/s61n8fxw (last access: 13 December 2019), 2012.
Lehner, M. and Whiteman, D. C.: Physical mechanisms of the thermally driven cross‐basin circulation, Q. J. Roy. Meteor. Soc., 140, 895–907, https://doi.org/10.1002/qj.2195, 2014.
Lilly, D. K. and Jewett, B. F.: Momentum and kinetic energy budgets of simulated supercell thunderstorms, J. Atmos. Sci., 47, 707–726, https://doi.org/10.1175/15200469(1990)047<0707:MAKEBO>2.0.CO;2, 1990.
Liu, L., Lin, Y.L., and Chen, S.H.: Effects of landfall location and approach angle of an idealized tropical cyclone over a long mountain range, Front. Earth Sci., 4, 1–14, https://doi.org/10.3389/feart.2016.00014, 2016.
Markowski, P. M. and Richardson, Y.: Mesoscale Meteorology in Midlatitudes, WileyBlackwell, 424 pp., https://doi.org/10.1002/9780470682104, 2010.
Moisseeva, N.: Dynamical analysis of sea breeze hodograph rotation in Sardinia, B.Sc. thesis, University of British Columbia, available at: http://hdl.handle.net/2429/46069 (last access: 18 December 2019), 2014.
Moisseeva, N. and Steyn, D. G.: Dynamical analysis of seabreeze hodograph rotation in Sardinia, Atmos. Chem. Phys., 14, 13471–13481, https://doi.org/10.5194/acp14134712014, 2014.
National Center for Atmospheric Research: WRFARW Code Version 3.8.1, available at: https://www2.mmm.ucar.edu/wrf/users/download/get_sources.html (last access: 6 June 2019), 2016.
Persson, P. O. G. and Warner, T. T.: The Nonlinear Evolution of Idealized, Unforced, Conditional Symmetric Instability: A Numerical Study, J. Atmos. Sci., 52, 3449–3474, https://doi.org/10.1175/15200469(1995)052<3449:TNEOIU>2.0.CO;2, 1995.
Pielke, R. A., Cotton, W. R., Walko, R. L., Tremback, C. J., Lyons, W. A., Grasso, L. D., Nicholls, M. E., Moran, M. D., Wesley, D. A., Lee, T. J., and Copeland, J. H.: A comprehensive meteorological modeling system – RAMS, Meteor. Atmos. Phys., 49, 69–91, https://doi.org/10.1007/BF01025401, 1992.
Potter, E. R., Orr, A., Willis, I. C., Bannister, D., and Salerno, F.: Dynamical drivers of the local wind regime in a Himalayan valley, J. Geophys Res.Atmos., 123, 13186–13202, https://doi.org/10.1029/2018JD029427, 2018.
Powers, J. G., Klemp, J. B., Skamarock, W. C., Davis, C. A., Dudhia, J., Gill, D. O., Coen, J. L., Gochis, D. J., Ahmadov, R., Peckham, S. E., Grell, G. A., Michalakes, J., Trahan, S., Benjamin, S. G., Alexander, C. R., Dimego, G. J., Wang, W., Schwartz, C. S., Romine, G. S., Liu, Z., Snyder, C., Chen, F., Barlage, M. J., Yu, W., and Duda, M. G.: The weather research and forecasting model: Overview, system efforts, and future directions, B. Am. Meteorol. Soc., 98, 1717–1737, https://doi.org/10.1175/BAMSD1500308.1, 2017.
RiosBerrios, R., Torn, R. D., and Davis, C.: An ensemble approach to investigate tropical cyclone intensification in sheared environments. Part I: Katia (2011), J. Atmos. Sci., 73, 71–93, https://doi.org/10.1175/JASD150052.1, 2016.
Sanders, F. and Emanuel, K. A.: The momentum budget and temporal evolution of a mesoscale convective system, J. Atmos. Sci., 34, 322–330, https://doi.org/10.1175/15200469(1977)034<0322:tmbate>2.0.co;2, 1977.
Satoh, M.: Conservative Scheme for the Compressible Nonhydrostatic Models with the Horizontally Explicit and Vertically Implicit Time Integration Scheme, Mon. Weather Rev., 130, 1227–1245, https://doi.org/10.1175/15200493(2002)130<1227:CSFTCN>2.0.CO;2, 2002.
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Duda, M. G, Huang, X.Y., Wang, W., and Powers, J. G.: A description of the advanced research WRF version 3, NCAR Tech. Note NCAR/TN475+STR, 113 pp., https://doi.org/10.5065/D68S4MVH, 2008.
Thompson, G., Field, P. R., Rasmussen, R. M., and Hall, W. D.: Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part II: implementation of a new snow parameterization, Mon. Weather Rev., 136, 5095–5115, https://doi.org/10.1175/2008MWR2387.1, 2008.
Trier, S. B., LeMone, M. A., and Skamarock, W. C.: Effect of threedimensional structure on the stormwide horizontal accelerations and momentum budget of a simulated squall line, Mon. Weather Rev., 126, 2580–2598, https://doi.org/10.1175/15200493(1998)126<2580:EOTDSO>2.0.CO;2, 1998.
Undén, P., Rontu, L., Järvinen, H., Lynch, P., Calvo, J., Cats, G., Cuxart, J., Eerola, K., Fortelius, C., GarciaMoya, J. A., Jones, C., Lenderlink, G., McDonald, A., McGrath, R., Navascues, B., Nielsen, N. W., Ødegaard, V., Rodriguez, E., Rummukainen, M., Rõõm, R., Sattler, K., Sass, B. H., Savijärvi, H., Schreur, B. W., Sigg, R., The, H., and Tijm, A.: HIRLAM5 scientific documentation, Swed. Meteorol. and Hydrol. Inst., 144 pp., available at: http://www.hirlam.org/ (last access: 9 August 2019), 2002.
Wicker, L. J. and Skamarock, W. C.: Timesplitting methods for elastic models using forward time schemes, Mon. Weather Rev., 130, 2088–2097, https://doi.org/10.1175/15200493(2002)130<2088:TSMFEM>2.0.CO;2, 2002.
Xue, M., Droegemeier, K. K., and Wong, V.: The Advanced Regional Prediction System (ARPS) – A multiscale nonhydrostatic atmospheric simulation and prediction tool. Part I: Model dynamics and verification, Meteor. Atmos. Phys., 75, 161–193, https://doi.org/10.1007/s007030070003, 2000.
Xue, M., Droegemeier, K. K., and Wong, V., Shapiro, A., Brewster, K. A., Carr, F. H., Weber, D. B., Liu, Y., and Wang, D.: The Advanced Regional Prediction System (ARPS) – A multiscale nonhydrostatic atmospheric simulation and prediction tool. Part II: Model physics and applications, Meteor. Atmos. Phys., 76, 143–165, https://doi.org/10.1007/s007030170027, 2001.
Zhang, D.L., Liu, Y., and Yau, M. K.: A multiscale numerical study of Hurricane Andrew (1992). Part III: dynamically induced vertical motion, Mon. Weather Rev., 128, 3772–3788, https://doi.org/10.1175/15200493(2001)129<3772:AMNSOH>2.0.CO;2, 2000.
 Abstract
 Introduction
 Model and numerical setup
 Methodology and results
 Tests on different cases or with different horizontal resolutions
 Discussion and summary
 Appendix A
 Appendix B
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model and numerical setup
 Methodology and results
 Tests on different cases or with different horizontal resolutions
 Discussion and summary
 Appendix A
 Appendix B
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References