Evaluation of a forest parameterization to improve boundary layer flow simulations over complex terrain

We thank the reviewers for their careful reading of our manuscript and constructive comments. Their comments helped us to reline our manuscript by clarifying several sections and by improving the presentation of the results. Pointby -point responses to individual questions and remarks for each reviewer are listed below. Note: The referee’s comments are reproduced in extenso in bold characters. The line numbers mentioned in our answers are those of the revised version of the manuscript.

1. Canopy aware LES (also sometimes called canopy resolving LES) has become reasonably mainstream in the last 10 years, with many community LES (i.e. WRF, DALES, PALM) having developed capabilities to represent canopy drag forces. It have been shown that such treatment, is capable of improved representation of turbulent statistics within and above the canopy (e.g. TKE, sigma w, sweeps/ ejections). It has been shown that a sufficient number of within canopy model layers are needed to do so (order 10 layers) and that results are best when model resolutions are isotropic (i.e dx = dy = dz).
Given the current setup of 2-3 model layers and dx/dy of 200 m in d03 (which would produce a grid aspect ratio of roghly 10:200), I am wondering to what extent the canopy drag formulation would be capable of improving representation of turbulent structures above the forest.
This brings me to my main question, which I could not find addressed in the manuscript: To what extent is the comparison between the forest parameterization and non-forest a true apples to apples comparison? The non-forest LES uses z0 from the land-use dataset (I assume this is the built in WRF one, that may not be very accurate at such high resolutions), while the forest parameterization uses a direct sink for momentum. A better comparison may be to apply a similar momentum sink to the first model layer only, or to first calibrate the model's z0 to observed values.
We appreciate the referees comments and we would like to comment as follows: • The ultimate goal of the exercise is to improve the "out of the box" model performance in a highresolution mesoscale application by including a simple forest parametrization. Here, we would like to clarify that although we are using LES techniques to deal with turbulence closure, we do not consider the simulations as truly LES due to the horizontal resolution used (even for the D04 domain). Since in our simulations (at best) we used a horizontal resolution of 40 m, the resolved motions in the grid do not belong to the Inertial subrange of the flow. To resolve the more energetic turbulent eddies, specially for the low level jet episodes observed during a stable nighttime boundary layer, a horizontal resolution of about O 5 m is then required (Cuxart, 2015).
Therefore, we believe that our simulations are closer to high-resolution mesoscale simulations with LES techniques for turbulence closure. This has been also clarified in the new version of the manuscript in line 61.
The aim of the present work is not to evaluate the performance of the model in terms of the reproduction of turbulent processes (we are not completely resolving turbulence in the simulations), instead, we are more interested in evaluating the response of the mean flow to the introduction of the forest parametrization.
The comparison with the default Zo is carried out because we need to chose a "control case" to assess the performance of our development. We agree with the referee that a tuning of the Zo for the simulated domain may improve the results of the non-forest simulation. Indeed, Wagner et al. (2019) observed that surface roughness lengths from CORINE land use data presented considerably smaller roughness lengths than they should be. In the double ridge, CORINE roughness lengths are of the order of 0.1 m, while in reality the hills were partially covered by eucalyptus trees with heights of about 20 to 25 m, which should be represented by roughness lengths of the order of 1 to 2 m. However, the ultimate goal of the exercise is to improve the model performance in such high resolution mesoscale simulation using the available tools/datasets.
One of the interesting ideas of the implemented forest parametrization is that at the end it only requires the leaf area index (LAI) as input. The LAI is already part of CORINE data base, which results in a rather inexpensive implementation process, especially when the evidenced improvement in the model performance is taken into account.
2. LLJ cases: The model evaluation against the LIDAR data is done as instantanious cross section plots. While I can see that there the forest parameterization exhibits more dampening of lee-waves and has somewhat lower near surface wind speeds, it is generally hard to make out finer differences in the plots. I would suggest that the authors think about how to better visualize these model results. For example, one might want to focus on the d04 domain and would also want to zoom in closer to the surface. Alternatively, difference plots would also be approciate. Lastly, I am wondering, given the 5 minute output interval, why comparisons are done for a single snapshot in time. Given the amount of output data, it would be interesting to see whehter there is a better quantiative comparison possible. Furthermore, Wagner et al 2019a presents a model evaluation paper with respect to d03 domain in non-forest mode. Compared to the model evaluation for d03 in the current manuscript, this evaluation appears to be much more detailed and I am wondering why this is the case.
• The idea behind the use of the cross section view is to fully observe the structure of the lee waves formed around the Perdigao double ridge. The comparison using fixed time Lidar observations has its origin in the fact that we are interested in assessing the performance of the model to reproduce the lee waves generated by the topography. The selection of the depicted times, were base on the selection of periods with quasi-steady structures (see Wagner et al. (2019)).
• We decreased the vertical limits of the graphs to focus more on the lee wave structure close to the topography to provide an additional quantitative comparison of the lee wave structure. Vertical profiles of the wind speed for both simulations and observations have been included in the new figure 7. The text has been rephrased accordingly in section 3.2.
• The D04 simulations are compared for longer periods in figures 9, 10 and 11. It worth noting that the Hovmöller diagrams (figures 10 and 12 in the new version of the manuscript) aim to test the D04 performance.
• We build our work on Wagner et al. (2019). There, the longer simulations without the forest parametrization are already validated with observations. The focus in the present work is to assess the "general" improvement of the simulations. For d03, the evaluation a "typical" day during the IOP is assessed. Therefore our focus on the average diurnal cycle through the IOP.
3. I am not familiar with the data from the field campaign that is available to the authors, but I am wondering whether there is a missed opportunity in not having any kind of turbulent quantities that the LES is evaluated against (e.g. from the towers or the LIDARs). Usually, the fact that LES are capable of partially resolving turbulence and to provide realistic profiles of turbulent statistics that can be compared against observations is seen as one of the crucial adavantages of the LES. In this manuscript, there is no consideration of turbulent quantities, which could be used to evaluate the actual LES.
• As mentioned above our focus is not on the turbulent structure. For the nighttime periods of interest, even the 40 m simulation corresponds rather to a high-resolution mesoscale simulation (Cuxart, 2015), and not a "true" LES which resolves motions down to the inertial-scale turbulence. For these stable boundary layers simulations higher horizontal resolutions are then required (O 5 m). Nevertheless, we do refer to the simulations as LES, because a LES-type turbulence closure is used.
• In addition to the fact that the focus of the present work is not the turbulence (as mentioned above), we would like to highlight that in the thermally driven low level jet episodes studied in the present manuscript (D04 simulations), the flow is not really turbulent. The potential temperature contours in the cross-section plots (figures 6, 8 and 9 of the manuscript's new version) suggest a rather laminar flow. As evidenced in the Hovmöller diagram for wind speed in figure 10 (new version of the manuscript), the flow becomes turbulent during the daytime between 09h00 and 18h00 UTC, period which is not really studied in the present work and it is left for future work.
4. Relationship to retracted Wager et al 2019b: The authors mention in the submission that this work builds on the retracted Wager et al 2019 b work. A cursory look at the discussion of this manuscript shows that two reviewers had issues with the WRF setup in the sense that they would have liked to see either a nudging or periodic restart for the long-term simulation. I would like to know whether these comments have been heeded in the present work. Similarly, reviewers of the W19b manuscript would have liked to see the full namelist of the WRF for added transparency/ reproducibility. I feel that something like this would be particularly valuable given the objectives of GMD. Lastly, since W19b was retracted, it should not be referenced in this manuscript.
• The comments done to ? have not been considered in the present work. We don't think that periodic restarts or nudging would have a significant positive effect on our simulations. The flow at the Perdigao site is dominated by mesoscale flows originating from the interaction of the synoptic flow with the nearby mountains and thermally forced flows. All of which is not well represented in the coarse resolution driving model.

Specific comments
1. P1L14: "The grids are now fine enough to fully resolve the atmosphere with techniques such as large eddy simulation (LES)" I take issue with this statement, as it is incorrect. LES partially resolves turbulence, but this is nowhere near fully resolving (especially given the model resolutions of 200 and 40 m in this paper). This needs to be rephrased.
It has been corrected in the updated version of the manuscript. The introduction has been modified according to the referee's comment.
Response to Reviewer # 2: To my opinion, the paper is worth publishing in GMD, although some modifications would be needed. The paper is somewhat uncompensated, as the 49-day run is very succintly described while the changes are relevant, while the high resolution runs are discussed extensively but more superficially. Also the Hovmöller plots are hard to follow as the details explained in the text are difficult to be seen. Sensitivity tests would have been interesting, varying the height of the canopy, and some discussion on the effects of the changes in the temperature profiles (not shown) would have been very welcome.
We thank the referee for his/her kind comments. We briefly described the 49-day simulations because an extensive analysis of the simulation has been carried out in Wagner et al. (2019) and only focus on the difference between the "out of the box" simulation and the simulation using the forest parametrization, which we do believe is missing from the later article.
Regarding the analysis of the high-resolution simulations, several modifications have been introduced to the manuscript (including new figures). The latter will be described in the replay to the referees' specific comments (see below). Finally, we discussed the effects of the implemented modifications in the temperature profiles in this replay. However, we concluded that it is not worth including in the manuscript (see below).
Referee's specific remarks 1. A discussion on the choice of the vertical resolution, especially the location of the lowest point at 10 m above the surface. Why the lowest level is at this height? Couldn't the first level be located closer to the surface as other studies in complex terrain do? Is there a computationally-related reason, such as numerical instabilities on steep slopes, or is it because of some choice related to the applicability of the similarity theory?
Previous work in terrain more complex than that presented in this work have found that an increase of the vertical resolution does not necessarily help the model to improve the results (Quimbayo The manuscript has been modified as follows: (Line 117) Previous work has shown that a first model level at 20 m is adequate for this kind of high resolution mesoscale simulations (Quimbayo et al., 2020, Umek et al., 2021.

2.
In fact no explicit explanation is given about the use of the similarity theory near the surface, and what are the equations employed.
The scheme used (see namelist in the archive alongside the manuscript) is the Eta Similarity (Monin-Obukhov-Janjic) scheme (sf sfclay physics = 2, in the namelist). The surface layer is formulated following the well known similarity theory by Monin and Obukhov (1954). It provides the lower boundary conditions for the Level 2.5 model (Janjic 1996). The Beljaars (1994) correction is applied in order to avoid singularities in the case of free convection and vanishing wind speed (and consequently u).
The manuscript has been modified as follows: (Line 121) The Eta Similarity (Monin-Obukhov-Janjic) scheme is used to couple the atmosphere and the surface. The scheme provides the lower boundary conditions for the Level 2.5 model (Janjic, 1996). The Beljaars (1995) correction is applied in order to avoid singularities in the case of free convection and vanishing wind speed (and consequently u).
3. Justify the choice of the Lalic-Mihailovic shape for the leaf-area density. Is this the only available possibility or are there others and this one has shown to be the best option?
The manuscript has been modified as follows: (Line:170) The profile proposed by Lalic and Mihailovic (2004) is chosen due to its flexibility to represent a broad range of forest canopies by only choosing the appropriated value for the z m parameter. Kolic (1978) forest classification, based on z m and h parameters, let us to divide the forest canopies into three groups: (1) zm = 0.2 h (oak and silver birch), (2) 0.2 h< z m <0.4 h (common maple), and (3) z m = 0.4 h (pine). The parametrization is adaptable for several cases. Following Kolic (1978) and Lalic and Mihailovic (2004), the empirical relation for the LAD can be applied for the broad range of forest canopies. The site at Perdigao presents an irregular vegetation coverage, made of no or low-height vegetation and patches of eucalyptus and pine trees (Vasiljevic et al., 2017). The section of z m , consisted with Mohr et al. (2014), is based on the idea of a forest canopy combination between pines and taller eucalyptus threes.

4.
A more elaborated discussion of what it means to use a canopy-drag compared to surface roughness would be appreciated. In the latter case it is costumarily assumed that a logarithmic profile is imposed between the lowest point and z0 and that the similarity expressions hold (even if they were originally derived for flat terrain). Instead introducing the LAD profile and the corresponding drag breaks this conceptual model. I think the paper would benefit of a reflection on how the surface layer is different in both approaches, perhaps showing comparison profiles. Also, is similarity theory still used in the lowest level below the canopy profile?
The manuscript has been modified as follows: (Line 139) The total drag can be modelled by adding the contribution of the orographic pressure (form) drag, the canopy drag and the (surface) frictional drag. In the standard WRF-LES model, three-dimensional roughness elements like trees are not explicitly considered and are only characterized by the roughness length Z0 (directly obtained from the landuse data set), implying that mainly the frictional drag is considered. The explicit treatment of forest drag in numerical models is of special importance for the realistic development of wind profiles, including inflexion points over forested areas enhancing flow separation.
In the present work, the forest parameterization proposed by Shaw and Schumann (1992) is implemented in the WRF model to study its impact on boundary layer flows over forested and complex terrain. The additional forest drag term F i (direct sink term in the momentum equation), which acts on the lowermost model levels is defined as where |V | is the magnitude of the three dimensional wind vector, u i is one of the three wind components, C d = 0.15 is a constant drag coefficient, and LAD is the leaf area density profile characterising the trees. The LAD depends on the tree type and the height of the trees meaning that the strength of F i varies as one move in the vertical inside the canopy. The tree type is defined by means of the leaf area index (LAI). Similarity expressions hold for the first model level, as the parametrization only introduce a momentum sink, it still needs to couple the atmosphere and the surface. It is worth noting that the flow momentum is mostly absorbed by the canopy, meaning that the details of the lower boundary condition becomes less important there (Ross and Vosper, 2005). The LAD-profile is computed according to Lalic and Mihailovic (2004) at all grid cells where trees are present as (Line 50) Different authors found a significant improvement in the performance of numerical simulations when using a canopy drag scheme. Ross and Vosper (2005) found that simulations using a forest canopy (within moderate complex terrain) performs reasonably better than simulations using the roughness-length parameterization when tested to experimental wind tunnel data. The roughness-length parameterization generally leads to significant underprediction in the pressure drag compared to the explicit representation of canopy using the canopy drag scheme. Finnigan and Belcher (2020) highlight the fact that mechanisms of "separated sheltering" in flow over low hills, which is increased with a canopy and is the dominant mechanism for drag over low hills, may become double the topographic drag if a roughness-length scheme is replaced by a forest canopy scheme.

5.
The statistics with the modified run are much better than the standard one for the along-valley wind, while the cross-valley is not significally improved. The latter is attributed by the authors to the use of a too high canopy and to the heterogeneous distribution of forested areas. This aspect is superficially commented and could have been complemented with some sensitivity runs of short duration (only two levels with canopy during one day for instance), as it is an important point in the results.
As noted in the manuscript, the problem may come from the representation of the forest canopy (tree height, canopy extension, tree distribution, etc.) at the north of the double ridge. The idea of a set of sensitivity test is appreciated, but the problem is that we need to test different variables as canopy deepness, tree height and canopy extension, which translates into several simulations. Each of these simulations needs to be long enough to complete several diurnal cycles to be comparable with the results reported in the manuscript ( figure 4 and figure 5). The inversion in time and resources for such a detailed investigation would be considerable and beyond the scope of the present work. We leave this as an open question for future research.
6. I find the plots of the short-term simulations far from clear when differences between using forest or not are commented. For instance, to my eye it is difficult to find the "southwest flow observed near the ground on the lee of the topography" (lines 195-196). Perhaps a circle on the figure would make the task easier for the reader. In general fig 6 is not good to follow the explanations given in the text, maybe vertical profiles at the points of interest comparing the two runs would be more informative.
As it is now, I find myself believing what the authors say. I suggest that the whole discussion concerning figures 6 to 8 is remade with more clear graphic evidence of the results that the authors indicate.
Vertical profiles of the wind speed for simulations and observations has been added at selected places to help the reader to understand figures 6, 7 and 8 (see new figure 7). We appreciate the referees comments and we believe the section is now more complete. Several modification has been introduced to section 3.2 of the manuscript. 7. when comenting figure 9 it is hard to see almost anything that is commented in text. I would suggest to show perhaps only the diferences between simulation and observation so that the relevant issues jump to the reader easily. Some graphical information on temperature -and the corresponding discussion-is missing too, not only here but all along the text. 2017-05-08T05:08 Figure I: Instantaneous potential temperature profiles from radio soundings (black lines), WRF simulations using the forest parametrization (green lines) and without the forest parametrization (blue lines). Colour lines represent the location where the radio sounding were launch (39.713680N, 7.736365W) in the D04 domain.
• We believe that the Hovmöller diagrams, as they are, can be an important tool to get a general view of the evolution in time of the wind in the area (creation and vanishing of the LLJ). We have clarified the text.
We believe a discussion on the effect of the forest canopy on temperature is unnecessary. As the forest canopy parametrization is only a sink term in the momentum equation, we do not expect important differences between the simulations. The latter can be observed in figure I, where instantaneous vertical profiles of potential temperature are plotted from radio soundings, and both WRF simulations. The results show that the model struggles to simulate the near surface thermal structure, and the implementation of the forest canopy does not help (neither makes it worse) to improve the results.
As the manuscript is aimed to show the influence of the forest parametrization on the results, we believe that a discussion of the thermal structure of the atmosphere does not add any beneficial information to the manuscript.
A note has been added on the small effect of forest parametrization on the thermal structure of the atmosphere as follows: The analysis and evaluation of the model focus on the performance of the model when simulating the wind field only. The simulations showed only small differences in the thermal structure of the atmosphere (not shown), and further investigation seems unnecessary.
8. Instead I find Figure 10 a good summary representation of what is going on, it would be nicer if it was given together with average profiles for wind and temperature in order to see in what parts of the column the effect of the changes is more evident for wind and temperature.
• We appreciate the referee's comment and we have added average profiles for wind speed to former figure 10 (figure 11 in the new version).