Reply on RC4

I suggest reconstructing the abstract to emphasize the unique features and applicability of WRFlux to differentiate this work from Chen et al. (2020), such as the retrieval of the resolved turbulence component, transformation to the easy-to-interpret Cartesian coordinate, etc., all of which are particularly useful for large-eddy simulations or simulations over nonuniform topography. Good point. To emphasize these features, we added a sentence and moved the last sentence to a more prominent place (L5-9). Introduction: It is important to mention that there have been some attempts to achieve a precise budget retrieval in the exact WRF model (Chen et al. 2020 is not the first), e.g., Lehner (2012, https://collections.lib.utah.edu/ark:/87278/s61n8fxw), Moisseeva and Steyn (2014, https://doi.org/10.5194/acp-14-13471-2014) and Potter et al. (2018, https://doi.org/10.1029/2018JD029427), etc. OK, we added these references in the introduction in L28-29. L35-38: I’m not sure how 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”? Can you please elaborate? We added a sentence there: “For instance, if the fluxes on two opposing sides of a grid box are equal, the flux form yields zero net tendency, while this is not necessarily the case for the advective form.” (L39-40) Lee et al. 2004 looked at the budget of a larger control volume and noted that this is due to “internal processes” in the control volume. We don’t want to explain this further in the paper, so we dropped the unclear mentioning of the “internal processes”. Chen et al. (2020) noted that a closed budget for w is challenging, partially because w equation is implicit, coupled with the geopotential tendency equation, and partially due to the variable’s inherent rapid variation on small scales. How does WRFlux overcome this issue for the w budget? It is also surprising to see that WRFlux performs the best for w compared to all the other variables and that the accuracy (in terms of NRMS and r99; Table 1) is much higher for the momentum than for the thermodynamic variables (θ and qv) in your case. Do you have any idea why? Is this result case-dependent or relevant to the physical design of WRFlux? Chen et al. found that closing the wbudget offline is difficult or impossible due to the mentioned issues, but the online budget calculation works fine for them. They also achieved a higher degree of budget closure for w than for horizontal momentum. Our budget calculations are also done online considering all budget components for w, including the important acoustic-step contributions. In contrast to Chen et al. (2020), our budget calculations are time-averaged. To check the effect of the time averaging on the budget closure, we did two small test simulations with and without time averaging. Without time averaging, the NRMSE scores in the terrain-following coordinate system become considerably better for qv and w and slightly worse for θ, u, and v. We did not investigate why the budget closure in the shown simulations is higher for momentum than for θ and qv. However, in some of our test simulations, it was the other way around, so it is definitely case-dependent (background wind speed, moisture content, physical parameterizations used, etc.). It would help readers follow Section 2 more easily if the connections between subsections are made more explicit and some terminologies are explained more precisely. For example, I find it confusing that “Eq. 11 is equivalent to Eq. 8” (L124-128), but “Eq. 8 cannot be closed numerically in contrast to Eq.11” (L128). This confusion may be related to how you define the term “close” here. Suppose you are referring to “budget closure”. In that case, my understanding is that it depends on the accuracy of the conducted budget calculation, not the format of the chosen equation to which the budget analysis is applied. Do you mean that because Eq. 11 is numerically consistent with the equation used in the model, the left-hand-side terms of Eq. 11 represent more accurately the actual tendency simulated by the model? And therefore, a budget analysis using Eq. 11 can achieve a higher degree of budget closure? Please clarify. Eq. 11 is similar to WRF’s governing equations because the coordinate metric zη appears within the derivatives. In the manuscript, we showed how the budget can be closed with Eq. 11. When using the product rule to transform Eq. 11 to Eq. 8, two terms cancel out analytically, but not numerically. This is because one of the canceling terms comes from a horizontal derivative and the other from a vertical derivative. Therefore, the right-hand side of Eq. 16 is not identical to the right-hand side of Eq. 14. Thus, we argue that Eq. 16 cannot be closed numerically in WRF. The left-hand sides of Eq. 14 and 16, instead, are almost identical in our simulations. To clarify this point, we added Eq. 12 and a bit of explanation in L167-170. Also, I suggest reminding readers at the beginning of Section 2.4 that Eq. 11 (i.e., Eq. 13 after discretization) will be used for the precise budget calculation. It would also be helpful to specify which terms in Eq. 11 can be decomposed into the mean advective and turbulent components (all righthand-side terms excluding S or all right-hand-side terms excluding only some partial S, i.e., the subgrid physics?) OK. We added a few introductory sentences in Section 2.4 with the suggested content (L184-189). L113-L114: “Numerical consistent” --> “Numerical consistency”. Also, placing this sentence (“Numerically... after discretization.”) here seems odd. You are neither showing the discretization nor discussing the budget closure in the following lines yet. That’s right. Nevertheless, in our opinion, this sentence is important for the logical train of thought. We added a reference to the discretization section in L116. L130: It is still unclear to me why recalculating the vertical velocity is desired. If Eqs. 6 and 9 are also used in WRF, doesn’t it make more sense to use WRF’s prognostic to be numerically consistent with the model? Unless this is relevant to the issue I mentioned in comment #4? Recalculating w is necessary to obtain a closed budget in the Cartesian coordinate system. Originally, there was a significant difference in the calculation of zηω in WRF‘s implementation of Eq. 6 and 9. We reduced this difference strongly by changing the discretization of Eq. 9 in WRF (L216-221, https://github.com/wrfmodel/WRF/pull/1338). Nevertheless, there is still a small difference remaining: In Eq. 6, zη is included in the coupled prognostic variable μdω, whereas in Eq. 9 it is calculated directly from the geopotential. To overcome this difference and obtain a closed budget, we recalculate w in a consistent way and use it in the vertical advection term. Whether this is appropriate is definitely debatable. But the effect on w itself is hardly noticeable. Only when adding up the large and counteracting budget components, the difference becomes noticeable. With the prognostic w, the budget closure (NRMSE) is about one order of magnitude worse (L356-359). We suppose the reviewer wanted to ask whether this is relevant to the issue mentioned in comment #5, not #4. These issues are not connected. The recalculation of w mentioned is only relevant for the transformation to the Cartesian coordinate system. It is not used as the prognostic variable in the conservation equation for w. We clarified this in L217-219. L191-194: Does this part follow exactly Chen et al. (2020)? The retrieval of subgrid-scale tendencies is different since we only extract the fluxes and then do the tendency calculation offline (as for the resolved fluxes). The retrieval of other forcing terms (physics, pressure gradient, damping, diffusion, Coriolis) was implemented very similarly to Chen et al. but independently from them. L195: “The fluxes and all... are averaged in time during model integration”. Averaged over the entire simulation? Or over a user-specified time window? If so, what are your suggested value and the physical reason behind it? Over a user-specified time window (we added this in L213). In this study, we use 30 min (L272) as is often used for Reynolds averaging (compromise between large sample for statistics and stationarity). We added a sentence with a reference on this in L272-274. L203-204: ”Since WRF does not actually solve the continuity equation...” This may be misleading. WRF does solve the continuity equation, just not in the same format as expressed in your Eq. (12), i.e., in terms of the 3D variable ρ but column dry air mass per unit area (2D μd). Note that μd is indeed advanced with time in the WRF dynamical solver. I suggest changing the sentence to something like “Since WRF does not directly solve for ρ but the integrated column dry air mass...” Actually, WRF's mass continuity equation in terms of μd (Eq. 2.12 in the technical notes version 4, http://www2.mmm.ucar.edu/wrf/users/docs/technote/v4_technote.pdf) does correspond to the rightmost term in our Eq. 13, except for a factor -g. In fact, μd and ρd are related by a diagnostic equation: μd = -g ρd zη. Ensuring that the last term of our Eq. 13 is exactly zero is difficult because the mass continuity equation in WRF is indeed solved but not integrated explicitly in the acoustic time steps, but rather μd is diagnosed from its definition (Eq. 3.20 in the technical notes) after vertically integrating the mass divergence (Eq. 3.22). We acknowledge that our original formulation was misleading, and we replaced it with the information given in the preceding sentence (L223-225). L244: “The setup leads to ... a very dry atmosphere, therefore moist processes are neglected.” Not sure if I overlook something but the initial setup of the moisture field is not mentioned? Yes. We added the specification of the moisture profile (constant RH=40%) in L263. L246: “We calculate full θ-tendencies and decompose them into resolved turbulence, subgrid-scale turbulence and mean advective.” What about the rest of the retrieved budget terms, such as “physical parameterizations and numerical diffusion and damping”? To clarify, it may help to move the sentence “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” (L261-263) here and mention that for general applications, other grid-resolved parameterized physics terms are possible and categorized as additional budget components. We moved the sentence as suggested and included the suggested additional sentence in L268-270. Figure 1: Panels a-c show the total turbulence (trb), which is the sum of the resolved and subgrid-scale components. I’m interested in their individual contributions. E.g., What is the relative magnitude between these two components? Do they have similar spatial distributions with the same signs, or do they offset each other? Considering that one of the distinctive features of this tool from previous WRF budget retrieval works is the decomposition into mean and turbulent components, I suggest strengthening the relevant discussion. The subgrid-scale component is mainly relevant very close to the surface and has a similar magnitude but opposite sign as the resolved turbulence component. We added a profile plot showing this decomposition (new figure 2) and some explanation in L294-300. L283: “...in the alternative form of the equation (Eq. 8), the correction term for the time derivative is almost negligible” I’m lost. This is not shown or can be inferred by figure 2? Yes, this is not shown (we added “not shown” in L316). The plot for Eq. 8 looks the same as Figure 3a. But the correction terms on the LHS of Eq. 8 and 11 are conceptually different (L134-136) since the prognostic variable is different in both equations (ρθ vs ρzηθ). L290-295: Have you checked if the sum of all the budget components in this alternative analysis is in close balance with the sum of your Eq. 11(13) in Fig. 1? Fig. 1 does not show the total tendency, only the turbulent and mean advective components, not the sum of both. The “total” in the plot refers to the sum of horizontal and vertical components. We added this clarification in the caption. The total tendency is shown in Fig. 3a. Close to the surface on the ridge the tendency is about 0.03 x 10 so on the order of 10 K s as written in L329 and visible in Fig. 3a. Technical corrections: L113, L124, ...: Replace “equation” with “Eq.” here. Please check the rest of the manuscript for consistency. Thanks, we fixed this at several locations. L134: “energy” --> “potential temperature” OK L154-155: “Since... to derive.” This sentence is confusing. Suggest changing to “Although the momentum variables are staggered differently from the thermodynamic variables, their discretized equations can be derived analogously...”

OK, we added these references in the introduction in L28-29. L35-38: I'm not sure how 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"? Can you please elaborate? We added a sentence there: "For instance, if the fluxes on two opposing sides of a grid box are equal, the flux form yields zero net tendency, while this is not necessarily the case for the advective form." (L39-40) Lee et al. 2004 looked at the budget of a larger control volume and noted that this is due to "internal processes" in the control volume. We don't want to explain this further in the paper, so we dropped the unclear mentioning of the "internal processes".
Chen et al. (2020) noted that a closed budget for w is challenging, partially because w equation is implicit, coupled with the geopotential tendency equation, and partially due to the variable's inherent rapid variation on small scales. How does WRFlux overcome this issue for the w budget? It is also surprising to see that WRFlux performs the best for w compared to all the other variables and that the accuracy (in terms of NRMS and r99; Table 1) is much higher for the momentum than for the thermodynamic variables (θ and q v ) in your case. Do you have any idea why? Is this result case-dependent or relevant to the physical design of WRFlux? Chen et al. found that closing the wbudget offline is difficult or impossible due to the mentioned issues, but the online budget calculation works fine for them. They also achieved a higher degree of budget closure for w than for horizontal momentum. Our budget calculations are also done online considering all budget components for w, including the important acoustic-step contributions.
In contrast to Chen et al. (2020), our budget calculations are time-averaged. To check the effect of the time averaging on the budget closure, we did two small test simulations with and without time averaging. Without time averaging, the NRMSE scores in the terrain-following coordinate system become considerably better for q v and w and slightly worse for θ, u, and v.
We did not investigate why the budget closure in the shown simulations is higher for momentum than for θ and q v . However, in some of our test simulations, it was the other way around, so it is definitely case-dependent (background wind speed, moisture content, physical parameterizations used, etc.).
It would help readers follow Section 2 more easily if the connections between subsections are made more explicit and some terminologies are explained more precisely. -For example, I find it confusing that "Eq. 11 is equivalent to Eq. 8" (L124-128), but "Eq. 8 cannot be closed numerically in contrast to Eq.11" (L128). This confusion may be related to how you define the term "close" here. Suppose you are referring to "budget closure". In that case, my understanding is that it depends on the accuracy of the conducted budget calculation, not the format of the chosen equation to which the budget analysis is applied. Do you mean that because Eq. 11 is numerically consistent with the equation used in the model, the left-hand-side terms of Eq. 11 represent more accurately the actual tendency simulated by the model? And therefore, a budget analysis using Eq. 11 can achieve a higher degree of budget closure? Please clarify.
Eq. 11 is similar to WRF's governing equations because the coordinate metric z η appears within the derivatives. In the manuscript, we showed how the budget can be closed with Eq. 11. When using the product rule to transform Eq. 11 to Eq. 8, two terms cancel out analytically, but not numerically. This is because one of the canceling terms comes from a horizontal derivative and the other from a vertical derivative. Therefore, the right-hand side of Eq. 16 is not identical to the right-hand side of Eq. 14. Thus, we argue that Eq. 16 cannot be closed numerically in WRF. The left-hand sides of Eq. 14 and 16, instead, are almost identical in our simulations.
To clarify this point, we added Eq. 12 and a bit of explanation in L167-170.
-Also, I suggest reminding readers at the beginning of Section 2.4 that Eq. 11 (i.e., Eq. 13 after discretization) will be used for the precise budget calculation. It would also be helpful to specify which terms in Eq. 11 can be decomposed into the mean advective and turbulent components (all righthand-side terms excluding S or all right-hand-side terms excluding only some partial S, i.e., the subgrid physics?) OK. We added a few introductory sentences in Section 2.4 with the suggested content (L184-189).
L113-L114: "Numerical consistent" --> "Numerical consistency". Also, placing this sentence ("Numerically… after discretization.") here seems odd. You are neither showing the discretization nor discussing the budget closure in the following lines yet.
That's right. Nevertheless, in our opinion, this sentence is important for the logical train of thought. We added a reference to the discretization section in L116.
L130: It is still unclear to me why recalculating the vertical velocity is desired. If Eqs. 6 and 9 are also used in WRF, doesn't it make more sense to use WRF's prognostic to be numerically consistent with the model? Unless this is relevant to the issue I mentioned in comment #4?
Recalculating w is necessary to obtain a closed budget in the Cartesian coordinate system. Originally, there was a significant difference in the calculation of z η ω in WRF's implementation of Eq. 6 and 9. We reduced this difference strongly by changing the discretization of Eq. 9 in WRF (L216-221, https://github.com/wrfmodel/WRF/pull/1338). Nevertheless, there is still a small difference remaining: In Eq. 6, z η is included in the coupled prognostic variable µ d ω, whereas in Eq. 9 it is calculated directly from the geopotential. To overcome this difference and obtain a closed budget, we recalculate w in a consistent way and use it in the vertical advection term. Whether this is appropriate is definitely debatable. But the effect on w itself is hardly noticeable. Only when adding up the large and counteracting budget components, the difference becomes noticeable. With the prognostic w, the budget closure (NRMSE) is about one order of magnitude worse (L356-359). We suppose the reviewer wanted to ask whether this is relevant to the issue mentioned in comment #5, not #4. These issues are not connected. The recalculation of w mentioned is only relevant for the transformation to the Cartesian coordinate system. It is not used as the prognostic variable in the conservation equation for w. We clarified this in L217-219.

L191-194: Does this part follow exactly Chen et al. (2020)?
The retrieval of subgrid-scale tendencies is different since we only extract the fluxes and then do the tendency calculation offline (as for the resolved fluxes). The retrieval of other forcing terms (physics, pressure gradient, damping, diffusion, Coriolis) was implemented very similarly to Chen et al. but independently from them. Over a user-specified time window (we added this in L213). In this study, we use 30 min (L272) as is often used for Reynolds averaging (compromise between large sample for statistics and stationarity). We added a sentence with a reference on this in L272-274.

L203-204: "Since WRF does not actually solve the continuity equation..." This may be misleading. WRF does solve the continuity equation, just not in the same format as expressed in your Eq. (12), i.e., in terms of the 3D variable ρ but column dry air mass per unit area (2D µ d ). Note that µ d is indeed advanced with time in the WRF dynamical solver. I suggest changing the sentence to something like "Since WRF does not directly solve for ρ but the integrated column dry air mass…"
Actually, WRF's mass continuity equation in terms of µ d (Eq. 2.12 in the technical notes version 4, http://www2.mmm.ucar.edu/wrf/users/docs/technote/v4_technote.pdf) does correspond to the rightmost term in our Eq. 13, except for a factor -g. In fact, µ d and ρ d are related by a diagnostic equation: µ d = -g ρ d z η . Ensuring that the last term of our Eq. 13 is exactly zero is difficult because the mass continuity equation in WRF is indeed solved but not integrated explicitly in the acoustic time steps, but rather µ d is diagnosed from its definition (Eq. 3.20 in the technical notes) after vertically integrating the mass divergence (Eq. 3.22). We acknowledge that our original formulation was misleading, and we replaced it with the information given in the preceding sentence (L223-225).
L244: "The setup leads to … a very dry atmosphere, therefore moist processes are neglected." Not sure if I overlook something but the initial setup of the moisture field is not mentioned?
Yes. We added the specification of the moisture profile (constant RH=40%) in L263.
L246: "We calculate full θ-tendencies and decompose them into resolved turbulence, subgrid-scale turbulence and mean advective." What about the rest of the retrieved budget terms, such as "physical parameterizations and numerical diffusion and damping"? To clarify, it may help to move the sentence "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" (L261-263) here and mention that for general applications, other grid-resolved parameterized physics terms are possible and categorized as additional budget components.
We moved the sentence as suggested and included the suggested additional sentence in L268-270.

Figure 1: Panels a-c show the total turbulence (trb), which is the sum of the resolved and subgrid-scale components. I'm interested in their individual
contributions. E.g., What is the relative magnitude between these two components? Do they have similar spatial distributions with the same signs, or do they offset each other? Considering that one of the distinctive features of this tool from previous WRF budget retrieval works is the decomposition into mean and turbulent components, I suggest strengthening the relevant discussion.
The subgrid-scale component is mainly relevant very close to the surface and has a similar magnitude but opposite sign as the resolved turbulence component. We added a profile plot showing this decomposition (new figure 2) and some explanation in L294-300. Yes, this is not shown (we added "not shown" in L316). The plot for Eq. 8 looks the same as Figure 3a. But the correction terms on the LHS of Eq. 8 and 11 are conceptually different (L134-136) since the prognostic variable is different in both equations (ρθ vs ρz η θ).

L290-295:
Have you checked if the sum of all the budget components in this alternative analysis is in close balance with the sum of your Eq. 11(13) in Fig.  1? Fig. 1 does not show the total tendency, only the turbulent and mean advective components, not the sum of both. The "total" in the plot refers to the sum of horizontal and vertical components. We added this clarification in the caption.
The total tendency is shown in Fig. 3a. Close to the surface on the ridge the tendency is about 0.03 x 10 -3 so on the order of 10 -5 K s -1 as written in L329 and visible in Fig. 3a.
Thanks, we fixed this at several locations.
Yes. However, since we included an additional equation (Eq. 12) in the revised manuscript, the labels are correct now. Yes Table 1: Same issue as above.