Numerically consistent budgets of potential temperature, momentum, and moisture 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 2 to 3 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 and, in comparison to existing similar tools, incorporates new capabilities. WRFlux transforms the budget equations from the terrain-following grid of the model to the Cartesian coordinate system, permitting a simplified interpretation of budgets obtained from simulations over non-uniform orography. WRFlux also decomposes the resolved advection into mean advective and resolved turbulence components, which is useful in the analysis of large-eddy simulation output. 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 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 1 to 3 orders of magnitude larger than with WRFlux.

1 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 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 result in an unclosed budget with a large residual, even if the relative error of each forcing term is small.

For the WRF model (Skamarock and Klemp, 2008), there have been several attempts to achieve precise budget calculations, e.g., Lehner (2012), Moisseeva and Steyn (2014), Potter et al. (2018), and Chen et al. (2020). Chen et al. (2020) developed a momentum and potential temperature budget analysis tool for the WRF model. 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 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 that it is numerically equivalent to the respective flux form (Xue and Lin, 2001), the advective form may hinder interpretation (Lee et al., 2004): if the fluxes on two opposing sides of a grid box are equal, the tendency component is zero for the flux form but not necessarily for the advective form of the equations.

By design, the budget computation method of Chen et al. (2020) cannot discriminate between tendencies caused by resolved-scale and subgrid-scale turbulence. This is important, especially when doing large-eddy simulations, which partially resolve 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 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 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.

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 us 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 to other, more simplified budget computation approaches.

2 Theory

2.1 Conservation equation transformations

The flux-form conservation equation for a variable ψ in the Cartesian coordinate system $\mathit{x}=(x_{\mathrm{0}},x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})=(t,x,y,z)$ reads

$$\begin{array}{}\text{(1)}& {\partial }_t\left(\mathit{\rho }\mathit{\psi }\right)=\sum _{i=\mathrm{1}}^{\mathrm{3}}-{\partial }_{x_i}\left(\mathit{\rho }{u}_i\mathit{\psi }\right)+S,\end{array}$$

where ρ is the air density, u_i is the ith component 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 ψ.

Equation (1) can be transformed from the Cartesian coordinate system to general curvilinear coordinates $\mathit{\xi }=(\mathit{\tau },{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}},{\mathit{\xi }}_{\mathrm{3}})$. Details about coordinate transformations, especially concerning coordinate systems with a generalized vertical coordinate, can be found for instance in Kasahara (1974), Pielke (1984), Byun (1999), and Liseikin (2010). We only give a short summary here.

The transformation of Eq. (1) yields

$$\begin{array}{}\text{(2)}& {\partial }_{\mathit{\tau }}\left(\mathit{\rho }\mathrm{|}\mathbf{J}\mathrm{|}\mathit{\psi }\right)=\sum _{i=\mathrm{1}}^{\mathrm{3}}-{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }\mathrm{|}\mathbf{J}\mathrm{|}{\mathit{\nu }}^i\mathit{\psi }\right)+\mathrm{|}\mathbf{J}\mathrm{|}S,\end{array}$$

where |J| is the determinant of the 4 × 4 Jacobi matrix of the transformation, ${\mathbf{J}}_{ij}={\partial }_{{\mathit{\xi }}_j}{x}_i$, and ${\mathit{\nu }}^i={\mathbf{J}}_{ij}^{-\mathrm{1}}{u}_j=u_j{\partial }_{x_j}{\mathit{\xi }}_i$ is the contravariant velocity in the new coordinate system with $\mathit{u}=(\mathrm{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 $\mathit{\xi }=(t,x,y,\mathit{\eta })$ 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

$$\begin{array}{}\text{(3)}& \mathrm{|}\mathbf{J}\mathrm{|}={\partial }_{\mathit{\eta }}z\equiv z_{\mathit{\eta }}.\end{array}$$

and the contravariant velocity vector is

$$\begin{array}{}\text{(4)}& \mathit{\nu }=(\mathrm{1},u,v,\dot{\mathit{\eta }})\end{array}$$

with

$$\begin{array}{}\text{(5)}& \dot{\mathit{\eta }}={\partial }_t\mathit{\eta }+u{\partial }_x\mathit{\eta }+v{\partial }_y\mathit{\eta }+w{\partial }_z\mathit{\eta }=:\mathit{\omega }.\end{array}$$

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 ${\mathit{\mu }}_{\mathrm{d}}=-{\mathit{\rho }}_{\mathrm{d}}g{z}_{\mathit{\eta }}$ in the model equations, where ρ_d is the dry-air density and g the acceleration due to gravity. All prognostic variables in WRF are coupled, i.e., multiplied, with μ_d.

Inserting Eqs. (3) and (4) into Eq. (2) yields

$$\begin{array}{}\text{(6)}& {\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\psi }\right)=\sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{u}_i\mathit{\psi }\right)\right]-{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\omega }\mathit{\psi }\right)+z_{\mathit{\eta }}S.\end{array}$$

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 τ, ξ₁, and ξ₂ in the equations, even though $(\mathit{\tau },{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}})=(t,x,y)$.

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

$$\begin{array}{}\text{(7)}& {\partial }_{{\mathit{\xi }}_i}A={\displaystyle \frac{\partial A}{\partial x_j}}{\displaystyle \frac{\partial x_j}{\partial {\mathit{\xi }}_i}}={\partial }_{x_i}A+z_{x_i}{\partial }_{z}A\end{array}$$

with $i=\mathrm{0},\mathrm{1},\mathrm{2}$, $j=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$,and $z_{x_i}:={\partial }_{{\mathit{\xi }}_i}z$.

Using Eq. (7) in Eq. (1) yields

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\partial }_{\mathit{\tau }}\left(\mathit{\rho }\mathit{\psi }\right)-z_t{\partial }_z\left(\mathit{\rho }\mathit{\psi }\right)=\\ & =\sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{u}_i\mathit{\psi }\right)+z_{x_i}{\partial }_z\left(\mathit{\rho }{u}_i\mathit{\psi }\right)\right]-{\partial }_z\left(\mathit{\rho }w\mathit{\psi }\right)+S.\end{array}\end{array}$$

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.

As will be pointed out in Sect. 2.4, Eq. (8) is not ideal for budget closure because the contained derivative terms cannot be discretized using numerical methods consistent with those for the governing equation (Eq. 6) in WRF. Therefore, we search for an alternative formulation which is equivalent to Eq. (8) and develop a discretization that is numerically consistent with Eq. (6). Starting with Eq. (6), we first replace ω with w using

$$\begin{array}{}\text{(9)}& w={\displaystyle \frac{\mathrm{d}z}{\mathrm{d}t}}={\mathbf{J}}_{\mathrm{3}j}{\mathit{\nu }}^j=z_t+z_{x}u+z_{y}v+z_{\mathit{\eta }}\mathit{\omega }.\end{array}$$

This equation is analogous to the geopotential tendency equation in WRF.

Solving for z_ηω, inserting in Eq. (6) and rearranging leads to

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\psi }\right)-{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{z}_t\mathit{\psi }\right)=\\ & =\sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{u}_i\mathit{\psi }\right)+{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{z}_{x_i}{u}_i\mathit{\psi }\right)\right]\\ & -{\partial }_{\mathit{\eta }}\left(\mathit{\rho }w\mathit{\psi }\right)+z_{\mathit{\eta }}S.\end{array}\end{array}$$

Dividing by z_η finally yields

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& z_{\mathit{\eta }}^{-\mathrm{1}}{\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\psi }\right)-{\partial }_z\left(\mathit{\rho }{z}_t\mathit{\psi }\right)=\\ & =\sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-z_{\mathit{\eta }}^{-\mathrm{1}}{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{u}_i\mathit{\psi }\right)+{\partial }_z\left(\mathit{\rho }{z}_{x_i}{u}_i\mathit{\psi }\right)\right]\\ & -{\partial }_z\left(\mathit{\rho }w\mathit{\psi }\right)+S.\end{array}\end{array}$$

Using the product rule and the commutativity of partial derivatives, one can show that Eq. (11) is mathematically equivalent to Eq. (8). For example, the horizontal flux divergence term in Eq. (10) can be expressed as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& -{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{u}_i\mathit{\psi }\right)+{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{z}_{x_i}{u}_i\mathit{\psi }\right)=\\ & =-z_{\mathit{\eta }}{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{u}_i\mathit{\psi }\right)-\mathit{\rho }{u}_i\mathit{\psi }{\partial }_{{\mathit{\xi }}_i}{z}_{\mathit{\eta }}+z_{x_i}{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{u}_i\mathit{\psi }\right)+\mathit{\rho }{u}_i\mathit{\psi }{\partial }_{\mathit{\eta }}{z}_{x_i}=\\ & =-z_{\mathit{\eta }}{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{u}_i\mathit{\psi }\right)+z_{x_i}{\partial }_{\mathit{\eta }}\left(\mathit{\rho }{u}_i\mathit{\psi }\right).\end{array}\end{array}$$

Dividing Eq. (12) by z_η gives the same expression of the horizontal flux divergence term as in Eq. (8). The left-hand side of Eq. (11) can be transformed analogously.

The correction terms in Eq. (11) (second term on the left-hand side and second term in square brackets) are conceptually 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.

Instead of Eq. (8), we select Eq. (11) as the budget equation because the coordinate metric z_η appears within the derivatives as in the WRF governing equation (Eq. 6), and so the associated budget analysis can be closed more precisely, consistent with the model dynamics (see Sect. 2.4). 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 model.

2.2 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 ${\mathit{\theta }}_{\mathrm{p}}=\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}$ with θ₀ = 300 K. Based on this decomposition, the potential temperature equation can be split up into advection of the perturbation and of the constant base state:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{z}_{\mathit{\eta }}S& ={\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\theta }\right)-\mathrm{\nabla }\cdot \left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\nu }\mathit{\theta }\right)=\\ & ={\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}{\mathit{\theta }}_{\mathrm{p}}\right)-\mathrm{\nabla }\cdot \left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\nu }{\mathit{\theta }}_{\mathrm{p}}\right)+{\mathit{\theta }}_{\mathrm{0}}\left[{\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\right)-\mathrm{\nabla }\cdot \left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\nu }\right)\right],\end{array}\end{array}$$

with the contravariant velocity ν. Due to the continuity equation, the last term on the right-hand side of Eq. (13) is identically zero. Equation (13) can be used to compute the components of the full-θ tendency with high numerical accuracy.

2.3 Advective form

So far we looked at the budget equations in flux form. This form corresponds to the tendency of ρψ. Often, however, we are interested in the tendency of ψ itself. This is particularly relevant for θ, which usually has a tendency opposing the one of ρ. The advective form in the terrain-following coordinate system can be obtained from Eq. (6) using the continuity equation and reads

$$\begin{array}{}\text{(14)}& {\partial }_{\mathit{\tau }}\mathit{\psi }=-\sum _{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\nu }}_i{\partial }_{{\mathit{\xi }}_i}\mathit{\psi }+S.\end{array}$$

To compute the advective form in a numerically consistent way, we can use the components of the flux-form equation and the mass tendencies from Eq. (13). The left-hand side of Eq. (14) can be computed as

$$\begin{array}{}\text{(15)}& {\partial }_{\mathit{\tau }}\mathit{\psi }=(\mathit{\rho }{z}_{\mathit{\eta }}{)}^{-\mathrm{1}}\left[{\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\psi }\right)-\mathit{\psi }{\partial }_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\right)\right]\end{array}$$

and the components of the right-hand side as

$$\begin{array}{}\text{(16)}& {\mathit{\nu }}_i{\partial }_{{\mathit{\xi }}_i}\mathit{\psi }=(\mathit{\rho }{z}_{\mathit{\eta }}{)}^{-\mathrm{1}}\left[{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{\mathit{\nu }}_i\mathit{\psi }\right)-\mathit{\psi }{\partial }_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{\mathit{\nu }}_i\right)\right].\end{array}$$

In the Cartesian coordinate system correction terms for the mass tendency components are introduced analogously to the correction terms in Eq. (11).

2.4 Discretization

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 Eq. (11) as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& z_{\mathit{\eta }}^{-\mathrm{1}}{\mathit{\delta }}_{\mathit{\tau }}\left(\mathit{\rho }{z}_{\mathit{\eta }}\mathit{\psi }\right)-{\mathit{\delta }}_z\left(\mathit{\rho }{z}_t{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)=\\ =& \sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-z_{\mathit{\eta }}^{-\mathrm{1}}{\mathit{\delta }}_{{\mathit{\xi }}_i}\left(\mathit{\rho }{z}_{\mathit{\eta }}{u}_i{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{x_i}\right)+{\mathit{\delta }}_z\left(\mathit{\rho }{z}_{x_i}{\stackrel{\mathrm{‾}}{u_i}}^{x_{i}z}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)\right]\\ & -{\mathit{\delta }}_z\left(\mathit{\rho }w{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)+S.\end{array}\end{array}$$

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 positive-definite, monotonic, and weighted essentially non-oscillatory options.

For Eq. (17) 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.

Although the momentum variables are staggered differently from the thermodynamic variables, their discretized equations can be derived analogously. We do not state them here for brevity.

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 precise budget equation derived above. This is because when applying the product rule in Eq. (12) to transform Eq. (11) into Eq. (8), the second and fourth term in the second line of Eq. (12) only cancel out analytically but not numerically: the flux ρu_iψ in the former originates from a horizontal derivative and thus must be discretized as ${\stackrel{\mathrm{‾}}{\mathit{\rho }{u}_i{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{x_i}}}^{x_i}$, while the latter comes from a vertical derivative and is discretized as ${\stackrel{\mathrm{‾}}{\mathit{\rho }{\stackrel{\mathrm{‾}}{u_i}}^{x_{i}z}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z}}^z$. To demonstrate the inconsistency, we compare Eq. (17) 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 and vertically to the grid of the vertical flux:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\mathit{\delta }}_{\mathit{\tau }}\left(\mathit{\rho }\mathit{\psi }\right)-z_t{\mathit{\delta }}_z\left(\mathit{\rho }{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)=\\ =& \sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-{\mathit{\delta }}_{{\mathit{\xi }}_i}\left(\mathit{\rho }{u}_i{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{x_i}\right)+{\stackrel{\mathrm{‾}}{z_{x_i}}}^{x_i}{\mathit{\delta }}_z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }{u}_i{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{x_i}}}^{x_{i}z}\right)\right]\\ & -{\mathit{\delta }}_z\left(\mathit{\rho }w{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)+S.\end{array}\end{array}$$

This is analogous to the implementation of subgrid-scale diffusion in WRF.

The second one adopts a different discretization in the horizontal correction term that is closer to the one in Eq. (17):

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\mathit{\delta }}_{\mathit{\tau }}\left(\mathit{\rho }\mathit{\psi }\right)-z_t{\mathit{\delta }}_z\left(\mathit{\rho }{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)=\\ =& \sum _{i=\mathrm{1}}^{\mathrm{2}}\left[-{\mathit{\delta }}_{{\mathit{\xi }}_i}\left(\mathit{\rho }{u}_i{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{x_i}\right)+{\stackrel{\mathrm{‾}}{z_{x_i}}}^{x_i}{\mathit{\delta }}_z\left(\mathit{\rho }{\stackrel{\mathrm{‾}}{u_i}}^{x_{i}z}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)\right]\\ & -{\mathit{\delta }}_z\left(\mathit{\rho }w{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^z\right)+S.\end{array}\end{array}$$

The impact of the approximate budget calculations (Eqs. 18 and 19) is discussed in Sect. 4.3.

2.5 Flux averaging and decomposition

After introducing the discretizations of the precise budget equation (Eq. 17) and of the alternative formulations (Eqs. 18 and 19) we now turn to the decomposition of the fluxes into mean advective, resolved turbulent, and subgrid-scale turbulent components. This decomposition often provides valuable insights for large-eddy simulations. It affects several terms in the conservation equations: the vertical flux, the two horizontal flux components and corresponding correction terms, and the subgrid-scale fluxes that are part of S. The subgrid-scale fluxes are taken as computed by the subfilter-scale model or planetary boundary layer scheme. Here, we show the decomposition of the resolved fluxes into mean advective and resolved turbulent. This decomposition requires averaging the fluxes over time and/or space (e.g.,  Schmidli, 2013). The averaging is regarded as an approximation of an ensemble average. Means and perturbations are defined by

$$\begin{array}{}\text{(20)}& \stackrel{\mathrm{̃}}{\mathit{\psi }}={\displaystyle \frac{〈\mathit{\rho }\mathit{\psi }〉}{〈\mathit{\rho }〉}},{\mathit{\psi }}^{\prime \prime }:=\mathit{\psi }-\stackrel{\mathrm{̃}}{\mathit{\psi }}.\end{array}$$

〈ψ〉 denotes the time and/or spatial block average, $\stackrel{\mathrm{̃}}{\mathit{\psi }}$ is the density-weighted average and ψ^′′ the perturbation thereof.

The decomposition of the resolved flux then reads

$$\begin{array}{}\text{(21)}& 〈\mathit{\rho }{u}_i\mathit{\psi }〉=〈\mathit{\rho }〉\stackrel{\mathrm{̃}}{u_i}\stackrel{\mathrm{̃}}{\mathit{\psi }}+〈\mathit{\rho }{u}_i^{\prime \prime }{\mathit{\psi }}^{\prime \prime }〉\text{ for }i=\mathrm{1},\mathrm{2},\mathrm{3},\end{array}$$

with the total resolved 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 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, the impact on the mean advective and resolved turbulent components is also hardly noticeable.

The correction flux used in the horizontal corrections in Eq. (11) can be decomposed as

$$\begin{array}{}\text{(22)}& 〈\mathit{\rho }{Z}_i\mathit{\psi }〉=〈\mathit{\rho }〉\stackrel{\mathrm{̃}}{Z_i}\stackrel{\mathrm{̃}}{\mathit{\psi }}+〈\mathit{\rho }{Z}_i^{\prime \prime }{\mathit{\psi }}^{\prime \prime }〉\text{ for }i=\mathrm{1},\mathrm{2},\end{array}$$

with $Z_i=z_{x_i}{u}_i$.

3 Implementation

We implemented the theoretical framework of the previous section in a diagnostic package for WRF: WRFlux. The main features of WRFlux are as follows.

• Budget components are retrieved for potential temperature, water vapor mixing ratio, and momentum, including tendencies from the acoustic time step, subgrid-scale diffusion (from all available subfilter-scale models and planetary boundary layer schemes), physical parameterizations, and numerical diffusion and damping.

• The subgrid-scale and resolved fluxes and all budget components except for advection are averaged in time during model integration over a user-specified time window. The optional spatial averaging and computation of advective tendencies with decomposition into mean advective and resolved turbulent is done in the post-processing. The resolved turbulent component is calculated using Eq. (21).

• 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. The recalculated w is only used in diagnosing the vertical advection term in the Cartesian coordinate system, not as the prognostic variable in the conservation equation for w. 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 geopotential.¹ 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 ${\mathit{\theta }}_{\mathrm{p}}=\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}$ and full θ, the last term in Eq. (13) needs to vanish. That means we have to close the mass continuity equation. This is difficult since the continuity equation in WRF is not integrated explicitly; rather μ_d is diagnosed from its definition after vertically integrating the mass divergence. To close the continuity equation anyway, we calculate the temporal term and horizontal divergence terms of the continuity equation explicitly and take the vertical term as the residual. Using the residual only has a marginal effect on the vertical component but avoids large residuals in the final θ budget.

• The mean advective tendencies, the total advective tendencies, and the final model tendency can be transformed to the advective form using the components of the time-averaged continuity equation in Eqs. (15) and (16). This can be done in the post-processing without changing the online part. The resolved turbulence tendency is left in flux form.

• 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.

• Before each update of WRFlux, an automated test suite is carried out that checks the output of WRFlux for consistency using idealized test simulations with a large number of different namelist settings. Details about the tests can be found in the documentation of WRFlux. The latest version of WRFlux, version 1.3, is based on WRF-ARW version 4.3. WRFlux is easy to install, and new releases of WRF are continuously integrated. The post-processing tool is written in Python.

4 Example of application and the effect of approximations

4.1 Simulation design

We demonstrate the capabilities of WRFlux with a simulation of the diurnal evolution of the convective boundary layer over mountainous terrain using WRFlux version 1.2.1. The model setup, i.e., the initial conditions, terrain specification, grid spacing, land surface properties, and the choice of the subfilter-scale model, follows Schmidli (2013).

WRF's hybrid terrain-following coordinate is used. There are 140 vertical levels with a vertical grid spacing ranging approximately from 8 m at the surface to 50 m at the model top at a height of 5 km. The horizontal grid spacing is Δx = 50 m. The model time step is Δt = 1 s. The domain size is 20 km in the x (cross-mountain) and 10 km in the y direction (along-mountain). The boundary conditions are periodic in both directions. The topography is a periodic two-dimensional cosine valley with a flat valley bottom and flat mountain ridge:

$$\begin{array}{}\text{(23)}& h\left(x\right)=\left\{\begin{array}{cc}{h}_{\mathrm{m}}& \mathrm{|}x\mathrm{|}\le x_{\mathrm{1}}\\ h_{\mathrm{m}}\left\{\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\cos \left[\frac{\mathit{\pi }}{x_{\mathrm{2}}-x_{\mathrm{1}}}\left(\mathrm{|}x\mathrm{|}-x_{\mathrm{1}}\right)\right]\right\}& x_{\mathrm{1}}<\mathrm{|}x\mathrm{|}\le x_{\mathrm{2}}\\ \mathrm{0}& \mathrm{|}x\mathrm{|}>x_{\mathrm{2}}\end{array}\right.,\end{array}$$

with the ridge height h_m = 1500 m and x₁ = 0.5 km, x₂ = 9.5 km. The valley-to-valley distance is thus 20 km, equal to the domain width. The distance in the x direction is defined to be 0 at the center of the domain.

Implicit Rayleigh damping (Klemp et al., 2008) 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 fifth-order in the horizontal and third-order in the vertical. Subgrid-scale diffusion follows 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

$$\begin{array}{}\text{(24)}& R_{\mathrm{n}}=S_{\mathrm{n}}+{\mathit{ϵ}}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathit{ϵ}}_{\mathrm{g}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},\end{array}$$

where R_n is the net radiation, S_n = 475 W m⁻² is the net shortwave flux, ϵ_a=0.725 and ϵ_g=0.995 are the emissivities of the 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 averaged sensible heat flux is roughly 150 W m⁻², similar to Schmidli (2013).

The model is initialized at rest with the lapse rate Γ = 3 K km⁻¹ and a constant relative humidity of 40 % and run for 4 h. 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, Coriolis force effects are not taken into account either.

Since no microphysics scheme is activated and the simplified radiation scheme only affects the surface energy balance, the heat budget in the atmosphere only consists of resolved advection and subgrid-scale diffusion. For general applications, other grid-resolved and parameterized physics terms are possible and categorized as additional budget components. We calculate full-θ tendencies and decompose them into resolved turbulence, subgrid-scale turbulence, and mean advective. The averaging in Eqs. (20)–(22) is over 30 min and in the y direction. The averaging interval of 30 min is often used to compute turbulence statistics as it typically provides a good compromise between obtaining a large sample for the statistics while still being able to assume stationarity (Stiperski and Rotach, 2016). 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. 17), with the two alternative formulations (Eqs. 18 and 19), and with second-order instead of the third- and fifth-order advection used by the model.

Figure 1 Cross sections of total turbulence (trb = resolved + subgrid-scale turbulence, panels a–c) and mean advective (panels d–f) θ tendency components in the Cartesian coordinate system for the averaging period between 3.5 and 4 h after initialization. The calculation is based on Eq. (11), discretized according to Eq. (17), and decomposed with Eqs. (21) and (22). The horizontal (a, d) and vertical (b, e) components are the flux derivatives in the cross-mountain and vertical direction, respectively. Panels (c) and (f) show the sum of the horizontal and vertical components. The units of the color bar are denoted in each panel. Note that a different color scale is used in panels (d) and (e). The contour lines represent mean potential temperature with a spacing of 0.5 K. Panel (f) also shows averaged wind vectors.

The budget components are divided by mean density to obtain tendencies of the form $\frac{〈{\partial }_t\mathit{\rho }\mathit{\theta }〉}{〈\mathit{\rho }〉}$ with units of Kelvin per second. This should 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:

$$\begin{array}{}\text{(25)}& \text{NRMSE}=\sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{(t-f{)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{{\left(t-\stackrel{\mathrm{‾}}{t}\right)}^{\mathrm{2}}}}}}.\end{array}$$

The averaging is over all grid points and 30 min averaging intervals. Following Chen et al. (2020), we also compute the 99th percentile of the absolute residual scaled by the 99th percentile of the absolute tendency:

$$\begin{array}{}\text{(26)}& r^{\mathrm{99}\text{th}}={\displaystyle \frac{p_{\mathrm{99}}(\mathrm{|}t-f\mathrm{|})}{p_{\mathrm{99}}\left(\mathrm{|}t\mathrm{|}\right)}}.\end{array}$$

4.2 Cross-valley circulation

We start with a short overview of the individual heat budget components in the example simulation for the averaging period between 3.5 and 4 h after initialization.

Figure 2 Profiles of turbulence θ tendency (subgrid-scale, resolved, and their sum) on the ridge at x=0, for the averaging period between 3.5 and 4 h after initialization.

Figure 3 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. 10) and the terrain-following coordinate system (left-hand side of Eq. 6), both divided by 〈z_ηρ〉.

Figure 4 Profiles of total θ tendency (resolved + subgrid-scale) on the ridge at x=0, for the averaging period between 3.5 and 4 h after initialization resulting from flux derivatives in the x (blue, lower x axis) and z (orange, lower x axis) directions and their sum (green, upper x axis) calculated with second-order advection (dashed) and third- (vertical) and fifth-order (horizontal) advection, consistent with the numerical model (solid).

Figure 1 shows cross sections of the total turbulence (resolved + subgrid-scale) and mean advective components of the heat budget in the Cartesian coordinate system (Eq. 17). The dynamics are driven by the surface sensible heat flux; vertical turbulent flux convergence in the layer close to the surface (Fig. 1b) induces upslope winds and compensatory return flows aloft and in the valley center (wind arrows in 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). Figure 2 shows the resolved and subgrid-scale turbulence components above the ridge separately. The subgrid-scale component is only relevant close to the ground. There, it causes a positive tendency (flux convergence) by diffusing the sensible heat flux from the ground, while resolved turbulent eddies cause a negative tendency (flux divergence) by transporting the heat further up. The two contributions largely balance out, but the subgrid-scale warming is of slightly lower magnitude than the cooling by the resolved turbulence. The resulting tendency from total turbulence is locally negative; this only occurs within the thermal plume at the ridge top, as visible in Fig. 1c. At all other locations along the slope, the heating by subgrid-scale heat flux convergence offsets the cooling operated by resolved turbulent transport. 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 essentially coincide with regions of mass convergence and divergence (not shown). The scale of the horizontal and vertical advective tendencies is 3 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 the latter is stronger than the turbulent cooling, leading to net warming close to the surface (Fig. 3a). 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 h this cooling zone spans almost the whole domain in the horizontal direction (Figs. 1f and 3a).

The total tendency in the mass-based terrain-following coordinate system shows somewhat different structures with stronger warming throughout the domain (Fig. 3b). The only difference between the total tendencies in the terrain-following and the 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; the coordinate layers expand as they are heated up. As we can see, this term has a considerable impact and is thus needed to close the budget in Eq. (11). In contrast, in the alternative form of the equation (Eq. 8), the correction term for the time derivative is almost negligible (not shown). With this formulation, Fig. 3a (with correction) and Fig. 3b (without correction) would look almost identical. This shows that the correction terms in Eqs. (11) and (8) are conceptually different, as mentioned in Sect. 2.1.

Since we use the equations in flux form, Figs. 1 and 3, in general, cannot be compared to Schmidli (2013), who used the 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 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 above the ridge are on the order of −1 and +1 K s⁻¹, respectively. Their sum is much smaller: on the order of −10⁻² K s⁻¹. The subgrid-scale diffusion is on the order of +10⁻² K s⁻¹. Adding the resolved advection and the subgrid-scale diffusion leads to a total tendency on the order of +10⁻⁵ K s⁻¹. When adding the large forcing terms, approximations in the budget calculation can lead to considerable errors, as we will demonstrate in the following.

4.3 Comparison of budget calculation methods

We compare the budget obtained with WRFlux (Eq. 17) with several alternative forms. The first alternative uses second-order advection in Eq. (17) instead of the advection order that is consistent with the model (third- and fifth-order). The differences between using second-order and the consistent advection order are largest on the ridge, where the vertical velocities are largest. The horizontal and vertical components of the total θ tendency (resolved + subgrid-scale) in Fig. 4 are both very close for the two calculation methods. But when adding these large and opposing components, the second-order calculation yields considerably different results close to the surface. Instead of the constant warming up to the entrainment layer, we see large oscillations of warming and cooling. The main errors derive from the vertical component (not shown). Only the deviation at 250 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ originates from errors in the horizontal turbulent entrainment. At other locations in the domain, where the up- and downdrafts are weaker, the differences are smaller.

Figure 5 Profiles of resolved turbulence θ tendency over the slope at x = −5 km, for the averaging period between 3.5 and 4 h after initialization resulting from flux derivatives in the x (blue) and z (orange) directions and their sum (green). The line styles indicate different formulations for the horizontal flux derivatives (explained in Sect. 2).

Figure 6 Scatterplots of the right-hand side (θ forcing) of Eq. (17), Eq. (17) with second-order advection, Eq. (18), and Eq. (19) against the respective left-hand side (θ tendency) for all model grid points and eight half-hourly (from initialization to 4 h later) values together with the corresponding NRMSE values (Eq. 25). For this plot, the data are only averaged temporally, not spatially. The color code indicates the vertical levels of the model grid points with a focus on the lowest five levels. The gray line is the 1:1 line that signifies a perfectly closed budget.

The different formulations for the horizontal flux derivatives (Sect. 2) differ significantly only where the grid elements are strongly tilted. Figure 5 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. (19), which uses the consistently discretized θ in the horizontal correction, yields very similar profiles as the reference WRFlux procedure (Eq. 17). In contrast, the formulation in Eq. (18), in which the horizontally de-staggered and then vertically staggered horizontal flux is used in the corrections, results in considerable errors close to the surface.

Table 1 NRMSE (Eq. 25) and r⁹⁹ (%, Eq. 26) values for all budget variables and budget calculation methods.

To quantify the differences, we plot the right-hand side (forcing) of Eq. (17), Eq. (17) with second-order advection, Eq. (18), and Eq. (19) against the respective left-hand side (tendency) and compute the normalized root-mean-square error (NRMSE, Eq. 25) in Fig. 6. To avoid averaging out large errors, we drop the spatial averaging for this plot and only use temporal averaging. For the WRFlux procedure (Eq. 17), the points lie close to the 1:1 line, indicating a good budget closure, quantified with an NRMSE of 9.17 × 10⁻³. The NRMSE increases by about 1 order of magnitude when using Eq. (19) and by about 2 orders when using Eq. (18) or second-order advection. The errors are largest at the lowest vertical levels. For the water vapor mixing ratio, the NRMSE of WRFlux is about 6 times smaller and for the wind speed components, it is about 15 times smaller (Table 1). For these variables, the other budget calculation methods lead to NRMSE values that are 2 to 3 orders of magnitude worse than the one of WRFlux.

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. 21). The effect on the budget closure is moderate. For both of these approximations, the NRMSE score is increased by about 1 order of magnitude.

To compare our results to Chen et al. (2020), we compute the 99th percentile of the absolute residual scaled by the 99th percentile of the absolute tendency (Eq. 26). For potential temperature, WRFlux reaches a value of r^99th≈ 1.2 %. For horizontal momentum (u wind speed), 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.

5 Conclusions

We developed a computational method to accurately diagnose the advective and turbulence components of the budgets of prognostic 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 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 2 orders of magnitude smaller than the individual components. Adding the subgrid-scale diffusion yields a total tendency that is another 3 orders of magnitude smaller.

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 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 the forcing terms are not trustworthy either 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 transformation of fluxes and flux divergence components to the Cartesian coordinate system. Possible applications of our budget analysis tool include the study of

• the reasons for the unclosed surface energy balance often reported in field studies (e.g.,  De Roo and Mauder, 2018) for which diagnostics in a layer close to the surface are required;

• the evolution of thermal updrafts in mountainous terrain that are subject to lateral and vertical turbulent entrainment (e.g.,  Kirshbaum, 2011, 2020);

• the exchange of heat and moisture between the boundary layer of a valley and the free troposphere (e.g.,  Rotach et al., 2015; Leukauf et al., 2015, 2017).

A conceivable extension for WRFlux is the inclusion of further budget variables, such as the mixing ratios of other water species or of a passive tracer.

Appendix A: Coordinate transformation

The Jacobian matrix of the transformation from the Cartesian coordinate system $\mathit{x}=(x_{\mathrm{0}},x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})=(t,x,y,z)$ to a coordinate system $\mathit{\xi }=(\mathit{\tau },{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}},{\mathit{\xi }}_{\mathrm{3}})=(t,x,y,\mathit{\eta })$ with generalized vertical coordinate η reads

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathbf{J}& {\displaystyle }=\left(\begin{array}{cccc}{\partial }_{\mathit{\tau }}t& {\partial }_{{\mathit{\xi }}_{\mathrm{1}}}t& {\partial }_{{\mathit{\xi }}_{\mathrm{2}}}t& {\partial }_{{\mathit{\xi }}_{\mathrm{3}}}t\\ {\partial }_{\mathit{\tau }}x& {\partial }_{{\mathit{\xi }}_{\mathrm{1}}}x& {\partial }_{{\mathit{\xi }}_{\mathrm{2}}}x& {\partial }_{{\mathit{\xi }}_{\mathrm{3}}}x\\ {\partial }_{\mathit{\tau }}y& {\partial }_{{\mathit{\xi }}_{\mathrm{1}}}y& {\partial }_{{\mathit{\xi }}_{\mathrm{2}}}y& {\partial }_{{\mathit{\xi }}_{\mathrm{3}}}y\\ {\partial }_{\mathit{\tau }}z& {\partial }_{{\mathit{\xi }}_{\mathrm{1}}}z& {\partial }_{{\mathit{\xi }}_{\mathrm{2}}}z& {\partial }_{{\mathit{\xi }}_{\mathrm{3}}}z\end{array}\right)=\text{(A2)}& {\displaystyle }& {\displaystyle }=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ {\partial }_{\mathit{\tau }}z& {\partial }_{{\mathit{\xi }}_{\mathrm{1}}}z& {\partial }_{{\mathit{\xi }}_{\mathrm{2}}}z& {\partial }_{\mathit{\eta }}z\end{array}\right),\end{array}$$

which yields the Jacobian determinant $\mathrm{|}\mathbf{J}\mathrm{|}={\partial }_{\mathit{\eta }}z$.

The inverse of J is given by

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathbf{J}}^{-\mathrm{1}}& {\displaystyle }=\left(\begin{array}{cccc}{\partial }_t\mathit{\tau }& {\partial }_x\mathit{\tau }& {\partial }_y\mathit{\tau }& {\partial }_z\mathit{\tau }\\ {\partial }_t{\mathit{\xi }}_{\mathrm{1}}& {\partial }_x{\mathit{\xi }}_{\mathrm{1}}& {\partial }_y{\mathit{\xi }}_{\mathrm{1}}& {\partial }_z{\mathit{\xi }}_{\mathrm{1}}\\ {\partial }_t{\mathit{\xi }}_{\mathrm{2}}& {\partial }_x{\mathit{\xi }}_{\mathrm{2}}& {\partial }_y{\mathit{\xi }}_{\mathrm{2}}& {\partial }_z{\mathit{\xi }}_{\mathrm{2}}\\ {\partial }_t{\mathit{\xi }}_{\mathrm{3}}& {\partial }_x{\mathit{\xi }}_{\mathrm{3}}& {\partial }_y{\mathit{\xi }}_{\mathrm{3}}& {\partial }_z{\mathit{\xi }}_{\mathrm{3}}\end{array}\right)=\text{(A4)}& {\displaystyle }& {\displaystyle }=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ {\partial }_t\mathit{\eta }& {\partial }_x\mathit{\eta }& {\partial }_y\mathit{\eta }& {\partial }_z\mathit{\eta }\end{array}\right).\end{array}$$