Carlos Román-Cascón
email: carlosromancascon@ucm.es

Marie Lothon

Fabienne Lohou

Oscar Hartogensis

Jordi Vila-Guerau De Arellano

David Pino

Carlos Yagüe

Eric R Pardyjak

Carlos Román

Surface representation impacts on turbulent heat fluxes in WRF (v.4.1.3)

The water and energy transfers at the interface between the Earth's surface and the atmosphere should be correctly simulated in numerical weather and climate models. This implies the need for a realistic and accurate representation of land cover (LC), including appropriate parameters for each vegetation type. In some cases, the lack of information and crude representation of the surface leads to errors in the simulation of soil and atmospheric variables. This work investigates the ability of the Weather Research and Forecasting (WRF) model to simulate surface heat fluxes in a heterogeneous area of southern France, using several possibilities for the surface representation. In the control experiments, we used the default LC database in WRF, which differed significantly from the actual LC. In addition, sub-grid variability was not taken into account since the model uses, by default, only the surface information from the dominant LC category in each pixel (dominant approach). To improve this surface simplification, we designed three new interconnected numerical experiments with three widely-used land-surface models (LSMs) in WRF. The first one consisted of using a more realistic and higher-resolution LC dataset over the area of analysis. The second experiment aimed at investigating the effect of using a mosaic approach, where 30-m sub-grid surface information was used to calculate the final grid fluxes, based on weighted averages from values obtained for each LC category. Finally, in the third experiment, we increased the model stomatal conductance for conifer forests, due to the large fluxes errors associated with this vegetation type in some LSMs. The simulations were evaluated with gridded area-averaged fluxes calculated from five tower measurements obtained during the Boundary Layer Late Afternoon and Sunset Turbulence (BLLAST) field campaign. The results from the experiments differed depending on the LSM and displayed a high dependency of the simulated fluxes on the specific LC definition within the grid cell, an effect that was enhanced with the dominant approach. The simulation of the fluxes improved using the more realistic LC dataset except for the LSMs that included extreme surface parameters for the coniferous forest. The mosaic approach produced fluxes more similar to reality and served to improve, especially, the latent heat flux simulation of each grid cell. Therefore, our findings stress the need to include an accurate surface representation in the model, including soil and vegetation sub-grid information with updated surface parameters for some vegetation types, as well as seasonal and man-made changes. This will improve the modelled heat fluxes and ultimately yield more realistic atmospheric processes in the model.

Copyright statement.

Introduction

. These changes can be extremely important because they modify the natural cycles of energy (Seneviratne et al.,

and higher resolution LC dataset from the Centre d'Etudes Spatiales de la Biosphère research laboratory (CESBIO, Inglada et al. (2017)). The results showed a high dependence of the flux on the specific LC categories, which motivated a second experiment including the sub-grid information of the surface, the so-called mosaic approach (e.g., [START_REF] Li | Development and evaluation of a mosaic approach in the WRF-Noah framework[END_REF]. Finally, an additional experiment was carried out due to the extreme biases found in the first two experiments over those pixels mostly covered by conifer trees in Noah-MP, with the aim of diminishing the biases by modifying some parameters associated with the transpiration processes of this LC category.

For the evaluation of the simulations, we took advantage of the large number of instruments deployed during the Boundary-Layer Late Afternoon and Sunset Turbulence (BLLAST) field campaign, carried out in 2011 in southern France (Lothon et al., 2014). The spatial density of eddy-covariance (EC) towers over different vegetation types facilitated the calculation of gridded area-averaged fluxes (AAF), as done in a similar way for the LITFASS experiment [START_REF] Beyrich | Areaaveraged surface fluxes over the LITFASS region based on eddy-covariance measurements[END_REF], where they obtained a good agreement with the fluxes measured from scintillometry. In the present work, the AAF were used to evaluate the results from the WRF model coupled with the three LSMs under the different conditions set in the experiments, and to analyse the results based on the different LC types.

The article is organised as follows: Section 2 provides information of the measurements taken during the BLLAST field campaign and explains how the area-averaged fluxes were calculated. Detailed information about the model configuration and the different experiments are also included in this section. In Section 3, we quantified the results from the different modelling experiments, including scientific discussion about them. Finally, a short summary and the main conclusions are provided in Section 4.

Evaluation data, WRF model and experiments design 2.1 Observational data for model evaluation

The surface turbulent heat fluxes simulated by the WRF model were evaluated during a period of the BLLAST field campaign (Lothon et al., 2014). This campaign took place from 14 June to 8 July 2011 on the Plateau of Lannemezan (southern France).

Its main objective was to better understand the turbulence decay observed during the afternoon transition, and the extensive instrumentation deployment included several surface-energy-balance towers installed over different surfaces, which were representative of the vegetation within the explored area: prairies, forests, wheat, corn and moor. The analysed period comprised from 09:00 UTC to 15:00 UTC 19 June 2011, corresponding to part of IOP (intensive observation period) 2 of the campaign. This IOP was characterised by fair weather, no clouds, and a typical development of the boundary layer up to 800-1000 m above ground level (agl). The IOP is representative of the general conditions of the rest of the IOPs of the campaign (the general meteorological conditions of the campaign can be found in Lothon et al. (2014)).

SH and Le were calculated uniformly, using the eddy-covariance (EC) method with the specifications indicated in [START_REF] De Coster | BLLAST-uniform processing of Eddy-Covariance data[END_REF] over five different LC types: grass, wheat, corn, moor and forest (conifers). The position above ground level (agl) of the instruments was set according to the vegetation height, which implied an homogeneous footprint for the five towers: 2 m for the moor, grass and wheat sites, 4 m for the corn site and 31 m for the forest site (approx. 6 m above the trees). These measurements showed that the SH was more sensitive to LC type (Fig. 1a andc) than Le (Fig. 1b andd). This was probably linked to the different surface properties of the vegetation types (albedo, thermal inertia, emissivity), with more influence on the surface temperatures. The highest SH was measured over the forest site, followed by the wheat site, with peaks of around 400 and 270 Wm -2 respectively. The measurements over the grass, moor and corn sites showed lower values with maximum SH of around 150 Wm -2 during the central hours of the day. However, Le exhibited values that were similar for all the vegetation types, with midday maxima around 300 Wm -2 . Nevertheless, the lowest Le was measured over the wheat field (this crop started to dry during the experiment), while the largest values were observed at the grass, forest and moor sites. Some differences were observed in the diurnal evolution of Le: the morning Le rise over the corn/wheat were delayed with respect to the measurements over other LCs, maybe linked to a delay in the transpiration processes associated with these crops (a full investigation of the reasons of it is out of the scope of this study). In any case, it should be highlighted that these fluxes values were representative of those of the different IOPs of the BLLAST field campaign (Lothon et al., 2014;[START_REF] Couvreux | Boundary-layer turbulent processes and mesoscale variability represented by numerical weather prediction models during the BLLAST campaign[END_REF], with no-water limitation due to the regular rain events observed in the area during the weeks before, which are typical conditions in these dates in the area.

Area-averaged fluxes (AAF) calculation

The EC measurements from towers were used for the calculation of AAF with resolution of 1 km in a region of 19x19 km around the central site of the field campaign: 43 • 07 27.15 N, 0 • 21 45.33 E; 600 m above sea level (asl). This approach can be considered as a mosaic observational method for flux computation. It was based on the multiplication of the values of the fluxes measured over each vegetation type by their respective vegetation-cover fraction within each 1x1 km pixel (see Fig. 2a,b). Finally, the contribution from each vegetation type was summed to provide the final values in each grid cell (area-averaged fluxes). High-resolution and realistic LC data were needed to construct the AAF, obtained from the 30-m LC database prepared by CESBIO, based on Landsat-5 data [START_REF] Inglada | Operational high resolution land cover map production at the country scale using satellite image time series[END_REF] (Fig. 2b). Hence, the main assumption of the AAF calculation was to consider the same SH and Le for sites with the same LC type at different locations, implying also the following simplifications (further discussion is included after the list):

1. There were no flux differences depending on the soil type.

2. There was no difference in activity for the same vegetation type in the area.

3. There was no horizontal variability of the soil moisture (SM).

4. The radiation and the wind forcing was the same for all the pixels, without altitude influences. 5. All the forests were considered coniferous (evergreen needle-leaf forest (ENF)).

6. Fluxes over urban and bare-soil surfaces were estimated using a simple model.

150

(1) Two soil types were the dominant ones in the evaluation area, based on the 1x1 km pixels from the soil database of WRF (United States Department of Agriculture, USDA): clay loam (81% of the area) and loam (19% of the area). Since the dominant soil type in the pixels where the measurements were taken was clay loam, those pixels with loam were not included in the evaluation of the model. Therefore, the limitation based on the soil type was ultimately not a problem. This area is shown in Fig. 2d with a black rectangle.

(2) The area analysed is relatively small to expect significant differences between the natural vegetation belonging to the same category. Although the limitation based on the possible different vegetation state can be more important in crop areas, seeding dates are normally coincident between fields in the area.

(3) Regarding the effect of possible SM differences within the area, two aspects should be noted. On the one hand, the potential SM horizontal variability due to possible inhomogeneous precedent rainfall over the area was not taken into account.

However, the SM input in the model is provided with a coarse resolution of 1 • and does not show any small-scale details.

This is a well-known limitation of mesoscale models which is sometimes addressed through the assimilation of SM data from satellites or with previous long simulations that serve to spin-up the surface in order to obtain the appropriate SM initial values [START_REF] De Rosnay | A simplified Extended Kalman Filter for the global operational soil moisture analysis at ECMWF[END_REF][START_REF] Angevine | Land surface spinup for episodic modeling[END_REF][START_REF] Santanello | Impact of soil moisture assimilation on land surface model spinup and coupled land-atmosphere prediction[END_REF]. In our case, this limitation of the mesoscale modelling is an advantage because it allows to perform a fairer model-observation comparison since this limitation also exist in the AAF.

The impact of applying a spin-up period in the model is investigated later.

On the other hand, part of the SM horizontal variability is caused by the different properties of the vegetation patches, affecting runoff, infiltration, evapotranspiration, etc., and, finally, the SM content associated with each vegetation type. This LC effect on SM may be implicitly included in the fluxes that were measured over each individual vegetation type. In any case, we investigated the effect of including the SM horizontal heterogeneity with a higher resolution in the model initialization (not shown). This was done by using high-resolution SM satellite data from the Disaggregation based on Physical And Theoretical scale Change (DISPATCH) [START_REF] Merlin | Self-calibrated evaporation-based disaggregation of SMOS soil moisture: An evaluation study at 3 km and 100 m resolution in Catalunya, Spain[END_REF][START_REF] Molero | SMOS disaggregated soil moisture product at 1 km resolution: Processor overview and first validation results[END_REF] product at 1 km of resolution, but conserving the original range of SM variability within the area evaluated. The effect of including this SM spatial heterogeneity was minimum in the fluxes simulations (in part due to the conservation of the range of SM values).

(4) The limitation due to the radiation and wind forcing, which were assumed as equal in all the pixels, is expected to have a small impact, due to the fair weather, with no clouds, and light wind conditions on the 19 of June (Lothon et al., 2014).

However, some differences could exist between some pixels situated at different altitude over the area (most of the pixels were in the range of altitude between 400 and 650 m asl, with minimum and maximum altitude of 282 m and 696 m asl in the whole evaluated area).

(5) The limitation of considering all forests as coniferous was due to the fact that the campaign only included measurements over this type of trees, and not over broadleaf deciduous forests (the other forest category in the area). Therefore, the measurements taken over conifers were extrapolated to those areas covered by deciduous for the calculation of the AAF. This limitation was addressed by converting all the forest in the WRF model to conifer, which allowed a fairer model-observations comparison at the expenses of loosing the analysis over one LC category (deciduous forest).

(6) For the urban surfaces, SH and Le were calculated following a simple Penman-Monteith type model [START_REF] De Bruin | A simple parameterization of the surface fluxes of sensible and latent heat during daytime compared with the Penman-Monteith concept[END_REF] shown in Eq. ( 1) and Eq. ( 2), respectively:

SH = β • R n -G 1 + β , (1)
Le = R n -G 1 + β , (2)
where R n (Eq. 3) is the net radiation, defined as:

R n = (1 -α) • SW ↓ +LW ↓ -σ T 4 s , (3)
and G (Eq. 4) is the ground heat flux, considered a fixed fraction of R n :

G = G f rac • R n . (4)
In these equations, SW ↓ and LW ↓ are the measured downward short-wave and long-wave radiative fluxes, respectively.

The albedo (α), emissivity ( ), the fraction of energy used for ground heat flux (G f rac ), and the Bowen ratios (β) were considered constant, with values of 0.15, 0.92, 0.3 (stable)/0.5(unstable) and 5, respectively, obtained from [START_REF] Lemonsu | Modeling the surface energy balance of the core of an old Mediterranean city: Marseille[END_REF] and [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF] for urban surfaces. These simplifications could have led to some overestimation of urban effect on the total fluxes, since the urban surfaces in the area were surrounded by gardens, prairies and trees (so-called urban diffuse in the CESBIO LC dataset). However, the 30-m resolution LC dataset should be appropriate to deal with the urban/vegetation patches.

Despite the high-quality data and the efforts carried out to reduce uncertainties in the AAF, the evaluation of models with surface observations is a necessary task which is not straightforward. The observational measurements are almost always linked to uncertainties that can add limitations to the evaluation process [START_REF] Bou-Zeid | The Persistent Challenge of Surface Heterogeneity in Boundary-Layer Meteorology: A Review[END_REF]. In our case, the measurements were taken at heights that implied homogeneous footprints, but even in this case, other uncertainties are always present [START_REF] Mauder | Surface-Energy-Balance Closure over Land: A Review[END_REF], especially in the case of the measurements taken over vegetation with tall canopies, where the heat storage can have an important role in the surface energy balance.

Performance scores

Three performance scores were used along the article to evaluate the modelled fluxes, using as benchmark the AAF. The first one is the Bias (Eq. 5), subsequently named as Bias:

Bias = N i=1 [F lux sim(i) -F lux obs(i) ] N , (5)
where F lux sim(i) is the simulated SH or Le for each grid cell and each analysed hour, being N the total number of data used, and F lux obs(i) is the respective value computed from the AAF (observations). The average was performed for all the analysed hours (from 09:00 UTC to 15:00 UTC) and for all the grid cells covered by certain land-cover category within the analysed area (note how in Table 3 these values are averaged for the whole area without land-cover distinction). Hence, although the bias provides an indication of the systematic overestimation or underestimation of the modelled flux, the averaging can lead to undesired conclusions, since positive and negative values can compensate in some cases (especially if contrasting LC categories 215 are merged).

We also include the root-mean-square error (RMSE), being the root mean of the squared differences between model and observed fluxes at each analysed hour and for all the grid cells with certain land-cover category:

RM SE = N i=1 [F lux sim(i) -F lux obs(i) ] 2 N . (6)
Finally, the standard deviation of the error (Eq. 7) is also included, which provides an indication of the model random error, 220 i.e., values close to zero will indicate that the model error tends to be systematic, while large values indicate that the error is variable and randomly distributed:

Std = 1 N -1 N i=1 | [F lux sim(i) -F lux obs(i) ] -Bias | 2 .
(7)

WRF model

The WRF-ARW (Weather Research and Forecasting-Advanced Research WRF) v.4. 1.3 (Skamarock et al., 2019) was used to run the different numerical experiments. The model ran for 60 hours (from 12:00 UTC of day 18 June 2011 to 00:00 UTC of day 21 June) but only the central hours of day 19 June were evaluated. The first 21 hours of the simulation were used as spin-up and the diurnal cycle of day 20 June was not analysed due to the presence of clouds in the area, which made it more difficult to carry out the strategy designed for the main objective [START_REF] Pedruzo-Bagazgoitia | Direct and diffuse radiation in the shallow cumulus-vegetation system: Enhanced and decreased evapotranspiration regimes[END_REF]. The simulations were configured with four nested domains of 27, 9, 3 and 1 km of resolution.

The inner domain covered an area of 120x120 km 2 , but only an area of 19x19 km 2 around the central point was evaluated, which had the same grid as the AAF used for the evaluation, i.e., the centre of the central pixel was located at the same location as the centre of the AAF used for the evaluation. More details about the WRF technical configuration common for all experiments are included in Table 1. In order to add robustness to the study, three different LSMs available in WRF were analysed. The information below is mostly extracted from the literature:

1. Noah [START_REF] Chen | Coupling an advanced land surface-hydrology model with the Penn State-NCAR MM5 modeling system. Part I: Model implementation and sensitivity[END_REF]. Noah is a widely used LSM resulting from the collaboration among many different institutions. It is the default option in WRF and it is used in many other models, with an important application in operational models from the National Centers for Environmental Prediction (NCEP). The model considers four soil layers, where it computes temperature and soil moisture. It takes into account the type of vegetation (LC category), monthly vegetation fraction and soil type to calculate the runoff, ET and root zone. Since the WRF v.3.6, there is a mosaic option available in the model [START_REF] Li | Development and evaluation of a mosaic approach in the WRF-Noah framework[END_REF] to deal with the sub-grid heterogeneity (this option is not activated by default).

2. Noah-MP (Noah Multi-physics) [START_REF] Niu | The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with localscale measurements[END_REF]. It is an extension of the Noah LSM that allows the use of multiple options for land-atmosphere processes (e.g., infiltration, run-off, etc.), resulting in a total of more than 4000 combinations (the default options are used in this work). This LSM contains a separate vegetation canopy with a two-stream radiation transfer approach, shading effects and complex physics for the snow/ice processes within the soil. This LSM uses a different set of parameters for each vegetation category than Noah, with more vegetation-dependent parameters.

3. RUC (The Rapid Update Cycle) [START_REF] Smirnova | Modifications to the rapid update cycle land surface model (RUC LSM) available in the weather research and forecasting (WRF) model[END_REF]. It uses 9 soil layers with higher density close to the surface. It has a complex treatment of snow processes. In the warm season, it corrects soil moisture in cropland areas to compensate for irrigation. This model also allows a mosaic approach for the sub-grid treatment of the cell heterogeneity, but it is different than in Noah, with albedo values that correspond only to those parameters associated with the dominant LC category. The vegetation parameters used are obtained from the same look-up table than in Noah.

For each LSM, we analyse their results for the different land-cover categories of the studied area, but we ignore the results associated with the urban category, for two main reasons: 1) this category is normally problematic for the LSMs without a specific and additional urban parameterization (which is not activated in this work), and; 2) their values were estimated with a simple model in the AAFs, i.e., no measurements were available. However, due to the difficulties removing this category in a real case study, the errors are shown along the article to illustrate these issues, but the analysis of the results is out of the scope of this work. In any case, the percentage of coverage of the urban category in the analysed area was relatively small (see Fig. 2a, b, and c).

Experiments design

We performed a set of four different modelling experiments aimed at checking the sensitivity of the model to different changes in the representation of the surface. These experiments are summarised in Table 2 and explained below:

Table 2. Summary of the simulations and the names used along the article. Note how some experiments were not possible (-) for some LSMs.

DEFAULT experiment

In this experiment we used the default configuration for the surface representation in WRF, i.e., the LC dataset from IGBP-MODIS (21 categories) and a dominant approach for the sub-grid heterogeneity, which means that the fluxes were calculated taking into account only the surface parameters of the dominant LC category, i.e., assuming that no sub-grid variability exists even when this information is available. The dominant LC category is the one with the highest percentage of coverage within each pixel. This method has an intrinsic dependency on the number of LC categories of the dataset. For example, a pixel with 40% of water, 30% of conifer forest and 30% of deciduous forest will be treated as water, since it is the dominant category despite both types of forest cover 60% of the total surface of the grid. On the contrary, the dominant category would be forest if both types of forest were merged into a single category.

It is expected that these simulations will be limited by the fact that the representation of the LC in the area by IGBP-MODIS is not totally correct in comparison with the reality (shown later). Besides, the dominant approach implies that the model does not take advantage of the higher resolution of the LC datasets.

NEW-LC experiment

In this experiment, we used a more realistic LC dataset than IGBP-MODIS, obtained from 30-m resolution data prepared by CESBIO [START_REF] Inglada | Operational high resolution land cover map production at the country scale using satellite image time series[END_REF]. Figure 2a shows a satellite image of the area evaluated, which serves to visually validate the CESBIO LC at 30-m resolution (Fig. 2b). The inaccurate representation of the pixel-dominant LC by IGBP-MODIS is revealed in Fig. 2c In order to incorporate the more realistic LC from CESBIO in the WRF model, first, the different categories of the new LC dataset (CESBIO, 17 categories) were transformed into the most appropriate ones of the WRF default LC dataset (IGBP-MODIS). That is, the same 21 LC categories of IGBP-MODIS were conserved, following a similar approach as done in [START_REF] Pineda | Using NOAA AVHRR and SPOT VGT data to estimate surface parameters: application to a mesoscale meteorological model[END_REF] and in [START_REF] Schicker | Influences of updated land-use datasets on WRF simulations for two Austrian regions[END_REF]. We used the transformations indicated in Appendix A, with some special considerations:

-The two types of forest distinguished by CESBIO (conifer and deciduous) were transformed to conifer (evergreen needleleaf forest, ENF), since the EC measurements used in the AAF were taken over this type of forest due to the lack of measurements over deciduous broadleaf forest (DBF). As commented before, the measurements taken over conifers were extrapolated to the areas with deciduous forests in the AAF calculation. By converting all the forest to conifer also in WRF, the observations-model comparison was fairer.

-The two possible types of crop in CESBIO (winter and summer crop type, i.e., wheat and corn respectively) were transformed to the single cropland category available in IGPB-MODIS. Hence, we did not take advantage of the differentiation of crop type by CESBIO, but it is important to note that wheat surfaces only covered 2.8% of the analysed area (Fig. 2b).

-The four different urban categories of CESBIO were transformed to the single urban definition of IGBP-MODIS. Most of the urban surfaces of the area evaluated were defined as diffuse urban (8.2%) and industrial zones (0.8%).

This new LC information was incorporated in the modelling system following the technical details included in Appendix B.

Note how even using a more realistic LC dataset over the area, the sub-grid information was not used because the dominant approach was maintained in this experiment.

MOSAIC experiment

The sub-grid heterogeneity is important in the case of the BLLAST area due to the small scale of the LC patches (see Fig. 2b). This makes the area very appropriate for investigating the use of a mosaic approach in the model, as done in the labelled MOSAIC experiment. The mosaic approach implies that the flux is computed for each LC category (tiles) within each grid cell based on the surface parameters tabulated for each one; and then averaged taking into account the percentage of coverage of each tile. That is, taking into account the sub-grid surface heterogeneity. In this context, Noah (see [START_REF] Li | Development and evaluation of a mosaic approach in the WRF-Noah framework[END_REF], for a complete description of the Noah mosaic approach) and RUC allow the possibility of using this approach in the WRF modelling system. The LC improvements of the NEW-LC experiment were also included here. More technical details about the model implementation of this experiment are included in Appendix B.

FOREST experiment

The last experiment was motivated by the large surface fluxes errors found over those pixels covered by conifer (ENF) in previous simulations using Noah-MP. Specifically, this LSM overestimated the SH and underestimated the Le. We hypothesised that this was caused by a parametrized resistance to transpiration that was to high. Hence, we designed the FOREST experiment for Noah-MP. We modified three parameters for ENF, with the values used in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF]:

1. The slope of the Ball Berry conductance (the inverse of the stomatal resistance to transpiration) equation, the so-called MP in Noah-MP and usually known as g 1 (the slope parameter). The Ball-Berry equation [START_REF] Ball | A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions[END_REF] linearly relates the stomatal conductance to the CO 2 assimilation rate:

g s = g 0 + g 1 A n • h s c s . (8)
In this equation, the slope (g 1 ) represents the sensitivity of the stomatal conductance to assimilation, CO 2 concentration, humidity and temperature. It is the parameter that most affects the plant's transpiration [START_REF] Cuntz | The impact of standard and hard-coded parameters on the hydrologic fluxes in the Noah-MP land surface model[END_REF], which increases with larger slope values. The g 0 parameter is the minimum conductance, h s and c s are the fractional relative humidity and the CO 2 concentration at the leaf surface, respectively, and A n is the net leaf CO 2 assimilation.

MP has a default tabulated value of 6 for ENF for Noah-MP, a value that is significantly different compared to the values assigned for other categories (9 for most of them, including, for example, the broadleaf evergreen forest (BDF)). The lower value for ENF limits the transpiration processes and could be the reason for part of the Le underestimation, leading to more energy being distributed to SH. This parameter was optimised in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF] and a value of 9 was also used for ENF in their study.

2. The minimum leaf conductance or the interception in the Ball-Berry equation (g 0 , indicated as BP in the Noah-MP scheme). This parameter was changed from 0.002 to 0.01 mol H 2 O m -2 s -1 , as in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF] for ENF.

3. The maximum carboxylation rate at 25 ºC (V cmax25 ), a photosynthetic parameter that was changed from 50 to 62.5 µmol m -2 s -1 , as used in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF], which value was also selected according to previous literature. The range of observational values in their work for ENF ranged from 48 to 72 µmol m -2 s -1 and a value of 62.5 µmol m -2 s -1 was finally selected.

Hence, these three parameters were modified in Noah-MP according to [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF], following the technical details included in Appendix B. Indeed, these changes allowed for more evapotranspiration in the ENF forest, which should improve the evaluation in our case study. It should be noted that the g 1 parameter (MP) was the one which most influenced the results, as indicated in [START_REF] Cuntz | The impact of standard and hard-coded parameters on the hydrologic fluxes in the Noah-MP land surface model[END_REF] and as observed in additional previous experiments performed in our case (not shown). In any case, the objective of this experiment was to demonstrate the high impact of the associated vegetation parameters on the surface fluxes, not the optimization of these values for the specific tree species present in the area here analysed.

Pre-experiments: soil moisture and initial conditions

The previous experiments were performed with initial conditions obtained from the NCEP-FNL data (NCEP, 2000) and using the soil moisture initial conditions included in these files (1 • of horizontal resolution and without small-scale details). The impact of these initial conditions has been investigated in two pre-experiments.

One the one hand, the most controversial issue was related to the soil moisture initial value used in the model. Indeed, each LSM used in WRF should be appropriately run for a long time in order to provide appropriate soil moisture values in accordance with the specific LSM dynamics [START_REF] De Rosnay | A simplified Extended Kalman Filter for the global operational soil moisture analysis at ECMWF[END_REF][START_REF] Angevine | Land surface spinup for episodic modeling[END_REF]. This should lead to more appropriate (and realistic) initial SM values. However, as commented before, the AAF used to evaluate the model were based on the assumption of homogeneous soil moisture in the area, similar to the SM field obtained with the 1 • resolution of NCEP-FNL. Therefore, in order to perform a fairer comparison of the model with the AAF, the spin up was not used in our study to avoid the SM horizontal heterogeneity. Besides, the simulated soil moisture is, to a large extent, driven by the simulated rainfall, which might be different depending on the LSM used. This can cause undesired differences in the rainfall between the LSM during the spinup time, which is not appropriate for the comparison performed. However, the impact of performing longer spin-up times was checked in the SPIN-UP pre-experiment (spin up of 1 month).

On the other hand, we also performed some pre-experiments aimed at quantifying the impact of using different initial conditions data. The results obtained with NCEP-FNL were compared to those obtained with ERA-INTERIM, both data with similar horizontal resolution.

The impact of these two experiments related to the previous initial conditions is summarised and discussed in the results section.

3 Results: Turbulent fluxes sensitivity to surface changes

The model skill simulating the AAF fluxes is quantified in this section for the different experiments listed in Table 2. The results are shown and commented individually for each experiment. Since some of the results and the associated discussion differed significantly depending on the LSM used, the results were subdivided for each LSM. As a general reference, Table 3 shows a summary of the total scores (in the whole evaluated area) obtained for SH and Le using the different LSMs and for all the experiments performed. 3a). For the Le, the pixels covered by SaW (Savanna Woody, which represent 33% of the area) are normally characterised by a remarkable Le underestimation, with biases values around -100 Wm -2 (Fig. 3d).

The definition of some of the LC categories existing in the DEFAULT experiment did not represent appropriately the real LC of the area (as seen in Fig. 2c), but they influenced notably the simulated model errors. This is the first indicator of the high dependency of the model results on the LC categories. Then, based on the main hypothesis, the results should improve using a more realistic LC representation (NEW-LC experiment).

The all-area averaged results for SH and Le from NEW-LC-Noah are shown in Table 3. In this experiment, the LC was modified towards a more realistic representation using the high-resolution data from CESBIO. This led to five LC categories that were present in the area of study as dominant LC, with evergreen needle-leaf forest (ENF) and grass (Gra) covering the 45% and 47% of the total area, respectively. The SH biases observed in the DEFAULT experiment were improved in the NEW-LC, especially for ENF and Gra, with slightly negative (close to zero) values. This led to a slightly negative SH bias in the whole area (Fig. 3b, see horizontal red dashed line). For the Le bias (Fig. 3e), the opposite is observed, with a slight Le overestimation mainly caused by the Le overestimation over ENF (too much ET simulated by the model in this forest type).

As it will be shown later, this result contrasts with the findings obtained with the rest of LSMs.

The bias is a good indicator of systematic errors for specific LC categories, but it is a poor indicator of the total model behaviour in the area (it can mix positive and negative values leading to neutral ones). Hence, the std and the RMSE have also been included. The std was also diminished in the NEW-LC experiment, with less variability of the error in the new experiment.

In the case of the RMSE (Fig. 3c andf), the results are shown by real LC categories (not model categories), whose fractions of coverage are very similar to those used in the NEW-LC experiment (slightly different due to small differences in the categories, e.g., corn and wheat in the real LC were merged into crop in the model).

The SH-RMSE (Fig. 3c) improved for the total area when using the improved LC dataset, from 56 Wm -2 (DEFAULT-Noah)

to 44 Wm -2 (NEW-LC-Noah). The vertical bars of Fig. 3c indicate the RMSE for each real LC category (the LC categories used in the AAF). They provide information about the types of vegetation associated with larger errors, i.e., vegetation types whose processes or parameters are not well represented by the model. In this case, the SH improvements in the NEW-LC experiment (red bars) are observed for all the pixels except for the only pixel where the dominant LC is wheat, associated with a larger SH flux than corn in the observations (see Fig. 2). In any case, wheat crops only covered 0.3% of the total area as dominant LC category (one pixel) and, therefore, the contribution to the total RMSE values was small. The results for Le-RMSE (Fig. 3f) also show significant improvements when using the improved LC dataset, except for those pixels covered by urban (with a small impact in the total values due to their scarce presence (2.4%) as dominant pixels and which results are out of the scope of this work). The Le improvement was more substantial in those pixels covered by grass than by forest, since Le is slightly overestimated in the forest pixels (see ENF in Fig. 3e). Notice how besides the better LC pixels distribution in the NEW-LC experiment, we also avoided using some LC categories associated with large errors in the DEFAULT-Noah experiment, such as the unrealistic SaW pixels related to large Le underestimation (Fig. 3d).

-

Noah-MP

Most of the LC categories showed a positive SH bias for DEFAULT-Noah-MP (Fig. 4a) and a slight negative Le bias (Fig. 4d), especially for those pixels whose dominant LC was some of the different types of forest (ENF, WiF, DBF, MixF or SaW). The std was the largest for SH and the smallest for Le among the different LSMs (see DEFAULT-Noah-MP in Table 3) This led to a large SH error and to the smallest error amongst the LSMs compared for Le (despite the slight underestimation), with RMSE values of 96 Wm -2 and 46 Wm -2 , respectively (see Fig. 4c and f and also Table 3).

Contrary to what happened with Noah, the results were, in general (all domain), worsened in the NEW-LC experiment, even for the case of SH where the RMSE was already high. The SH bias remained positive, mainly influenced by the significantly high SH bias in ENF pixels (more than +150 Wm -2 ), while the bias over the grass and crop pixels were close to zero (Figure 4b). For the Le, important negative biases were observed in those pixels mostly covered by ENF (around -120 Wm -2 ), but with lower errors for grass or crop pixels. The contrasting differences among the biases for the two more common LC categories (grass and forest) led to high values of std. The urban pixels were characterised by important biases but their small proportion as dominant LC category led to a small contribution to the total error. This error variability among categories and the fact that the 45% of the total area was covered by ENF led to final RMSE values of NEW-LC-Noah-MP (horizontal red dashed line in Fig. 4c) significantly higher than those of DEFAULT-Noah-MP (shown in blue), increasing the total SH-RMSE up to 116 Wm -2 and the Le-RMSE up to 90 W m -2 , with a significant worsening for the forest and urban pixels. However, the bars in Fig. 4c and f also illustrate the improvement in the crop, moor and grass pixels for SH and grass for Le, but they were not enough to compensate for the large errors found over the ENF pixels. This was not observed in the Noah experiment (previous subsection), which uses a different set of vegetation parameters than Noah-MP for each LC category. This suggests that the biased results for Noah-MP could be partially caused by some vegetation parameters that were inappropriate for the ENF LC category, hypothesis investigated latter in the FOREST experiment.

-

RUC

The SH and Le biases obtained when using the DEFAULT-RUC simulation are the highest ones among the compared LSMs (see comparative Table 3), while the std values are similar than for the other LSMs. Hence, RMSE values of 97 Wm -2 and 81 Wm -2 for SH and Le, respectively, were obtained for the DEFAULT-RUC simulation. There is a systematic SH overestimation for all the LC categories (Fig. ??a) for DEFAULT-RUC, which is especially aggravated by those pixels unrealistically represented by WiF (Wild forest), with a bias of more than 100 Wm -2 which covered an important fraction of the total analysed area (27%).

The Le is systematically underestimated (Fig. ??d) for all the LC categories, which shows a tendency towards too little ET in this LSM in this studied case, especially in those LC categories with shorter vegetation (SaW, Sav and Gra).

These biases, the std and the RMSE values were very similar in the NEW-LC-RUC experiment, using the more realistic LC, with SH-RMSE of 90 Wm -2 and Le-RMSE of 85 Wm -2 . The values obtained for each LC category (Fig. ??b and e) were within the same order of magnitude than those obtained in the DEFAULT-RUC experiment, or even higher, as it is the case of the SH simulated in the grass surfaces, with a remarkable overestimation of around +75 Wm -2 . This led to very similar (and high) RMSE values for the DEFAULT and NEW-LC experiments (Fig. ??c andf). Although the SH simulation improved for the grass and corn pixels, it worsened for the forest pixels. For the Le, the contrary was observed: the biases of the forest and crop pixels improved but worsened for the grass ones, which were associated with too low ET. Note how the RMSE in the NEW-LC experiment over the grass surfaces was higher than 100 Wm -2 in RUC while they were around 40 Wm -2 in the rest of the analysed LSMs, where normally the grass surfaces were associated with improvements.

Although it is not shown in this work, we detected a higher sensitivity of this LSM to the soil type in comparison with the other LSMs. Thus, the parameters associated with clay loam could be inappropriate for this LSM, leading to a wrong partitioning of the net radiation in the model (too low Le and too high SH). In contrast, the sensitivity of RUC to the LC categories was the lowest one among the compared LSMs. In any case, the analysis of the relative contribution of the soil type on the surface fluxes is out of the scope of this study. Besides, the soil moisture values used to initialise this LSM should be corrected due to its different SM dynamics and baseline (see the discussion included in the pre-experiments shown in Section 3.5).

- the more the percentage of forest, the higher the SH, since forests were associated with the highest values of observed SH.

The opposite is observed for the grass pixels, with SH diminishing for increasing percentages since grass EC measurements showed the lowest SH values. These dependencies are slightly observed for the Le (Fig. 5h (Le)) due to the similar Le values measured over all the surface types with the EC towers (see Fig. 1b andd). Coming back to Fig. 5, one can deduce that Noah is the LSM that provides the most realistic values in comparison with the AAF from a visual comparison, especially for Le.

The other LC category associated with large biases was the urban one. In this case, the effect on the scores of the total area was minimum, since only the 2.4% of the pixels have urban as dominant LC category. As commented before, a fair discussion about the absolute values of the errors found over this type of LC is out of the scope of this article due to the relatively high uncertainties in the AAFs and the lack of a specific urban parameterization in the model, with too extreme values for Noah and Noah-MP: 0 Wm -2 for Le, see grey circles in Fig. 5e andf. Furthermore, the dependence of the AAF with the percentage of coverage of the dominant LC in each pixel is well observed in Fig. 5d andh, which is inherent to the methodology used in the AAF calculation. However, this dependence is hardly seen from the model outputs (as expected using the dominant approach, note the small slope of the scatter plots for each coloured category). As commented before, it should be noted that, by default, these LSMs use a dominant approach for each grid cell, using some surface information (LAI, roughness length, or albedo) only from the dominant LC [START_REF] Li | Development and evaluation of a mosaic approach in the WRF-Noah framework[END_REF]. This led to well-differentiated fluxes values for each category (Fig. 5).

Hence, from the NEW-LC experiment, we can conclude that there is a high dependence of the fluxes on the LC type in WRF, which can lead to important biases if the parameters and processes associated with specific categories (conifer forest in this case) are not appropriate. These results agree with those found by [START_REF] Couvreux | Boundary-layer turbulent processes and mesoscale variability represented by numerical weather prediction models during the BLLAST campaign[END_REF], where twelve IOPs of the same field campaign were simulated with ARPEGE [START_REF] Courtier | A global numerical weather prediction model with variable resolution: Application to the shallow-water equations[END_REF], AROME [START_REF] Seity | The AROME-France convective-scale operational model[END_REF], and ECMWF models [START_REF] Simmons | The ECMWF medium-range prediction models development of the numerical formulations and the impact of increased resolution[END_REF], with grid sizes of 10 km, 2.5 km and 16 km, respectively. In this work, the SH was systematically overestimated over the area analysed by the three models. Indeed, a larger SH overestimation was found for the ARPEGE model in two of the three evaluated pixels close to the area, which surface was represented as forest. This also highlights the issues found over this type of vegetation cover, which was enhanced with a dominant approach (as in WRF in this work and The scores obtained using the mosaic approach in Noah and RUC LSMs are shown in Table 3, with significant improvements in comparison to the NEW-LC experiment in Le. For Noah, a reduction of almost 20 Wm -2 in the Le-RMSE (from 63 Wm -2 to 45 Wm -2 ) was observed, result of a decrease in the bias and in the std. For RUC, the improvement in the Le-RMSE was of 9 Wm -2 . The SH scores are almost no modified in Noah and worsened in RUC-MOSAIC.

For the case of Noah, Fig. 6 shows the simulated SH and Le obtained from the NEW-LC and the MOSAIC experiments. The use of the mosaic approach (panels b and e) caused a merging effect among values from different categories on the simulated fluxes in comparison with the dominant one (panels a and d), with results that better resemble to those figures obtained from the AAF (panels c and f). This is particularly clear for SH, where the stronger dependence of the fluxes with the type of dominant LC (Fig. 6a) is removed in the mosaic approach (Fig. 6b). The impact was larger on those pixels where the dominant LC was forest, which were associated with more extreme fluxes values. Note also the near zero values of Le in the dominant approach in urban pixels (panel d) and the Le values that depend on the fraction of urban within each pixel (panel e). The correlations of the fraction of coverage of the dominant LC category and the surface fluxes are included in the legend for each LC category in Fig. 6. This correlation depends on the strength of the fluxes values associated with the specific categories. For example, the anti-correlation observed in Fig. 6e for the urban pixels is due to the strong impact when increasing the percentage of cover of this LC category. Thus, if a pixel with 100% urban were present, its Le value would tend to 0 Wm -2 using a mosaic approach, as observed when using the dominant method in Fig. 6d.

However, this merging effect was only slightly observed for the case of the MOSAIC experiment in RUC (Fig. 7). This was probably caused by the fact that the mosaic approach used in RUC is applied for different variables than in Noah. In the case of RUC, only average emissivity, LAI, roughness length and plant resistance are calculated based on the percentage of each land cover. In RUC, the averaged albedo, which is the parameter that has the highest impact on the net radiation of each grid cell, subsequently affecting SH and Le, is not used. Figure 8 shows the albedo differences between NEW-LC and MOSAIC used by the Noah (upper figures) and RUC models (bottom figures). While the values for MOSAIC-Noah (Fig. 8b) consisted of a weighted average from the different LC of each grid cell, it was not the case for MOSAIC-RUC (Fig. 8d), which diminished the impact of the mosaic approach application in the RUC model.

Adding a mosaic approach caused a change in the LC percentage technically used for some LC categories. This is well observed in the comparison of the percentages shown in Fig. 1b andd. For example, 2.2% of the pixels were characterised as urban with the dominant approach (panel d) while the coverage using a mosaic approach increased up to 9% (panel b).

In this case, the mosaic approach increases the percentage of urban fluxes contributing to the averaged values, although in a much more diffuse way (concentrated in more pixels). However, the extreme effect of those pixels fully considered as urban are removed. This is also the case for the crop pixels, that cover 9.9% of the area (wheat and corn) shown in Fig. 1b (taking into account the subgrid variability), but only 3.9% as dominant category in 1x1 km pixels (Fig. 1d).

In any case, the mosaic approach might be always more realistic than the dominant one; in the latter, the sub-grid information is not used and a unique LC type is defined, even when the combination of secondary and similar LC categories have a higher percentage than the dominant one (e.g., 40% of conifer, 30% of dense shrub and 30% of open shrub). If the resulting dominant category is associated with inappropriate parameters (e.g. conifer in the example and in our case), the error will be greater in the dominant approach. For this reason, it is also crucial to be cautious with some LC definitions that could lead to too extreme fluxes. This can have an important impact on the results of the simulations, either with a dominant or a mosaic approach, as stated in [START_REF] Mallard | Effects of Mosaic Land Use on Dynamically Downscaled WRF Simulations of the Contiguous United States[END_REF]. Furthermore, appropriate fluxes measurements for model evaluation over urban surfaces still remain a challenge, and the LSMs present issues without an additional urban parameterization included, aspect that should be better analysed in future studies. . Same as in Fig. 5 but comparing the results from Noah-NEW-LC (using a dominant approach) and from Noah-MOSAIC (using a mosaic approach). Panels c and f corresponds to the area-averaged fluxes (AAF) from observations, considered as the reference. The correlation between the fluxes and the fraction of coverage of the dominant LC category is included in the legend for each LC category within each pixel, with the number of pixels used in brackets.

FOREST experiment

The previous experiments revealed the large biases in simulated surface fluxes associated with conifer forests (ENF) for Noah-MP, with too low simulated values of Le and too high SH in comparison with the observed values (Fig. 5). This was also the case for RUC (Fig. 5d andi), but in this case the biases were not only limited to the ENF pixels but to all the categories.

Since the biases were only observed in ENF for Noah-MP, we hypothesised that this was caused by a too high resistance to the transpiration parameterised for these type of trees. We explained in Section 2.3.4 the motivation for this experiment and we justified the changes applied, based on the parameter values used in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF] for ENF forest. Specifically, we changed the values of the g 1 , g 0 and V cmax25 parameters. These modifications should lead to increased transpiration, and therefore, also to decreased SH. Indeed, this is observed in the results, the RMSE decreases to 77 Wm -2 for SH and to 56

Wm -2 for Le, which are lower errors compared to the NEW-LC-Noah-MP experiment (see Table 3). This was the result of an improvement in the systematic error found in the ENF pixels, due to the reduction in the extreme biases found in ENF.

Specifially, the ENF mean bias reduced from +156 Wm -2 to +85 Wm -2 for SH and from -110 Wm -2 to -32 Wm -2 for Le comparing the NEW-LC and the FOREST experiments (these results are not shown in figures), while the std values remained similar, illustrating how we improved the systematic bias associated with this LC type due to incorrect parameters related to the evapotranspiration of these trees in this case study.

The results are consistent with the fluxes shown in Fig. 9, where the SH and the Le observed at the forest site (black line)

were compared to the simulated values obtained with the central pixel of the domain completely (100%) covered by ENF (conifer, blue), or by ENF with these parameters changed (red line). These two additional simulations shown in Fig. 9 would correspond to the NEW-LC Noah-MP and the FOREST Noah-MP experiments, respectively. However, in these graphics the effect of the change is analysed in a clearer way for the whole diurnal cycle. Again, it is demonstrated how the tuning of these parameters towards less resistance to transpiration provided better results than the original parameters tabulated for ENF, both for SH and Le. Specially, the slope of the Ball-Berry equation (g 1 ) is the parameter that most influenced the results, which does not seem appropriate in this area for this type of trees, since it limits transpiration too much. As asserted in [START_REF] Medlyn | Reconciling the optimal and empirical approaches to modelling stomatal conductance[END_REF], g 1 is the key parameter for plant stomatal conductance, being quite variable among species in areas with different environmental conditions. The authors suggested that g 1 should increase with increasing temperatures, as might be the case 3.4 Analysis of the simulated evaporative fraction (EF)

The previous analyses provided information about the biases and the RMSE of SH and Le independently for the different experiments and LSMs. However, these analyses did not take into account the total energy available to be partitioned into the atmospheric fluxes, which may be different depending on the surface type. Hence, in Fig. 10 In Noah (Fig. 10a), the DEFAULT experiment (blue) provided a large variation in the EF of each real LC category. This resulted from the unrealistic and varied LC representation in this experiment. The NEW-LC experiment (red) partially corrected the large EF variability, especially for the grass and forest pixels. Therefore, in general, not only the SH and Le biases were corrected, but also the fluxes partitioning, except for the moor and urban pixels, with scarce presence. Besides, the MOSAIC experiment (in orange) served to improve the simulated values of EF, which were very close to those obtained from the AAF.

The DEFAULT experiment in Noah-MP (Fig. 10b) systematically underestimated the EF (blue) in comparison to the observations (dark grey). The NEW-LC experiment (red) served to correct this underestimation in all LC categories except in the forest and urban grid cells. As discussed before, the Le was significantly underestimated in the forest (ENF) pixels, while the SH was overestimated, which led to underestimation of EF. This was partially corrected in the FOREST experiment (represented with green in Fig. 10b) by applying modified parameters that made easier the transpiration from these plants. However, even with these changes, the EF values obtained from the model were far from those from the observations. RUC (Fig. 10c) exhibited the smallest sensitivity to the modifications performed in each experiment in comparison to the other LSMs. As commented before, the LC changes in RUC had a lower impact on the results, which where more impacted by the soil type (not shown) than by LC. Hence, RUC showed a general tendency to underestimate EF with systematic Le underestimation and SH overestimation for all the LC categories, as it was previously observed. to -21 Wm -2 for SH and Le, respectively), but at the expenses of an increase of the std, which led to similar RMSE values. This LSM showed biases with different sign depending on the LC category (see Figure 4), especially for forest (SH overestimation and Le underestimation) and grass (the contrary), which led to a mean bias value that should be analysed with caution: the value in the whole area, as shown in Table 4, includes compensating errors the can conduct to erroneous conclusions. In the case of RUC, the SPIN-UP experiment improved the biases for both fluxes (from 77 Wm -2 to 47 Wm -2 and from -70 Wm -2 to -46 Wm -2 for SH and Le, respectively) while the standard deviation of the error is the same, leading to a substantial reduction of the RMSE. Indeed, this was suspected in Reviewer 1 (2021) since the RUC model has a different baseline and dynamics for the soil moisture than Noah, Noah-MP and the NCEP-FNL data (based on Noah). In this case, leading the model to auto-spin up some time is convenient, since the fluxes simulated by this LSM improved with their RUC-specific SM values. Hence, the SM initial value for RUC should be higher than for Noah (in this case study), which corrects the Le underestimation and the SH overestimation observed in Fig. ??.

ERA-INTERIM experiment

Not only the SM can have an important effect on the model results, the rest of the atmospheric variables used to initialise the model can also impact the whole simulation. This possibility was studied using a different database to initialise the model (ERA-INTERIM versus NCEP-FNL). The comparison of the scores obtained with the ERA-INTERIM experiment slightly affect the evaluation: in Noah the bias changed from -17 Wm -2 to -16 Wm -2 for SH and from 34 Wm -2 to 45 Wm -2 for Le, leading to similar RMSE values. In the case of Noah-MP, the bias varied from 66 Wm -2 to 75 Wm -2 for SH and from -42 Wm -2 to -39 Wm -2 for Le, also with similar RMSE values. Therefore, although there is an impact of the initial conditions on the evaluation of the model, in this case the effects were not very large. Nonetheless, using other datasets with higher resolution (as ERA-5) could lead to improvements due to the improved resolution used at the model initialization, but additional experiments aimed at investigating the impact of the resolution of the data used to initialise the model are out of the scope of this article.

Summary and conclusions

The changes in the LC of the Earth's surface trigger varied and sometimes unpredictable consequences at different spatiotemporal scales, affecting biophysical processes in the soil and in the atmosphere. Hence, it is crucial to know more about the impacts of the LC changes on all these processes. In this work, we investigated the sensitivity of turbulent heat fluxes simulated by the WRF model to the manner in which the surface is represented in it.

To this aim, different sensitivity experiments were performed for a case study over an heterogeneous area in the south of France. They were evaluated with gridded area-averaged fluxes (AAF), computed from tower measurements installed over five vegetation types (forest, corn, wheat, grass and moor) during the BLLAST field campaign. In order to add robustness to the study and to detect differences, the experiments were carried out using three LSMs available in WRF: Noah, Noah-MP and RUC.

First, a control experiment was performed with the default options in WRF: LC from the IGBP-MODIS database and a dominant sub-grid approach, i.e., the model used the tabulated surface parameters of the LC category with the highest percentage of coverage in each pixel. We hypothesised that these simulations were limited because of the large differences between the LC representation and the actual surface, and because of the loosing of LC information at the sub-grid scale.

Thus, a new experiment was designed (named NEW-LC), which improved the surface representation by adding the LC information from the 30-m resolution CESBIO maps, which were much more accurate and realistic in the area than the default IGBP-MODIS dataset. The observed changes in the surface fluxes were dependent on the LSM used, due to their differences in the parameters associated with each vegetation type, and also to their different representation of the surface processes. The improvement was clear for Noah for all the LC categories. RUC was the LSM that showed the weakest response of the fluxes to the LC categories, without substantial changes in the scores. Noah-MP showed some improvements in those pixels covered by crops or grass, but they also exhibited an important SH overestimation and Le underestimation in the surface fluxes simulated over those pixels mainly covered by conifer forest (ENF). The ENF biases contributed significantly to the total model error due to their relatively high percentage of coverage (45%) in the analysed area. The NEW-LC experiment revealed the need for a correct representation of LC in the analysed area, in part due to the high dependency of the fluxes on the LC categories. In addition, the appropriate characterization of surface parameters associated with some LC categories (e.g., conifer) still needs to be improved, as it was also discussed in previous works (e.g., [START_REF] Li | Development and evaluation of a mosaic approach in the WRF-Noah framework[END_REF][START_REF] Cuntz | The impact of standard and hard-coded parameters on the hydrologic fluxes in the Noah-MP land surface model[END_REF].

In the second experiment (named MOSAIC), the sub-grid heterogeneity (below 1 km) was taken into account with a mosaic approach in Noah and RUC, meaning that the fluxes in each grid were calculated as weighted averages from individual surface fluxes obtained from each tile or LC category. The mosaic approach caused more homogeneity among the surface fluxes simulated in the analysed pixels, which corresponded better to the AAF used as benchmark data. This improvement in the surface representation led to improved scores for Noah (especially for Le) while smaller changes were observed in RUC. This was because in RUC the mosaic approach did not include the use of pixel-averaged albedo values based on the percentages of each LC category, as done for the surface roughness, LAI and emissivity. In Noah, the albedo was also averaged, significantly contributing to the improvements, since the albedo is the parameter that seems to have a larger impact on the net radiation available to be partitioned into SH and Le.

Finally, a last experiment (named FOREST) was motivated by the issues found in the conifer pixels for Noah-MP. The modifications in the FOREST experiment were conducted to reduce the resistance of conifer trees to transpiration, using updated parameters as used in [START_REF] Bonan | Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil-plant-atmosphere continuum[END_REF]. The effect of these changes was to facilitate the transpiration processes, coinciding better with the observations and improving the scores.

Two additional pre-experiments aimed at checking the model sensitivity to the initial conditions were performed. In the first one, a spin-up period of 1 month was applied for each LSM in order to obtain more appropriate SM values at each pixel. This was very important for the case of RUC since the SM in this LSM has a different dynamics and range of values than Noah, Noah-MP and the SM from NCEP-FNL, leading to improved scores with a corrected SM value. In the second pre-experiment the sensitivity to the initial conditions was checked, comparing the NCEP-FNL results with those using ERA-INTERIM, with a slight impact on the simulation scores.

This work demonstrates again the importance of a correct representation of LC in the area which is evaluated, as it was also shown in previous works [START_REF] Cheng | Investigation of the effects of different land use and land cover patterns on mesoscale meteorological simulations in the Taiwan area[END_REF][START_REF] Schicker | Influences of updated land-use datasets on WRF simulations for two Austrian regions[END_REF][START_REF] Jiménez-Esteve | Land use and topography influence in a complex terrain area: A high resolution mesoscale modelling study over the Eastern Pyrenees using the WRF model[END_REF]. This can considerably affect the simulation of the fluxes that will drive the associated boundary-layer processes. Besides, it is worth using a mosaic approach to benefit from the sub-grid surface information that is normally available.

In any case, the particular conditions of the region and dates analysed make that some specific conclusions might not be applicable to other regions or even to the same region under different conditions. For example, it is also possible that in our case we observed more ET over the conifer trees than in the model due to the possible particularities of the area. Hence, the parameters adjusted in the model for the conifers could be due to the differences among species belonging to the same LC category [START_REF] Granier | Variation of hydraulic conductance of some adult conifers under natural conditions[END_REF], tree density, tree age [START_REF] Sellin | Hydraulic and stomatal adjustment of Norway spruce trees to environmental stress[END_REF], or to the specific surface conditions of this case study The possible uncertainties in the EC measurements used to evaluate the model should be also taken into account, especially over those vegetation types where it is, somehow, more difficult to have accurate or representative high-frequency measurements, as in the case of the forest. Furthermore, the calculation of the AAF consists on important assumptions based on the spatial extrapolation of EC data that can add uncertainty to the data used for the evaluation. In our case, they were also constructed using a simple estimation of the fluxes for urban surfaces, due to the lack of measurements. This could be also associated with errors (although with a small percentage of coverage in the analysed area). All these necessary simplifications that were done highlight again the importance of having extensive measurements over a wide variety of surface types [START_REF] Cuxart | Evapotranspiration over Land from a Boundary-Layer Meteorology Perspective[END_REF] and including atmospheric, soil measurements, and those related to the plant physiology and status.

Appendix A

Table A1 shows the LC transformation performed in the NEW-LC experiment. [START_REF]NCEP FNL Operational Model Global Tropospheric Analyses, continuing from[END_REF](NCEP, 2000)), which were used to initialise the WRF model. Thanks to CESBIO for the land cover data and for the soil moisture data from the DISPATCH soil-moisture product used to test some initial conditions, especially Olivier Merlin and Nitu Ojha. Finally, we would like to thank Wayne Angevine (reviewer 1) and the anonymous reviewer 2 for

Figure 1 .

 1 Figure 1. a) Observed sensible heat flux (SH) on 19 June 2011 during the BLLAST campaign over different vegetation types. b) Idem for the latent heat flux (Le). c, d) Same as in a and b, but with the time series smoothed for a better visualization of the main differences between vegetation types. The fluxes values were calculated using the eddy covariance technique with averaging windows of 30 minutes.

  and conifer transpiration increased

Figure 2

 2 Figure 2. a) Satellite picture of the area analysed (From ©Google Earth, image Landsat/Copernicus) b) CESBIO land cover map from 30-m resolution data. c) Dominant land cover (1 km) used in the default experiment (DEFAULT), from the IGBP-MODIS database. d) Dominant land cover (1 km) used in NEW-LC experiment, obtained from the CESBIO land cover map shown in b. The black rectangle indicates the area with soil type dominated by loam, which was not used for the evaluation. The rest of the area was characterised by clay loam (dominantly), as the pixels including the sites where the EC measurements were taken (maximum 4 km away from the central point of the BLLAST field campaign, indicated with a black x symbol). The complete names corresponding to each LC category abbreviation are included in Appendix A.

3. 1

 1 DEFAULT experiment versus NEW-LC experiment 3.1.1 NoahThe bias of the simulated SH and Le by the DEFAULT-Noah simulation relative to the observations (the AAF) is shown in Fig.3a and drespectively. The results are presented for each LC category, taking into account the dominant LC of each grid cell (1 km) in the model. Note that the DEFAULT experiment includes ten LC categories in the area, based on the IGBP-MODIS database. The total bias values for SH and Le (represented with blue dashed horizontal lines) are close to 0 Wm -2 due to the compensation effect when merging cells with positive and negative biases, but the std has significant values due to these differences among LC categories, with some surfaces that systematically present noticeable biases. Those pixels represented as WiF (Wild forest) in the model (27% of the total area) normally present an overestimation of the SH, while the Mixed forest

Figure 3

 3 Figure 3. a) Sensible heat fluxes (SH) bias boxplots by model LC category (x-axis) for the default Noah simulation (DEFAULT-Noah).Values are calculated for all the pixels with the specific LC category as dominant in the entire analysed area from 09:00 UTC to 15:00 UTC (hourly values). The standard deviation (std) is indicated with stars for each LC category. Horizontal coloured dashed line indicates the mean bias in the area analysed and the black pointed line for the std (these values are also indicated in the legend). The horizontal black solid line indicates the line of 0 Wm -2 bias. The boxplots indicate the 50% of the distribution (from percentile 25 to percentile 75), the upper and lower whiskers the rest of the distribution and the red crosses the outliers. Note that these figures include the variability due to the number of pixels with the same LC category as dominant (variable depending on their coverage) and also the variability due to the seven hours analysed. b) Idem for the results obtained from the Noah simulation with the improved LC (NEW-LC-Noah). c) SH root-mean-square error (RMSE) for DEFAULT-Noah and NEW-LC-Noah for each real LC in the area. Horizontal dashed lines indicate the mean RMSE for the whole analysed area. d, e, f) Idem than a, b, c but for the latent heat flux (Le). The percentages included in the x-labels refer to those covered by each category with respect to the total evaluated area.

Figure 4 .

 4 Figure 4. Same as in Fig. 3 but for Noah-MP

Figure 5 .

 5 Figure 5. Same as in Fig. 3 but for RUC

Figure 6 .

 6 Figure 6. SH (upper figures) and Le (bottom figures) simulated at 12:00 UTC by each LSM: Noah (a and e), Noah-MP (b and f), and RUC (c and g). The x-axis shows the total fraction covered by the dominant LC category in each respective pixel. The results are provided for the different LC categories with colours. Panels d and h represent the same but for the gridded area-averaged fluxes (AAF) calculated from the EC measurements, which are used for the model evaluation.

  Figure7. Same as in Fig.5but comparing the results from Noah-NEW-LC (using a dominant approach) and from Noah-MOSAIC (using

Figure 8 .

 8 Figure 8. Same as in Fig. 6 but for NEW-LC-RUC (dominant) and MOSAIC-RUC (mosaic).

Figure 9 .

 9 Figure 9. Model albedo over the analysed area. Comparison of the grid values between the dominant (left) and the mosaic (right) approaches for Noah (a, b) and for RUC (c, d). Note the absence of change in the mosaic approach using RUC.

Figure 10

 10 Figure 10. a) Sensible heat flux (SH) from two simulations with the evaluated pixel completely covered by the original ENF (conifer, blue) and ENF-mod (conifer with parameters modified, red). The observations taken over the forest site are indicated with black line. b) Idem for the latent heat flux (Le). The model results were obtained using the Noah-MP land-surface model for day 19 June 2011. The biases and root-mean-square error (RMSE) are included.

  we represented the evaporative fraction (EF = Le / (Le+SH)) for each LSM (panels), sorted by real LC categories of each grid cell for the different experiments performed (DEFAULT in blue, NEW-LC in red, MOSAIC in orange and FOREST in green). The values obtained from the AAF are also included as reference with dark grey colours in the different panels. The highest observed values were observed in pixels with grass as dominant LC, while in the forest ones the EF was around 0.5. The values observed from the EC towers (representing an homogeneous surface) are indicated with black circles and would represent what it would be expected in a pixel completely covered by each LC category. For example, EF is 0.73 for grass surfaces and 0.45 for forest ones.

Figure 11 .

 11 Figure 11. Evaporative fraction (EF = Le/(SH+Le)) for the different LSMs used: a) Noah; b) Noah-MP; and c) RUC. The results are ordered by real dominant LC type of the pixels (vertical subdivisions) and for the different experiments performed: DEFAULT (blue), NEW-LC (red), MOSAIC (orange, only in Noah and RUC), and FOREST (green, only in Noah-MP). The values obtained from the observations (AAF) are indicated with dark-grey boxplots. Black circles indicate the EF obtained from the EC towers for each corresponding (homogeneous) surface.

  (for example with relatively high values of SM). All these aspects open an interesting new line of research with the objective of improving the parameters associated with each vegetation type, which could be achieved, among others, by including leaf-level measurements of stomatal conductance in experimental campaigns, as recently stated in Vilà-Guerau de Arellano et al. (2020).

Table 1 .

 1 Details about the WRF model configuration common to all simulations. *National Centers for Environmental Protection (NCEP)

	FNL (final) Operational Model Global Tropospheric Analyses, continuing from July 1999 (NCEP, 2000).
	Model version	WRF-ARW v.4.1.3
	Number of domains	4
	Resolution of domains	27/9/3/1 km
	Initial and boundary data	NCEP-FNL* data (1°), each 6 hours
	PBL scheme	Yonsei University (YSU, Hong et al. (2006))
	Surface-layer scheme	MM5 similarity (Jiménez et al., 2012)
	Land-surface models	Noah / Noah-MP / RUC
	Microphysics scheme	WRF Single-Moment 3-class (Hong et al., 2004)
	Long-wave radiation scheme Rapid radiative transfer model (RRTM, Mlawer et al. (1997))
	Short-wave radiation scheme Dudhia (Dudhia, 1989)
	Number of vertical levels	40
	Time step	90/30/10/3.3 s
	Model initial date	18 June 2011 at 12:00 UTC
	Period analysed	19 June 2011 (09:00 -15:00 UTC)
	Leading time (spin-up)	21 h
	2.2.1 WRF land surface models (LSMs)

Table 3 .

 3 Summary of the scores in Wm -2 calculated for the evaluated area for each LSM (columns) configured for each experiment (raws).The first value refers to the sensible heat flux (SH), and the second one to the latent heat flux (Le). Some experiments were not possible (-) for some LSMs. The used scores are: the bias (systematic error), standard deviation (std) (random error) and the root mean square error (RMSE).

		Bias (Wm -2 ) -SH/Le	Std (Wm -2 ) -SH/Le	RMSE (Wm -2 ) -SH/Le
	Experiment	Noah	Noah-MP	RUC	Noah	Noah-MP	RUC	Noah	Noah-MP	RUC
	DEFAULT	-4 / -8	69 / -21	75 / -68 59 / 78	69 / 43	62 / 52 56 / 76	96 / 46	97 / 81
	NEW-LC	-17 / 34	66 / -42	77 / -70 42 / 54	93 / 81	48 / 54 44 / 63	116 / 90	90 / 85
	MOSAIC	-16 / 18	-	84 / -61 44 / 42	-	50 / 51 43 / 42	-	97 / 76
	FOREST	-	38 / -11	-	-	66 / 54	-	-	77 / 56	-

(MixF) pixels (14%) underestimate it (Fig.

Table A1 .

 A1 LC categories of CESBIO and respective transformation to the IGBP-MODIS LC categories. The respective codes of each dataset are included in the central columns. * All crop types were transformed to the single cropland category available in IGBP-MODIS; however, most of the crop types in the area were summer crop (mainly corn). ** Deciduous forests were transformed to the LC category of evergreen needleleaf forest (conifers), even when a deciduous broadleaf forest category is available in IGBP-MODIS. This was done due to the lack of measurements over deciduous trees, which made that the area-averaged maps used to evaluate the model were constructed with data from the conifers; this strategy allowed a fairer model-observation comparison. Abbreviations used along the article in the text and figures are indicated with brackets. FR infrastucture Reasearch. Thanks also to the OCASO project (ref. 0223_OCASO_5_E) for the access to the server used to run the WRF model. Thanks to the National Centers for Environmental Protection (NCEP) for the NCEP FNL (final) Operational Model Global Tropospheric Analyses data, continuing from July

	leans on ACTRIS-
	Name and code in CESBIO		Code and name in IGBP-MODIS
	Summer crop*	11	12	Cropland (Cro)
	Winter crop	12	12	Cropland (Cro)
	Deciduous broadleaf forest**	31	1	Evergreen needleleaf forest (ENF)
	Evergreen needleleaf forest	32	1	Evergreen needleleaf forest (ENF)
	Grass	34	10	Grassland (Gra)
	Moor	36	6	Closed shrublands (CShr)
	Dense urban	41	13	Urban (Urb)
	Diffuse urban	42	13	Urban (Urb)
	Industrial and commercial areas	43	13	Urban (Urb)
	Roads	44	13	Urban (Urb)
	Mineral surfaces	45	16	Barren or sparsely vegetated
	Beaches and dunes	46	16	Barren or sparsely vegetated
	Water	51	17	Water
	Glaciers and snow	53	15	Snow and ice
	Prairies	211 10	Grassland (Gra)
	Orchards	221	8	Woody savanna (WSa)
	Vineyards	222	8	Woody savanna (WSa)

their comments and suggestions during the revision of this manuscript.

Acknowledgements. Carlos Román-Cascón work was funded through a Postdoctoral Grant from the Centre National d'Études Spatiales (CNES). The BLLAST field experiment was made possible thanks to the contribution of several institutions and supports: INSU-CNRS (Institut National des Sciences de l'Univers, Centre national de la Recherche Scientifique, LEFE-IDAO program), Météo-France, Observatoire Midi-Pyrénées (University of Toulouse), EUFAR (EUropean Facility for Airborne Research) and COST ES0802 (European Cooperation in the field of Scientific and Technical). The field experiment would not have occurred without the contribution of all participating European and American research groups, which all have contributed in a significant amount. BLLAST field experiment was hosted by the instrumented site of Centre de Recherches Atmosphériques, Lannemezan, France (Observatoire Midi-Pyrénées, Laboratoire d'Aérologie). The 60m tower is partly supported by the POCTEFA/FLUXPYR European program. BLLAST data are managed by SEDOO, from Observatoire Midi-Pyrénées and maintained by the French national data center Data Terra/AERIS. The French ANR (Agence Nationale de la Recherche) supported BLLAST analysis in the BLLAST_A project. The authors acknowledge the support of the U.S. National Science Foundation grant number PDM-1649067 and the Spanish governments project CGL2015-65627-C3-3-R and CGL2016-75996-R (MINECO/FEDER). This

Code availability.

Data availability.

Code and data availability. The source code of the Weather Research and Forecasting model (WRF v4.1.3) is available at https://github.com/wrfmodel/WRF/releases (last access: January 2021). The initial and boundary data used for the specific analysed period can be downloaded

in AROME in [START_REF] Couvreux | Boundary-layer turbulent processes and mesoscale variability represented by numerical weather prediction models during the BLLAST campaign[END_REF]). This simplified approach does not take advantage of the available sub-grid surface information and motivated the MOSAIC experiment, where the sub-grid heterogeneity was taken into account in the Noah and RUC models.

As commented in the data section, the impact of the data used to initialise the model has also been investigated and quantified.

These are named pre-experiments because their impact is mainly limited to the initial values used to initialise the model (i.e., they do not alter surface parameters that affect the whole simulation as in the experiments analysed before). Two preexperiments were designed: 1) SPIN-UP, in which the WRF model was run with each LSM for a longer time (1 month) before the analysed period, allowing to obtain SM values more appropriate for the specific dynamics of each LSM, and; 2) ERA-INTERIM, in which the WRF model was initialised with a different set of initial and boundary conditions data from ERA-INTERIM [START_REF]European Centre for Medium-Range Weather Forecasts: ERA-Interim Project[END_REF]. Note how in this case the boundary conditions can have a small impact every six hours in the simulation, not only at the starting time (boundary conditions). The rest of the model configuration of these experiments was the same than in the NEW-LC experiments; therefore, their scores were compared in Table 4 to those obtained with the NEW-LC experiment (using initial conditions from NCEP-FNL and no spin-up time).

Table 4. Summary of the scores in Wm -2 calculated in the whole evaluated area, for each LSM (columns) configured for each pre-experiment (raws, SPIN-UP and ERA-INTERIM). The scores of the NEW-LC experiments are included to quantify the impact of the change in the initial conditions. The first value refers to the sensible heat flux (SH), and the second one to the latent heat flux (Le). Some experiments were not possible (-) for some LSMs. The used scores are: the bias (random error), standard deviation (std) (systematic error) and the root mean square error (RMSE).

SPIN-UP experiment

With the SPIN-UP experiment, we have checked the LSMs sensitivity to auto-spin up their SM values performing longer simulations (1 month before the model output are used and analysed). As discussed before, in our case this was not appropriate to allow a fairer model-observation comparison (the AAF used to evaluate the model were based on the assumption of homogeneous SM, which does not exist after the SPIN-UP time). Besides, the inherent differences between the LSMs cause variations in the rainfall simulated by WRF in each LSM, which make difficult the comparison of the LSMs scores (reply AC1, 2021).

The analysis of the scores obtained from the SPIN-UP experiment in Table 4 reveals that the biases were worsened for Noah (from -17 Wm -2 to -35 Wm -2 and from 34 Wm -2 to 59 Wm -2 for SH and Le, respectively), while the std values were similar, leading to a larger RMSE for both fluxes. For Noah-MP the bias was corrected (from 66 Wm -2 to 54 Wm -2 and from -42 Wm -2

NEW-LC experiment

Two variables were modified in the geo_em_d04.nc file (the fourth-domain output from the geogrid.exe program of the WRF preprocessing system (WPS)). On the one hand, the LANDUSEF variable was modified, including the new percentages for each LC category in all the 1-km grid cells of the fourth domain. The same was done for the rest of the domains, but only in the area covered by the inner domain. On the other hand, the dominant LC category in each pixel was calculated based on the new information, which served to modify the LU_INDEX variable of the same files than before.

These files with the modified LANDUSEF and LU_INDEX variables were re-incorporated to the WPS system and the rest of the preprocessing programs were executed, i.e., ungrib.exe and metgrid.exe, obtaining the final met_em files used to run the model. Note how in order to use these modified files in the model simulations, the surface_input_source parameter in the WRF namelist.input file was set to 3.

MOSAIC experiment

To activate the mosaic approach, some options should be included in the namelist.input file of WRF. In the case of Noah, the sf_surface_mosaic option should be set to 1. Besides, the mosaic_cat option indicates the maximum number of tiles to be used, which was set to 19, the maximum possible in our area. In the case of RUC, the mosaic_lu and mosaic_soil should be included and set to 1.

FOREST experiment

The three parameters modified in the FOREST experiment were changed in the MPTABLE.TBL file of WRF. This file contains the vegetation parameters tabulated of the different LC categories for the two LC datasets available in WRF. Since our simulations used the IGBP-MODIS LC dataset, we changed these parameters in its corresponding section within the file.

Specifically, we modified MP from its original value (6) to 9, BP from 0.002 to 0.01 mol H 2 O m -2 s -1 , and V cmax25 from 50 to 62.5 µmol m -2 s -1 for the LC type 1 (ENF). This was done before running the WRF simulation with the Noah-MP LSM. Competing interests. The authors declare that they have no conflicts of interest.