Numerically consistent budgets of energy, momentum and mass in Cartesian coordinates: Application to the WRF model

Numerically accurate budgeting of the forcing terms in the governing equations of a numerical weather prediction model is hard to achieve. Because individual budget terms are generally two to three orders of magnitude larger than the resulting tendency, exact closure of the budget can only be achieved if the contributing terms are calculated consistently with the model numerics. We present WRFlux, an open-source software that allows precise budget evaluation for the WRF model, as well as transforma5 tion of the budget equations from the terrain-following grid of the model to the Cartesian coordinate system. The theoretical framework of the numerically consistent coordinate transformation is also applicable to other models. We demonstrate the performance and a possible application of WRFlux with an idealized simulation of convective boundary layer growth over a mountain range. We illustrate the effect of inconsistent approximations by comparing the results of WRFlux with budget calculations using a lower-order advection operator and two alternative formulations of the coordinate transformation. With 10 WRFlux, the sum of all forcing terms for potential temperature, water vapor mixing ratio and momentum agrees with the respective model tendencies to high precision. In contrast, the approximations lead to large residuals: The root mean square error between the sum of the diagnosed forcing terms and the actual tendency is one to three orders of magnitude larger than with WRFlux. Furthermore, WRFlux decomposes the resolved advection into mean advective and resolved turbulence components, which is useful in the analysis of large-eddy simulation output. 15


Introduction
Budget analysis for variables of a numerical weather prediction model is a widely used tool when examining physical processes in the atmospheric sciences. Energy and mass budgeting has been used, for instance, to understand the governing dynamics of thermally-driven circulations in the mountain boundary layer (Rampanelli et al., 2004;Lehner and Whiteman, 2014;Potter et al., 2018). Other examples, such as Lilly and Jewett (1990), Kiranmayi and Maloney (2011) and Huang et al. (2018) are 20 listed in Chen et al. (2020). In a budget analysis, the relative weight and the spatial or temporal patterns of individual forcing terms are assessed. For potential temperature, for instance, forcing terms include resolved advection, subgrid-scale diffusion and diabatic processes, such as radiative heating and latent heat release.
Large competing forcing terms adding up to a relatively small total tendency make the budget calculation particularly error-prone (Chen et al., 2020): If approximations inconsistent with the model numerics are made, the sum of all forcing terms can 25 result in an unclosed budget with a large residual, even if the relative error of each forcing term is small. To obtain reliable budget calculations, Chen et al. (2020) developed a momentum and potential temperature budget analysis tool for the widely used WRF model (Skamarock and Klemp, 2008). Their study focuses on the horizontal momentum budget in idealized 2D simulations of slantwise convection and squall lines, and demonstrates that very small residual values can be achieved by retrieving all relevant forcing terms during the runtime of the model. Chen et al. (2020) compare their results with other ways of 30 approximating the momentum budget: neglection of grid staggering, use of a lower-order advection operator and use of the advective form instead of the flux-form of the equations. The authors demonstrate that these approximations strongly deteriorate the budget closure.
Even if the equations in the numerical model are cast in flux-form, some authors use the advective form in budget analyses as they find it easier to interpret (e.g., Umek et al., 2021). Although it is possible to discretize an advective-form equation in a way 35 that it is numerically equivalent to the respective flux form (Xue and Lin, 2001), the advective form may hinder interpretation when internal processes dominate the budget of the control volume, which can be a single grid cell or a larger volume (Lee et al., 2004).
By design, the budget computation method of Chen et al. (2020) cannot discriminate between tendencies caused by resolvedscale and subgrid-scale turbulence. This is important, especially when doing large-eddy simulations, which partially resolve 40 the turbulence spectrum. A budget analysis tool capable of estimating tendencies from resolved turbulence must rely on online computation of turbulence statistics during model integration.
Furthermore, budget analysis is more intuitively carried out in the Cartesian coordinate system, while numerical weather prediction models generally adopt a curvilinear terrain-following system. For budget diagnostics in simulation domains with non-uniform orography, accurate computation of the coordinate transformation between the terrain-following and the Cartesian 45 system is therefore mandatory. This is mainly an issue when tendencies resulting from flux derivatives in a particular spatial direction, such as the vertical derivative of the resolved turbulent flux, are inspected. Some numerical weather prediction models, e.g. WRF, adopt a mass-based vertical coordinate. Because the atmospheric mass in a model column generally varies during integration, the height of the vertical levels changes with time. Thus, time derivatives on constant model levels and at constant height are not equal. This also needs to be accounted for if one wishes to compute 50 the total model tendency in the Cartesian coordinate system accurately. When looking at the instantaneous tendencies between individual model time steps, this effect can usually be neglected. However, if the budget is averaged over a time interval, the distance between the vertical levels can change considerably.
The decomposition into mean and turbulent components and the coordinate transformation to the Cartesian coordinate system were implemented, e.g., by Schmidli (2013) and Umek et al. (2021). However, neither of the two studies aim at a closed budget.

55
In this study, we present WRFlux, an open-source budget calculation tool for WRF that yields a closed budget, a consistent transformation to the Cartesian coordinate system and decomposition into mean and turbulent components. WRFlux allows to output time-averaged resolved and subgrid-scale fluxes and other tendency components for potential temperature, water vapor mixing ratio and momentum for the Advanced Research WRF (ARW) dynamical core.
The paper is organized as follows. First, we summarize the theoretical foundation of the approach in Sect. 2. This is relevant not only for the WRF model but for any hydrodynamic model in flux-form that utilizes a generalized vertical coordinate. In Sect. 3, details about the implementation of WRFlux are given, followed by the results of an example simulation in Sect. 4. The purpose of the example simulation is to illustrate a possible application of WRFlux, show its performance and compare it 65 to other, more simplified budget computation approaches.

Conservation equation transformations
The flux-form conservation equation for a variable ψ in the Cartesian coordinate system x = (x 0 , x 1 , x 2 , x 3 ) = (t, x, y, z) reads: where ρ is the air density, u i are the components of the wind speed vector and ψ is a prognostic variable (e.g., u i , potential temperature θ or mixing ratio, e.g. of water vapor). The first term on the right-hand side is the advective tendency, the second term contains all other forcing terms for ψ.
The transformation of Eq. 1 yields: where |J| is the determinant of the 4x4 Jacobi matrix of the transformation, J ij = ∂ ξj x i and ν i = J −1 ij u j = u j ∂ xj ξ i is the 80 contravariant velocity in the new coordinate system with u = (1, u, v, w). Following Liseikin (2010), we use a four-dimensional coordinate system since the coordinate transformation can be time-dependent.
Many atmospheric models use a coordinate system of the form ξ = (t, x, y, η) with the generalized vertical coordinate η . η can be a function of space and time and must possess a monotonic relationship to height z (Kasahara, 1974). The Jacobian matrix for the transformation to such a coordinate system and its inverse are given in appendix A. The Jacobian determinant reads: and the contravariant velocity vector is: Examples of generalized vertical coordinates include terrain-following coordinates and pressure-based coordinates. WRF, for instance, uses a hybrid terrain-following vertical coordinate based on hydrostatic pressure (Klemp, 2011). In WRF, η is a function of space (x, y and z) and time. The coordinate metric |J| appears as part of the dry air mass µ d = −ρ d gz η in the model equations, where ρ d is the dry-air density and g the acceleration due to gravity .

95
Inserting Eq. 3 and 4 into Eq. 2 yields: This form of the conservation equations is typically used in numerical weather prediction models. The horizontal and temporal derivatives are taken on constant η-levels. To make this clear, we continue using τ , ξ 1 and ξ 2 in the equations, even though (τ, ξ 1 , ξ 2 ) = (t, x, y).

100
For a budget calculation tool, taking the derivatives on constant η-levels is convenient since it avoids interpolation of the model output to constant height levels. However, we would like to have the individual tendency terms as in the Cartesian coordinate system (Eq. 1) for improved interpretability and for comparison with measurements. To attain both of these requirements we can transform the derivatives in Eq. 1 to be on constant η-levels. Derivatives with respect to x i and ξ i in a coordinate system with a generalized vertical coordinate are related (Kasahara, 1974;Byun, 1999) as: with i = 0, 1, 2, j = 0, 1, 2, 3 and z xi := ∂ ξi z.

110
The second term on the left-hand side and the second term in square brackets are the correction terms that account for the derivatives being natively computed on constant η instead of on constant z-levels.
One can show that equation 8 is numerically not consistent with Eq. 6. Numerically consistent means that the budget closes not only analytically but also after discretization. Therefore, we search for an alternative formulation which is equivalent to Eq. 8 and numerically consistent with Eq. 6. Starting with Eq. 6, we first replace ω with w using: This equation is analogous to the geopotential tendency equation in WRF.
Solving for z η ω, inserting in Eq. 6 and rearranging leads to: 120 Dividing by z η finally yields: Using the product rule and the commutativity of partial derivatives, one can show that equation 11 is equivalent to Eq. 8.
However, the correction terms in Eq. 11 (second term on the left-hand side and second term in square brackets) are conceptually 125 different from the correction terms in Eq. 8. While the latter only correct for the derivatives being taken on constant η-levels, the former also correct for z η being used in the temporal and horizontal derivatives.
In contrast to Eq. 8, Eq. 11 can be closed numerically. For this, we need to close the model equation in the terrain-following coordinate system (Eq. 6) and the geopotential equation (Eq. 9) and z η ω needs to be numerically equivalent in both equations.
The latter two requirements can be achieved by recalculating w based on Eq. 9 instead of using the prognostic value of the 130 model.

The θ-budget
For numerical reasons, WRF uses potential temperature perturbation as a prognostic variable. The perturbation is computed with respect to a constant base state, as θ p = θ − θ 0 with θ 0 = 300 K. Based on this decomposition, the energy equation can be split up into advection of the perturbation and of the constant base state: with the contravariant velocity ν. Due to the continuity equation, the last term on the right-hand side of Eq. 12 is identically zero. Equation 12 can be used to compute the components of the full-θ-tendency with high numerical accuracy.
WRF uses C grid staggering (Arakawa and Lamb, 1977) to discretize the governing equations due to its favorable conservation properties. For the thermodynamic variables (potential temperature and mixing ratio), we discretize equation 11 as: On the right-hand side, the operator δ denotes central finite differences of the staggered fluxes while overbars denote spatial averaging to the correct location. The averaging operation for ψ depends on the type and order of the advection operator. While the even-order advection operators are spatially centered, the odd-order operators consist of an even-order operator and an upwind term (Skamarock et al., 2019;Chen et al., 2020). In addition to the standard advection operators, WRF offers positivedefinite, monotonic and Weighted Essentially Non-Oscillatory options.

150
For Eq. 13 to be numerically consistent with the conservation equation used in the model (Eq. 6), all terms need to use the same advection operator as in the numerical model. The correction terms derive from the vertical advection term and thus must be discretized in the same way as the vertical advection.
Since the momentum components are staggered themselves the discretization differs from Eq. 13, but still straightforward to derive. We do not state them here for brevity.

155
In Eq. 8, we introduced a form of the conservation equation that follows immediately from the Cartesian conservation equation but is numerically not consistent with the terrain-following formulation (Eq. 6) used by the model. To demonstrate the inconsistency, we compare Eq.13 with two different discretizations of Eq. 8.
In the first one the corrections for the horizontal derivatives are built by taking the horizontal flux and staggering it horizontally 160 and vertically to the grid of the vertical flux: The second one adopts a different discretization in the horizontal correction term to avoid double averaging: The impact of the approximate budget calculations (Eq. 14 and 15) is discussed in Sect. 4.3.

Flux averaging and decomposition
The decomposition of the resolved fluxes into mean advective and resolved turbulent components requires averaging the fluxes over time and/or space (e.g., Schmidli, 2013). The averaging is considered as an approximation of an ensemble average.
Means and perturbations are defined by: 175 ψ denotes the time and/or spatial block average, ψ is the density-weighted average and ψ the perturbation thereof.
The decomposition of the resolved flux then reads: with the total flux on the left-hand side and the mean advective and resolved turbulent fluxes on the right-hand side.
We use the density-weighted average, also known as Hesselberg or Favre averaging (Hesselberg, 1926;Favre, 1969), because 180 WRF is a compressible model. Other studies using density-weighted averaging include Kramm et al. (1995), Greatbatch (2001) and Kowalski (2012). The budget closure is insensitive to whether or not density-weighted averaging is applied. In fact, the latter only affects the partitioning between the mean advective and resolved turbulent fluxes, but not the total flux itself. For typical atmospheric applications, also the impact on the mean advective and resolved turbulent components is hardly noticeable.
The correction flux used in the horizontal corrections in Eq. 11 can be decomposed as:

Implementation
We -The vertical velocity is recalculated with Eq. 9. The last term in Eq. 9 is formulated to be consistent with the vertical advection of the budget variable in the terrain-following coordinate system. To achieve a close match of this recalculated velocity and the prognostic one, we removed an unnecessary double-averaging of ω in the vertical advection of geopo-200 tential 1 . Except for this modification, which leads to about 10 % stronger updrafts and downdrafts, WRFlux does not change the dynamics of the WRF model.
-To close the budget for both the perturbation θ p = θ − θ 0 and full θ, the last term in Eq. 12 needs to vanish. Since WRF does not actually solve the continuity equation but instead integrates it vertically, this is not trivial. Therefore, we calculate the temporal term and horizontal divergence terms of the continuity equation explicitly and take the vertical 205 term as the residual. Using the residual has only a marginal effect on the vertical component.
-Dry θ-tendencies can be output even when the model is configured to use moist θ (Xiao et al., 2015) as the prognostic variable.
-Map-scale factors are taken care of as described in Skamarock et al. (2019). Thus, WRFlux is also suited for real-case simulations.
225 with the ridge height h m = 1500 m and x 1 = 0.5 km, x 2 = 9.5 km. The ridge-to-ridge distance is thus 20 km, equal to the domain width. The distance in x-direction is defined to be 0 at the center of the domain.
Implicit Rayleigh damping  is used above 4 km. We verified that the damping layer is sufficiently deep to dissipate vertically propagating gravity waves before they could reach the model top, be reflected and cause numerical instabilities. The advection scheme is 5th-order in the horizontal and 3rd-order in the vertical. Subgrid-scale diffusion follows 230 Deardorff (1980), with different eddy diffusivities for the horizontal and vertical to account for the anisotropic grid.
The boundary layer evolution is driven by a simplified radiation scheme as in Schmidli (2013). The radiative balance at the surface is given by: where R n is the net radiation, S n = 475 Wm −2 is the net shortwave flux, a = 0.725 and g = 0.995 are the emissivities of the 235 atmosphere and the surface, respectively, σ is the Stefan-Boltzmann constant, T a is the air temperature averaged over the lowest two model levels and T s is the surface temperature. The remaining components of the surface energy balance-the surface heat and moisture fluxes and the ground heat flux-and the resulting surface temperature T s are calculated with the NOAH land surface model (Tewari et al., 2004). The surface layer is parametrized with the revised MM5 similarity theory scheme (Jiménez et al., 2012). The soil type is sandy loam and the roughness length is 0.02 m. With these settings, the spatially and temporally 240 averaged sensible heat flux is roughly 150 Wm −2 , similar to Schmidli (2013).
The model is initialized at rest with the lapse rate Γ = 3 K km −1 . Random initial perturbations of potential temperature drawn from a uniform distribution between -0.5 and 0.5 K are added to the lowest five model levels. The setup leads to negligible latent heat fluxes and a very dry atmosphere, therefore moist processes are neglected. Due to the small domain size and zero background wind, also Coriolis force effects are not taken into account. 245 We calculate full θ-tendencies and decompose them into resolved turbulence, subgrid-scale turbulence and mean advective. The averaging in Eq. 16-18 is over 30 minutes and in the y-direction. Due to the y-averaging and the periodic boundary conditions, the flux derivatives in the y-direction are almost zero and therefore not shown. The budget calculation is carried out with the WRFlux procedure (Eq. 13), the two alternative formulations (Eq. 14 and Eq. 15) and with 2nd-order instead of the 3rd and 5th-order advection used by the model.
The budget components are divided by mean density to obtain tendencies of the form ∂tρθ ρ with units Kelvin per second. This shall not be confused with tendencies in advective form ∂ t θ .
We quantify the budget closure with the root-mean-square error of the sum of all forcing terms f with respect to the actual model tendency t normalized by the standard deviation of t: The averaging is over all gridpoints and 30-min averaging intervals. Following Chen et al. (2020), we also compute the 99thpercentile of the absolute residual scaled by the 99th percentile of the absolute tendency:

Cross-valley circulation
We start with a short overview of the individual heat budget components in the example simulation for the averaging period   Fig. 1f). Above the slope wind layer, vertical turbulent entrainment leads to cooling (Fig. 1b).
Above the ridge, a convective core develops that transports heat from the surface to higher levels. Lateral turbulent entrainment cools the convective core and warms the surrounding air (Fig. 1a). The mean advective tendency shows regions of horizontal θ-flux divergence on the slope and convergence on the ridge, and vice versa for the vertical (Fig. 1d and e). These regions 270 essentially coincide with regions of mass convergence and divergence (not shown). The scale of the horizontal and vertical advective tendencies is three orders of magnitude larger than for the corresponding turbulence tendencies ( Fig. 1a and b).
However, the respective sums of the horizontal and vertical components are of comparable magnitude ( Fig. 1c and f). Close to the surface, the net mean advective tendency shows cooling of the slope wind layer by the mean upslope wind and warming of the convective core due to horizontal mass convergence (Fig. 1f). The former is weaker than the turbulent heating, while 275 the latter is stronger than the turbulent cooling, leading to net warming close to the surface (Fig. 2a). Away from the surface, the mean advective tendency leads to cooling zones that propagate with time from the ridge towards the valley center. After 4 hours this cooling zone spans almost the whole domain in the horizontal direction ( Fig. 1f and 2a).
The total tendency in the mass-based terrain-following coordinate system shows somewhat different structures with stronger warming throughout the domain (Fig. 2b). The only difference between the total tendencies in the terrain-following and the 280 Cartesian formulation is the second term on the left-hand side in Eq. 11, which accounts for the height of the vertical levels being time-dependent. As we can see, this term has a considerable impact and is thus needed to close the budget in Eq. 11. In  Figure 2. Cross-sections of total θ-tendency for the averaging period between 3.5 and 4 h after initialization for the Cartesian (left-hand side of Eq. 11 divided by ρ ) and the terrain-following coordinate system (left-hand side of Eq. 6 divided by zηρ ). contrast, in the alternative form of the equation (Eq. 8), the correction term for the time derivative is almost negligible. This shows that the correction terms in Eq. 11 and 8 are conceptually different as mentioned in Sect. 2.1.
Since we use the equations in flux-form, Fig. 1 and 2, in general, cannot be compared to Schmidli (2013), who used the 285 advective form. However, as Schmidli (2013) points out, under the Boussinesq approximation the total turbulence tendency is equivalent in both formulations. In fact, the total turbulence tendency in Fig. 1c is of comparable magnitude and shows very similar spatial patterns as the one in Fig. 6c in Schmidli (2013).
As shown above, a budget equation typically consists of large competing forcing terms that add up to a relatively small total tendency. To illustrate this, instead of looking at the decomposition into total turbulence and mean advective tendencies as in 290 Fig. 1, we consider the two budget components as they are calculated by the model (not decomposed): resolved advection and subgrid-scale diffusion. The horizontal and vertical components of the resolved advection close to the surface are on the order of −1 K s −1 and +1 K s −1 , respectively. Their sum is much smaller, on the order of −10 −2 K s −1 . The subgrid-scale diffusion is on the order of +10 −2 K s −1 . Adding the resolved advection and the subgrid-scale diffusion leads to a total tendency on the order of +10 −5 K s −1 . When adding the large forcing terms, approximations in the budget calculation can lead to considerable 295 errors, as we will demonstrate in the following.

Comparison of budget calculation methods
We compare the budget obtained with WRFlux (Eq. 13) with several alternative forms. The first alternative uses 2nd-order advection in Eq. 13 instead of the advection order that is consistent with the model (3rd/5th-order). The differences between using 2nd order and the consistent advection order are largest on the ridge, where the vertical velocities are largest. The 300 horizontal and vertical components of the total θ-tendency (resolved + subgrid-scale) in Fig. 3 are both very close for the two   Figure 4 shows profiles of the resolved turbulence θ-tendency over the slope at x = −5 km for the three different formulations. The calculation of the vertical component is identical for all three formulations. The formulation in Eq. 15, which uses the consistently discretized θ in the horizontal correction, yields very similar profiles as the reference WRFlux procedure 310 (Eq. 13). In contrast, the formulation in Eq. 14, in which the horizontally destaggered and then vertically staggered horizontal flux is used in the corrections, results in considerable errors close to the surface.
To quantify the differences, we plot the sum of all forcing terms for each budget calculation method against the actual model tendency and compute the normalized root-mean-square error (NRMSE, Eq. 21) in Fig. 5. To avoid averaging out large errors, we drop the spatial averaging for this plot and only use temporal averaging. For the WRFlux procedure (Eq. 13), the points The lower row shows the results for the alternative formulations in Eq. 14 and Eq. 15. For this plot, the data is only averaged temporally, not spatially. The color code indicates the vertical levels of the model gridpoints with a focus on the lowest five levels. The gray line is the 1:1 line that signifies a perfectly closed budget.
lie close to the 1:1 line, indicating a good budget closure, quantified with an NRMSE of 9.17 · 10 −3 . The NRMSE increases by about one order of magnitude when using Eq. 15 and by about two orders when using Eq. 14 or 2nd-order advection. The errors are largest at the lowest vertical levels. For water vapor mixing ratio, the NRMSE of WRFlux is about 6 times smaller and for the windspeed components, it is about 15 times smaller (Table 1). For these variables, the other budget calculation methods lead to NRMSE values that are two to three orders of magnitude worse than the one of WRFlux. 320 We also tested two other approximations: using WRF's prognostic vertical velocity in the resolved vertical flux instead of the one recalculated with Eq. 9 and not including the density in the time-averaging of the total resolved flux (left-hand side of Eq. 17). The effect on the budget closure is moderate. For both of these approximations, the NRMSE score is increased by about one order of magnitude.
percentile of the absolute tendency (Eq. 22). For potential temperature, WRFlux reaches a value of r 99th ≈ 1.2%. For horizontal momentum (u windspeed), we reach a score of r 99th ≈ 0.07 %, similar to the value of 0.1 % that Chen et al. (2020) state for their simulations.

Conclusions
We developed a computational method to accurately diagnose the advective and turbulence components of the budgets of prog-330 nostic variables in a numerical weather prediction model. The method is based on a numerically consistent implementation of the transformation from a coordinate system with a generalized vertical coordinate, such as a terrain-following coordinate system, to the Cartesian coordinate system. The partitioning of the advective tendency into horizontal and vertical components is different in the two coordinate systems and thus the coordinate transformation is helpful when investigating the horizontal and vertical components separately. We illustrated this by assessing the local heat budget in a simulation of a convective 335 boundary layer over an idealized 2D mountain ridge. The slope flow layer is subject to vertical resolved and subgrid-scale turbulent heating from the ground and turbulent cooling due to vertical entrainment. Close to the surface, the sum of the potential temperature tendencies due to resolved horizontal and vertical advection is about two orders of magnitude smaller than the individual components. Adding the subgrid-scale diffusion yields a final tendency which is another three orders of magnitude smaller.

340
The circumstance of large and counteracting budget components adding up to a relatively small total tendency makes the budget calculation sensitive to approximations. While the sum of all forcing terms in WRFlux agrees to very high precision with the actual model tendency, we could show that approximations based on a lower-order advection operator or a numerically inconsistent formulation of the coordinate transformation lead to large residuals in the budget and noticeable differences in the tendency profiles. When looking at cross-section plots, the differences between the budget calculation methods are hardly 345 noticeable. Nevertheless, a budget analysis tool that yields large residuals is unreliable. In general, if the residual is large, we do not know whether the individual forcing terms are more or less reliable and only the sum is erroneous or whether also the forcing terms are not trustworthy due to approximations or software bugs. Therefore, a closed budget as achieved by WRFlux is essential. This requires the budget calculations to be consistent with the model numerics.
WRFlux expands the approach of Chen et al. (2020) by the computation of resolved turbulence tendencies and the transforma- the evolution of thermal updrafts in mountainous terrain that are subject to lateral and vertical turbulent entrainment The Jacobian matrix of the transformation from the Cartesian coordinate system x = (x 0 , x 1 , x 2 , x 3 ) = (t, x, y, z) to a coordinate system ξ = (τ, ξ 1 , ξ 2 , ξ 3 ) = (t, x, y, η) with generalized vertical coordinate η reads: which yields the Jacobian determinant |J| = ∂ η z.

370
The inverse of J is given by: