Impact of Noah-LSM Parameterizations on WRF Mesoscale Simulations: Case Study of Prevailing Summer Atmospheric Conditions over a Typical Semi-Arid Region in Eastern Spain

: The current study evaluates the ability of the Weather Research and Forecasting Model (WRF) to forecast surface energy ﬂuxes over a region in Eastern Spain. Focusing on the sensitivity of the model to Land Surface Model (LSM) parameterizations, we compare the simulations provided by the original Noah LSM and the Noah LSM with multiple physics options (Noah-MP). Furthermore, we assess the WRF sensitivity to different Noah-MP physics schemes, namely the calculation of canopy stomatal resistance (OPT_CRS), the soil moisture factor for stomatal resistance (OPT_BTR), and the surface layer drag coefﬁcient (OPT_SFC). It has been found that these physics options strongly affect the energy partitioning at the land surface in short-time scale simulations. Aside from in situ observations, we use the Meteosat Second Generation (MSG) Spinning Enhanced Visible and Infrared Imager (SEVIRI) sensor to assess the Land Surface Temperature (LST) ﬁeld simulated by WRF. Regarding multiple options in Noah-MP, WRF has been conﬁgured using three distinct soil moisture factors to control stomatal resistance ( β factor) available in Noah-MP (Noah, CLM, and SSiB-types), two canopy stomatal resistance (Ball–Berry and Jarvis), and two options for surface layer drag coefﬁcients (Monin–Obukhov and Chen97 scheme). Considering the β factor schemes, CLM and SSiB-type β factors simulate very low values of the latent heat ﬂux while increasing the sensible heat ﬂux. This result has been obtained independently of the canopy stomatal resistance scheme used. Additionally, the surface skin temperature simulated by Noah-MP is colder than that obtained by the original Noah LSM. This result is also highlighted when the simulated surface skin temperature is compared to the MSG-SEVIRI LST product. The largest differences between the satellite data and the mesoscale simulations are produced using the Noah-MP conﬁgurations run with the Monin–Obukhov parameterization for surface layer drag coefﬁcients. In contrast, the Chen97 scheme shows larger surface skin temperatures than Monin–Obukhov, but at the expense of a decrease in the simulated sensible heat ﬂuxes. In this regard, the ground heat ﬂux and the net radiation play a key role in the simulation results.


Introduction
Accurate and reliable information of surface energy fluxes and atmospheric variables is of critical importance for different application fields, not only for atmospheric and climate modeling but also for agricultural, ecological, and hydrological purposes, among others. One of the efforts focused on improving the simulation of surface energy fluxes and atmospheric variables has been dedicated to an accurate simulation of land surface and atmosphere interactions [1][2][3][4][5][6][7]. All these studies have shown essential synergies between the land surface and the atmosphere. In fact, land surface has been found to be a key component of the Earth's climate system due to water and energy exchanges between these components of the Earth environment. Partition of net radiation into latent and sensible heat fluxes are determined by these exchanges [8]. Therefore, a key for this interaction between the land surface and the atmosphere is the soil-vegetation-atmosphere (SVA) feedback [1], since the parameterization of the land surface impacts a wide range of physical processes, such as cloud formation, precipitation, and evaporation.
Interactions and feedbacks between the soil, vegetation, and the atmosphere are quantified using Land Surface Models (LSMs). Parameterizations of land surface processes have evolved in the last decades based on a close multidisciplinary cooperation in the fields of hydrology, plant and soil physics, and meteorology [9]. As a result, a comprehensive representation of ecosystem and water dynamics have been incorporated based on additional feedback processes between water, terrestrial ecosystems, and climate [10]. LSMs can be applied both offline [8, [11][12][13][14][15][16][17][18] and coupled to Numerical Weather Prediction (NWP) and climate models. In this latter case, LSMs provide bottom boundary conditions for the corresponding atmospheric model by means of the calculation of energy and mass transfer at the ground interface [1,19].
The current study is focused on two LSMs available within the Weather Research and Forecasting (WRF) mesoscale numerical model [20,21], namely, the Noah scheme [9,22,23] and the Noah-MP scheme [12,13]. On the one hand, the Noah LSM is one of the most widely used land surface schemes for research studies, as well as for weather and regional climate modeling [24]. In this regard, Noah LSM has not only been coupled to the WRF model but also to the operational North American Mesoscale (NAM) Forecast System, the Global Forecast System (GFS), and other National Centers for Environmental Prediction (NCEP) modeling systems [19,25,26]. On the other hand, the Noah-MP model is an evolutional version of the original Noah LSM. Even though most of the schemes of the Noah LSM have been retained, multiparameterization options are available for selected key processes, introducing a large range of alternatives when using the corresponding LSM [12,13]. For instance, the introduction of a simple groundwater model, the dynamic treatment of the vegetation canopy, and the separation of vegetated and bare cell portions. In this regard, Noah_MP represents an augmented version of the original Noah LSM. The interaction between the land surface and the atmosphere, comparing the SVA feedbacks, are modeled by the corresponding LSM models implemented within the WRF modeling system. Both Noah and Noah-MP coupled to the Weather Research and Forecast model have been used by atmospheric and hydrological researchers for different purposes [5,6,15,16,24,[27][28][29]. In addition, these LSMs have been compared to each other for a wide range of applications. For instance, LSMs, and particularly Noah/Noah-MP LSMs, have been used to simulate energy fluxes in response to near-surface atmospheric forcing based on their representation of the land surface processes [5,6,[28][29][30][31][32][33][34]. As such, it is of significant importance to accurately document the model's performance for the simulation of surface energy fluxes [15]. Furthermore, we should also consider the interplay between surface energy fluxes and surface and near-surface magnitudes, such as temperature and moisture.
Considering all these mentioned points, the present study aims: (1) to compare the original Noah LSM performance with the augmented Noah-MP and (2) to perform sensitivity tests using different physical processes implemented in the Noah-MP model in order to assess the ability of this model to properly simulate surface energy fluxes and closely related variables. Based on these three research items, we try to identify patterns produced by distinct WRF configurations regarding Noah LSM schemes over semi-arid locations and under different dominant summer atmospheric conditions over the study area. It is important to evaluate these models over different surfaces as well as at different scales [15]. In this regard, the current work highlights the strengths and limitations of different Noah LSM configurations coupled to the WRF model under the specific conditions mentioned above. WRF is a highly used model, extensively used to test different parameterization schemes. The sensitivity experiments conducted in the current study regarding Noah and Noah-MP LSMs parameterizations are designed with the aim of improving previous errors detected using numerical mesoscale modeling over the area of study as well as determining optimal options and configurations of these LSMs that could be used as a reference in future studies.
The organization of the present paper is the following. Methodology and methods, together with the experimental design, modeling framework, and evaluation procedure, are presented in Section 2, whereas the simulation results are presented in Section 3. Finally, a discussion of these results, as well as some final conclusions and remarks, are drawn in Section 4.

Numerical Setup
The Advanced Research core of the WRF model (WRF-ARW) [20,21], version 3.6.1 (Boulder, CO, USA), was used to simulate the selected forecasting period: 6 to 11 July 2011. The model setup is based on that used in previous studies [6,28] conducted over the study area. WRF was set up using the Dudhia shortwave radiation [35] and the rapid radiative transfer model (RRTM) longwave radiation [36] as radiation options, whereas the YSU scheme [37] was used as the PBL parameterization. WRF used a nested grid configuration with an outer domain of 48 km horizontal resolution and an inner grid of 3 km, as well as a 12 km grid in between, as shown in Figure 1. Regarding the vertical structure, 45 vertical levels were used, with 24 levels below 2000 m and 9 levels below 300 m.  In relation to the period of study, a western synoptic advection dominates on 6-8 July 2011, with low atmospheric moisture and high wind speeds over the area of study. In contrast, mesoscale circulations were well established during the 9 and 10 July 2011, favoring the transition between day and nighttime winds. Finally, an eastern synoptic advection was developed on 11 July 2011. Therefore, the simulation period covers the typical summer meteorological conditions prevailing over the area of study. Independently of the dominant atmospheric situation, clear skies were observed for the complete period of study. We decided to simulate this period of study following the results obtained in previous works [2,6,28]. We selected the summer season with the aim of improving the simulation of surface energy fluxes and meteorological parameters, considering some previous errors detected using numerical mesoscale modeling over the area of study. Additionally, focusing on clear sky conditions allows the study of In relation to the period of study, a western synoptic advection dominates on 6-8 July 2011, with low atmospheric moisture and high wind speeds over the area of study. In contrast, mesoscale circulations were well established during the 9 and 10 July 2011, favoring the transition between day and nighttime winds. Finally, an eastern synoptic advection was developed on 11 July 2011. Therefore, the simulation period covers the typical summer meteorological conditions prevailing over the area of study. Independently of the dominant atmospheric situation, clear skies were observed for the complete period of study. We decided to simulate this period of study following the results obtained in previous Sustainability 2021, 13, 11399 4 of 17 works [2,6,28]. We selected the summer season with the aim of improving the simulation of surface energy fluxes and meteorological parameters, considering some previous errors detected using numerical mesoscale modeling over the area of study. Additionally, focusing on clear sky conditions allows the study of situations where the incoming solar radiation is important and to determine optimal configurations of the model that could be used as a reference in future studies and applications. Although the description of the different meteorological conditions is included here to highlight the whole picture, the analysis of results will be conducted considering the whole simulation period. For each of these days, a 36 h simulation was performed beginning the previous day at 12:00 UTC.NCEP FNL dataset (1 • × 1 • horizontal resolution and 6 h temporal resolution) was used as boundary and initial conditions. Hourly outputs are saved, leaving out the first 12 h as a spin-up time. In this regard, only the remaining 24 h are considered in the analysis of results. This approach is based on the one previously followed in [38,39].

Numerical Experiments
WRF sensitivity to LSM choices and simulated processes is investigated ( Table 1). The Noah LSM model [9,23] was used in the control WRF simulation (WRF_CON hereafter). In addition to this control simulation, six sensitivity experiments were conducted, updating the original Noah LSM with the more recent and advanced Noah LSM with multiparameterization options (Noah-MP) [12,13] (Table 1). Noah LSM was first implemented in WRF in 2004, corresponding to the early release WRF 2.0, whereas Noah-MP was first implemented in the WRF version 3.4 (Boulder, CO, USA. 2012). Even though these LSMs are not the only advanced schemes implemented in WRF, the Noah LSM has been comprehensively developed, evaluated, and applied for more than three decades, and it is probably the most widely used LSM among the WRF community; the Noah-MP model is an evolution of the original Noah LSM that includes multiparameterization and different structural changes. A full description of the Noah LSM parameterization is included in [9,23], whereas a full description of the Noah-MP scheme may be found in [12]. The Noah LSM surface layer is a combination of the soil and vegetation surfaces, where energy fluxes are resolved. This configuration does not permit explicitly calculating the canopy temperature, water, carbon and energy fluxes, and photosynthetically active radiation (PAR). Surface skin temperature plays a key role in computing the energy fluxes, which are calculated separately over the bare ground and the vegetation. As a result, the net flux considers the fluxes over the canopy and bare ground, weighted by the vegetation fraction. Noah-MP presents a framework where multiple options can be selected in order to parameterize diverse physical processes [12], such as stomatal resistance, runoff and groundwater, soil permeability, precipitation partitioning between snow and rain, as well as ground surface albedo. Augmentations to the Noah LSM (Noah-MP) include a separate computation of the canopy temperature and the ground surface temperature. This is accomplished through a canopy layer for vegetation. Additionally, in order to calculate shaded and sunlit leaves fractions and related parameters, a modified two-stream scheme has been incorporated.
Three distinct physical processes implemented in Noah-MP are evaluated in the current study: the soil moisture threshold for transpiration (OPT_BTR), canopy stomatal resistance (OPT_CRS), and surface layer exchange (OPT_SFC). These three surface physical processes coupled in Noah-MP were found to be the most sensitive in the study conducted by [18], and in this order starting from the highest. That is, the soil moisture factor for stomatal resistance (OPT_BTR) has been shown to be the most sensible physical process in Noah-MP. Additionally, these process parameterizations were chosen based on their significant importance [13,17,[40][41][42][43]. We adopted the other multiparameterization options as the default option sets suggested for Noah-MP. WRF_MP1 in the current study is the same as the EXP4 configuration in Niu et al. (2011) [11] and Yang et al. (2011) [13], whereas WRF_MP4 is the same as EXP5 in these previous studies, as well as [29,30]. WRF_MP1 has also been used by different authors [15,19,44,45]. Furthermore, WRF_MP6 is similar to the configuration used by Milovac et al. (2016) [1], in terms of OPT_CRS, OPT_CRS and OPT_SFC, while Pilotto et al. (2015) [8] used a similar configuration to WRF_MP1 but with the scheme used in earlier versions of Noah [45] for the surface exchange coefficients.
The first experiment performed in the current work using the Noah-MP LSM (Table 1) WRF_MP1 is based on the combination of the following schemes: a canopy stomatal resistance following the Jarvis scheme and a soil moisture factor for stomatal resistance (also known as the β factor) based on the Noah scheme, TOPMODEL with groundwater for groundwater and runoff, and the Monin-Obukhov (M-O) parameterization for the surface layer drag coefficient calculation. TOPMODEL simulates downslope lateral transport of water both on the surface and subsurface within saturated regions of the soil, that is, a groundwater redistribution provoked by the gravity-driven downslope flow. Two additional sensitivity experiments, WRF_MP2 and WRF_MP3, are focused on evaluating the impact of soil moisture limitation on transpiration in the WRF model results (β factor). The soil moisture threshold for transpiration determines the efficiency of plant transpiration for a given soil moisture amount. Therefore, the β factor, also known as the drought stress factor, is critical for terrestrial ecosystem dynamics and its interactions with climatic and hydrological processes [11]. There are three options for this soil-water stress function in Noah-MP (OPT_BTR): Noah LMS's version (WRF_MP1), the Community Land Model (CLM), and the Simplified Simple Biosphere (SSiB) β factor. The Noah-type β (WRF_MP1) factor is essentially a liquid volumetric capacity model, calculated as a function of the soil moisture at wilting point of the corresponding soil layer and reference soil moisture parameters [21]. WRF_MP2 and WRF_MP3 were set up similarly to the WRF_MP1 configuration, but with the CLM and SSiB-type factor scheme instead of the Noah-type factor scheme, respectively. CLM (WRF_MP2) is defined as a function of the soil layer matric potential, the wilting matric potential and the saturated matric potential. SSiB (WRF_MP3) computes the β factor using the wilting matric potential and the soil layer matric potential. Both β-type approaches use liquid soil moisture to calculate water potential. More information about the CLM-type β factor may be consulted in [46], whereas the SSiB-type β factor scheme may be consulted in [47]. Noah-type β factor soil moisture depends on soil type while SSiB includes, in the OPT_BTR formulation, a slope factor that depends on the vegetation type. In the study conducted by [18], they found that Noah-type β factor parameterization was their best option for the corresponding physical process in comparison to CLM and SSiB parameterizations. In this regard, the Noah β-type showed larger variations with soil moisture than the CLM and SSiB β-types [12], due to the soil moisture overestimation produced by CLM and SSiB. A better simulation of latent heat fluxes (LH) and sensible heat fluxes (SH) is produced by Noah β-type because higher transpiration values implemented in the corresponding schemes could result in the over evaluation of LH [18].
Another Noah-MP physical parameterization scheme investigated in the current study is the canopy stomatal resistance (OPT_CRS). Two options of this parameter are available in Noah-MP, namely Jarvis [48] and Ball_Berry [49]. Regarding these options, the Jarvis scheme proposes a parameterization of the canopy resistance based on different parameters: the minimal stomatal resistance, the Leaf Areas Index (LAI), and different stress functions that account for the effects of air temperature, vapor pressure deficit, soil moisture stress, and solar radiation, [24,48]. In contrast, the Ball-Berry type scheme links stomatal resistance to rates of photosynthesis per unit LAI of sunlit and shaded leaves [49]. In the current study, two sensitivity experiments were performed using the Ball-Berry canopy resistance parameterization. On the one hand, this scheme is selected together with the Noah βtype option (WRF_MP4) and the CLM β-type option (WRF_MP5). Finally, an additional experiment (WRF_MP6) was conducted using the Chen97 option for exchange in the surface layer (OPT_SFC). Noah-MP implements two different computations of the surface exchange coefficient of water and heat. The first one is that used in Noah version 3.0 (Boulder, CO, USA) [45]. This scheme uses a different roughness length for momentum and heat. Another option implemented in Noah-MP for the parameterizations of heat surface exchange coefficient uses the M-O similarity theory [44], where the same roughness length for momentum and heat is used, but a zero-displacement height is included when computing the surface exchange coefficient for heat. Chen97 and the more general M-O scheme take the same stability correction functions under unstable and stable conditions. However, the zero-displacement height used in [44] is not considered in Noah version 3.0 [45].

Observational Datasets
WRF-simulated LH and SH were assessed using observational datasets of these magnitudes over an anchor FLUXNET station ( Figure 1) located over El Bonillo (BON). The FLUXNET observations over BON are sampled every ten minutes. However, instantaneous hourly data are used in the current study to match the WRF time output resolution. Two-meter temperature and relative humidity measured over this weather station are used, as well in the model assessment. The BON site covers 12,872 ha. This area is one of the main semi-arid juniper woodlands in Spain. J. thurifera is found in the more unfavorable soils, whereas Quercus ilex is found in the best soils. Lithic leptosol is the dominant soil type, according to FAO (1988). The USGS (United States Geological Survey) 30-second global 24-category vegetation (land-use) data is used in the WRF simulations. In this regard, the dominant land-use category within the WRF grid cell where the simulation results are compared to the observations over BON corresponds to the Cropland/Woodland Mosaic.
Additionally, WRF-simulated LH and SH are also presented over the portable weather station located in Barrax (BRX; Figure 1). The soil in the BRX study area has a silty-clay loam texture (13.4% sand, 48.9% silt, and 37.7% clay) and is cataloged as Petrocalcic Calcixerepts. Additionally, this is a flat agricultural plot with different croplands. The dominant land-use category used in WRF for the corresponding BRX grid cell is the same as that over BON (Cropland/Woodland Mosaic).
Even though no LH and SH observations are available over this latter weather station, measures of 2-m relative humidity and temperature are used as a reference as well. Based on these magnitudes, the specific humidity is computed as well. Additionally, measures of shortwave and longwave downward radiation, as well as net radiation recorded over BON and BRX, are also used in the current study. BON observed net radiation together with LH and SH are used to compute the ground energy flux over this anchor FLUXNET station. Finally, besides in situ observations, Meteosat Second Generation (MSG) Spinning Enhanced Visible and Infrared Imager (SEVIRI) (MSG-SEVIRI) Land Surface Temperatures (LSTs) are also used [50,51]. MSG-SEVIRI is characterized by an approximate 3 × 4.5 km footprint of the region. The MSG-SEVIRI LST product used in the current study presents a 15 min temporal resolution and is generated by the EUMETSAT LSA SAF (available at https://landsaf.ipma.pt/en/, accessed on 10 October 2021). With the aim of extracting the corresponding LST for BON and BRX locations, the spatial "nearest point" interpolation method is used to compare the MSG-SEVIRI LST with the simulation outcomes. The model results are validated over BON and BRX using different statistical metrics, taking into account the whole simulation period (Tables 2-4). These statistical measures include the Mean Bias Error (MBE) and the Root Mean Square Error (RMSE), as defined by the following equations:   Additionally, we use the Correlation Coefficient (R) and the Index of Agreement (IoA) in the evaluation procedure: where F and O represent time average simulated and observed values, respectively, while F represents the simulated value, O the observation value, and N the number of observations included in the calculation [52,53]. MBE is defined as the average of the simulated value minus the observed value and quantifies the systematic error of the model, while RMSE is the square root of the individual differences between the simulated and observed values; it quantifies the accuracy of the model. R measures the correlation between forecasts and observations. The model correlates with observations if the value of R is near 1. Finally, the IoA is a modified correlation coefficient that measures the degree to which a model's prediction is free of error. A value of 0 means complete disagreement, whereas a value of 1 implies a perfect score. Figure 2 shows the SH and LH results over BON and BRX, whereas the statistical scores for these parameters over the BON FLUXNET station are presented in Table 2. WRF_CON produces the lowest MBE scores for LH over BON, with a value of −5 W m −2 , an RMSE of 18 W m −2 and a correlation coefficient and IoA of 0.782 and 0.821, respectively ( Table 2). WRF_MP1 produces an overestimation of the observed LH, as shown in Figure 2c. In general, WRF_CON tends to underestimate the observations of this magnitude, while WRF_MP1 shows the opposite trend. WRF_MP2 shows a clear overestimation of the observed SH, as shown by an MBE of 50 W m −2 ( Table 2). These results obtained for SH are related to a significant underestimation of the observed LH (Figure 2c), with an MBE of −20 W m −2 . In addition, the highest LH simulated by WRF_MP2 over BON is around 6 W m −2 .

Observations Versus Model Comparison
Sustainability 2021, 13, x FOR PEER REVIEW 9 of 19 respectively ( Table 2). WRF_MP1 produces an overestimation of the observed LH, as shown in Figure 2c. In general, WRF_CON tends to underestimate the observations of this magnitude, while WRF_MP1 shows the opposite trend. WRF_MP2 shows a clear overestimation of the observed SH, as shown by an MBE of 50 W m −2 ( Table 2). These results obtained for SH are related to a significant underestimation of the observed LH ( Figure 2c), with an MBE of −20 W m −2 . In addition, the highest LH simulated by WRF_MP2 over BON is around 6 W m −2 . Confronting WRF_MP3 with WRF_MP2, no significant differences are obtained for LH ( Table 2). The negative MBE of −20 W m −2 produced by both WRF_MP2 and WRF_MP3 is related to a more positive SH MBE (50 W m −2 ). This MBE score is larger than those found using other WRF-MP configurations, such as WRF_MP1 or WRF_MP4, but similar to that produced by WRF_MP5 ( Table 2). The β factor represents the effect of soil moisture limitation on surface humidity and controls the evapotranspiration term. In this sense, Noah-type (WRF_MP1) increases the transpiration rate, leading to a moister environment than that simulated by the CLM-type formulation (WRF_MP2).
Additionally, Figure 2c also highlights a relevant similarity between WRF_CON and WRF_MP6 in terms of LH, which is also reflected in Table 2 by means of similar statistical scores. In terms of SH, however, WRF_MP6 tends to underestimate the observed field, with an MBE of −12 W m −2 , instead of the positive MBE of 14 W m −2 found in the WRF_CON simulation (Table 2). WRF_MP6 produces the lowest SH (Figure 2a), while still reproducing a similar LH pattern to the one produced by WRF_CON.
Focusing on the BRX weather station (Figure 2b,d), the WRF_CON and WRF_MP Confronting WRF_MP3 with WRF_MP2, no significant differences are obtained for LH ( Table 2). The negative MBE of −20 W m −2 produced by both WRF_MP2 and WRF_MP3 is related to a more positive SH MBE (50 W m −2 ). This MBE score is larger than those found using other WRF-MP configurations, such as WRF_MP1 or WRF_MP4, but similar to that produced by WRF_MP5 ( Table 2). The β factor represents the effect of soil moisture limitation on surface humidity and controls the evapotranspiration term. In this sense, Noah-type (WRF_MP1) increases the transpiration rate, leading to a moister environment than that simulated by the CLM-type formulation (WRF_MP2).
Additionally, Figure 2c also highlights a relevant similarity between WRF_CON and WRF_MP6 in terms of LH, which is also reflected in Table 2 by means of similar statistical scores. In terms of SH, however, WRF_MP6 tends to underestimate the observed field, with an MBE of −12 W m −2 , instead of the positive MBE of 14 W m −2 found in the WRF_CON simulation (Table 2). WRF_MP6 produces the lowest SH (Figure 2a), while still reproducing a similar LH pattern to the one produced by WRF_CON.
Focusing on the BRX weather station (Figure 2b,d), the WRF_CON and WRF_MP outcomes show similar patterns for both SH and LH to those simulated over BON.
A characteristic that appears comparing WRF_MP1 and WRF_MP2 (similar for WRF_MP3; not shown) is a difference of around 50 W m −2 for SH and around 100 W m −2 for LH, producing a significant reduction of the simulated LH over the two weather stations. Moreover, the maximum LH is produced around 12 UTC by WRF_MP1, while the maximum LH simulated by WRF_MP2 is simulated at 14 UTC, thus delaying this magnitude for about 2 h (Figure 2c,d).
Considering the net radiation (Figure 3a,b), similar results are found between the distinct WRF configurations over BON and BRX (Tables 3 and 4  Moreover, the maximum LH is produced around 12 UTC by WRF_MP1, while the maximum LH simulated by WRF_MP2 is simulated at 14 UTC, thus delaying this magnitude for about 2 hours (Figure 2c,d).
Considering the net radiation (Figure 3a,b), similar results are found between the distinct WRF configurations over BON and BRX (Table 3 and Table 4  Based on these results, WRF_MP6 provides the lowest net radiation to be partitioned into sensible, latent, and ground heat fluxes than the other WRF-MP experiments. WRF_MP6 produces around 50 and 350 W m −2 for LH and SH, respectively ( Figure 2). Therefore, plausible reasons for the differences found between WRF_MP6 and the other WRF-MP simulations could be related to both the available net radiation and the part of this magnitude partitioned into ground heat flux. In this regard, the main difference in the selection of physical options is the application of the surface exchange Based on these results, WRF_MP6 provides the lowest net radiation to be partitioned into sensible, latent, and ground heat fluxes than the other WRF-MP experiments. WRF_MP6 produces around 50 and 350 W m −2 for LH and SH, respectively ( Figure 2). Therefore, plausible reasons for the differences found between WRF_MP6 and the other WRF-MP simulations could be related to both the available net radiation and the part of this magnitude partitioned into ground heat flux. In this regard, the main difference in the selection of physical options is the application of the surface exchange coefficient for heat proposed by Chen et al., (1997) [45] in WRF_MP6 (the Chen97 scheme). In contrast, the other WRF-MP simulations use the more general M-O similarity theory [44].
No significant differences arise comparing the shortwave and longwave downward radiation simulated by WRF experiments (Figure 4), as shown by the Mann-Whitney-Wilcoxon test (p-value > 0.5). An MBE of 4-5 W m −2 with a RMSE of 40 W m −2 is simulated by WRF_CON and WRF-MP for the shortwave downward radiation over BON (Table 3). In addition, an MBE around −10 W m −2 , as well as a RMSE of 13-14 W m −2 , are simulated by WRF_CON and WRF-MP for the longwave downward radiation over BON. Higher errors of both the shortwave and longwave downward radiation are simulated over BRX (Table 4). Still, the model tends to underestimate the observations of these two magnitudes both over BON and BRX. Considering these results, the simulated differences for the surface energy fluxes among the WRF configurations do not seem to be related to these radiation magnitudes. scheme). In contrast, the other WRF-MP simulations use the more general M-O similarity theory [44]. No significant differences arise comparing the shortwave and longwave downward radiation simulated by WRF experiments (Figure 4), as shown by the Mann-Whitney-Wilcoxon test (p-value > 0.5). An MBE of 4-5 W m −2 with a RMSE of 40 W m −2 is simulated by WRF_CON and WRF-MP for the shortwave downward radiation over BON (Table 3). In addition, an MBE around −10 W m −2 , as well as a RMSE of 13-14 W m −2 , are simulated by WRF_CON and WRF-MP for the longwave downward radiation over BON. Higher errors of both the shortwave and longwave downward radiation are simulated over BRX (Table 4). Still, the model tends to underestimate the observations of these two magnitudes both over BON and BRX. Considering these results, the simulated differences for the surface energy fluxes among the WRF configurations do not seem to be related to these radiation magnitudes.  Figure 5 shows similar 2-m temperatures for WRF_CON and WRF_MP1, although slightly lower in the latter case. This result is related to the moister field presented in Figure 6 in terms of both the specific and relative humidities. However, notable differences, up to 10 °C, are obtained between WRF_CON and WRF_MP1 considering the skin temperature (Figure 5c,d). Regarding the LST parameter, the current study shows significant differences between WRF_CON, WRF_MP6, and WRF_MP1-5. In this regard, WRF_CON really captures the evolution of the SEVIRI LST product over BRX, whereas WRF_MP6 matches the LST field provided by SEVIRI over BON. The other WRF_MP simulations (1)(2)(3)(4)(5) show a clear tendency to underestimate the LST parameter derived from the SEVIRI instrument. If we focus on the ground energy flux simulated by WRF_CON and WRF_MP1 (Figure 3c,d), WRF produces ground values around 20-30 W m −2 larger than WRF_MP1. The increased skin temperature simulated by WRF_CON compared to that simulated by WRF-MP seems to be related to this larger ground heat flux produced by WRF_CON (Figure 3c,d). Additionally, there is a shift in the hour where the maximum ground flux is reached. In this regard, WRF_CON presents the  Figure 5 shows similar 2-m temperatures for WRF_CON and WRF_MP1, although slightly lower in the latter case. This result is related to the moister field presented in Figure 6 in terms of both the specific and relative humidities. However, notable differences, up to 10 • C, are obtained between WRF_CON and WRF_MP1 considering the skin temperature (Figure 5c,d). Regarding the LST parameter, the current study shows significant differences between WRF_CON, WRF_MP6, and WRF_MP1-5. In this regard, WRF_CON really captures the evolution of the SEVIRI LST product over BRX, whereas WRF_MP6 matches the LST field provided by SEVIRI over BON. The other WRF_MP simulations (1)(2)(3)(4)(5) show a clear tendency to underestimate the LST parameter derived from the SEVIRI instrument. If we focus on the ground energy flux simulated by WRF_CON and WRF_MP1 (Figure 3c,d), WRF produces ground values around 20-30 W m −2 larger than WRF_MP1. The increased skin temperature simulated by WRF_CON compared to that simulated by WRF-MP seems to be related to this larger ground heat flux produced by WRF_CON (Figure 3c,d). Additionally, there is a shift in the hour where the maximum ground flux is reached. In this regard, WRF_CON presents the steepest curve for this magnitude, reaching its maximum peak earlier in time and with larger values than WRF-MP. These issues could be related to the different schemes used for the surface exchange coefficient for heat. While the M-O parameterization for the surface exchange coefficient is used by WRF-MP (with the exception of WRF_MP6), WRF_CON used the original parameterization implemented in Noah (Chen97). According to Niu et al., (2011) [12], Noah-MP produced a larger diurnal amplitude than Noah for the ground heat flux, better capturing the observations. The results found in [12] contrast with those found here over the region of study comparing Noah-MP (WRF_MP1) and Noah (WRF_CON). Nevertheless, Niu et al., (2011) [12] highlighted that "the better simulations may not hold for all combinations of uncertain parameters" and that evaluation of models' performance should continue in the future, as pointed out by Gulden et al., (2008) [54], and it is performed in the current study.  [12], Noah-MP produced a larger diurnal amplitude than Noah for the ground heat flux, better capturing the observations. The results found in [12] contrast with those found here over the region of study comparing Noah-MP (WRF_MP1) and Noah (WRF_CON). Nevertheless, Niu et al. (2011) [12] highlighted that "the better simulations may not hold for all combinations of uncertain parameters" and that evaluation of models' performance should continue in the future, as pointed out by Gulden et al. (2008) [54], and it is performed in the current study.  The relation between the ground heat flux and the skin temperature is observed, as well, comparing the WRF_CON and WRF_MP6 simulations. Larger ground fluxes, simulated by WRF_MP6 (Figure 3c,d), produce higher skin temperatures (Figure 4c,d) than WRF_CON, both over BON and BRX. Furthermore, WRF_MP6 significantly increases the simulated 2-m temperatures compared to the other six experiments (Figure 5a,b), while WRF_CON produces similar 2-m temperatures to WRF-MP, running with the M-O surface layer drag coefficient parameterization. The WRF_MP6 2-m temperature overestimation may be caused by errors in the initial soil moisture. However, no information on this parameter is available over the area of study and the simulation period. Therefore, and according to the results found using WRF_MP6, the corresponding daytime 2-m temperature biases obtained here should be investigated in the future, including suitable soil moisture measurements. Higher differences between skin temperatures and 2-m temperatures, such as in the case of WRF_CON or WRF_MP6, should be related to larger SH values. However, this is not the case here. In contrast, WRF_MP6 produces the lowest SH (Figure 2a,b) compared to the other WRF-MP experiments, with differences up to 200 W m −2 . Therefore, there should be another reason for these low SHs. Figure 7 shows the surface exchange coefficient simulated by WRF-MP. Even though similar values are obtained for this magnitude for all WRF-MP simulations, WRF_MP6 is an exception that leads to significantly low C h values. This seems to be the most plausible reason for the lowest SH simulated by WRF_MP6 in comparison to those simulated by the other WRF-MP configurations. In the current study, the M-O scheme produces cold biases in terms of the skin temperature (Figure 5c,d). In contrast, the Chen97 scheme (WRF and WRF_MP6) improves the simulation of this magnitude over the area and period of study by means of lower surface exchange coefficients (Figure 7), which reduces ventilation and increases the land surface heating during the daytime summer period. The relation between the ground heat flux and the skin temperature is observed, as well, comparing the WRF_CON and WRF_MP6 simulations. Larger ground fluxes, simulated by WRF_MP6 (Figure 3c,d), produce higher skin temperatures (Figure 4c,d) than WRF_CON, both over BON and BRX. Furthermore, WRF_MP6 significantly increases the simulated 2-m temperatures compared to the other six experiments ( Figure  5a,b), while WRF_CON produces similar 2-m temperatures to WRF-MP, running with the M-O surface layer drag coefficient parameterization. The WRF_MP6 2-m temperature overestimation may be caused by errors in the initial soil moisture. However, no information on this parameter is available over the area of study and the simulation period. Therefore, and according to the results found using WRF_MP6, the corresponding daytime 2-m temperature biases obtained here should be investigated in the future, including suitable soil moisture measurements. Higher differences between skin temperatures and 2-m temperatures, such as in the case of WRF_CON or WRF_MP6, should be related to larger SH values. However, this is not the case here. In contrast, WRF_MP6 produces the lowest SH (Figure 2a,b) compared to the other WRF-MP experiments, with differences up to 200 W m −2 . Therefore, there should be another reason for these low SHs. Figure 7 shows the surface exchange coefficient simulated by WRF-MP. Even though similar values are obtained for this magnitude for all WRF-MP simulations, WRF_MP6 is an exception that leads to significantly low Ch values. This seems to be the most plausible reason for the lowest SH simulated by WRF_MP6 in comparison to those simulated by the other WRF-MP configurations. In the current study, the M-O scheme produces cold biases in terms of the skin temperature ( Figure  5c,d). In contrast, the Chen97 scheme (WRF and WRF_MP6) improves the simulation of this magnitude over the area and period of study by means of lower surface exchange coefficients (Figure 7), which reduces ventilation and increases the land surface heating during the daytime summer period.

Different Additional Parameters Investigated
Looking for a plausible explanation regarding the differences found between WRF_MP1-5, we have also investigated different additional parameters. It seems that WRF_MP2 and WRF_MP3 find more resistance to the release of LH than WRF_MP1. The range of variation of LH simulated by WRF-MP (Figure 2c,d) shows decreasing values starting from WRF_MP1, WRF_MP4, WRF_MP5, and WRF_MP2 (as well as WRF_MP3). The transpiration field shows a similar distribution to that found for LH (not shown). Focusing on the soil moisture content (Figure 8), the reversed trend as that found for LH and transpiration is produced by WRF_MP. Therefore, the lowest LH is related to the highest soil moisture content, as expected. The Ball-Berry scheme links sunlit and shaded leaves stomatal resistance to photosynthesis [49]. Using this physical parameter, Noah-MP takes into account the stomatal resistance per unit LAI of shaded and sunlit leaves [12]. Both resistances are included in the parameterization of LH due to vegetation (vegetation canopy). The values of both stomatal resistances are significantly greater for WRF_MP2 and WRF_MP3 than for the other WRF-MP simulations (Figure 9). Therefore, WRF_MP2 and WRF_MP3 produce, in fact, more resistance to the release of LH, which is compatible with larger soil moisture contents (Figure 8).

Different Additional Parameters Investigated
Looking for a plausible explanation regarding the differences found between WRF_MP1-5, we have also investigated different additional parameters. It seems that WRF_MP2 and WRF_MP3 find more resistance to the release of LH than WRF_MP1. The range of variation of LH simulated by WRF-MP (Figure 2c,d) shows decreasing values starting from WRF_MP1, WRF_MP4, WRF_MP5, and WRF_MP2 (as well as WRF_MP3). The transpiration field shows a similar distribution to that found for LH (not shown). Focusing on the soil moisture content (Figure 8), the reversed trend as that found for LH and transpiration is produced by WRF_MP. Therefore, the lowest LH is related to the highest soil moisture content, as expected. The Ball-Berry scheme links sunlit and shaded leaves stomatal resistance to photosynthesis [49]. Using this physical parameter, Noah-MP takes into account the stomatal resistance per unit LAI of shaded and sunlit leaves [12]. Both resistances are included in the parameterization of LH due to vegetation (vegetation . The values of both stomatal resistances are significantly greater for WRF_MP2 and WRF_MP3 than for the other WRF-MP simulations ( Figure 9). Therefore, WRF_MP2 and WRF_MP3 produce, in fact, more resistance to the release of LH, which is compatible with larger soil moisture contents (Figure 8).
starting from WRF_MP1, WRF_MP4, WRF_MP5, and WRF_MP2 (as well as WRF_MP3). The transpiration field shows a similar distribution to that found for LH (not shown). Focusing on the soil moisture content (Figure 8), the reversed trend as that found for LH and transpiration is produced by WRF_MP. Therefore, the lowest LH is related to the highest soil moisture content, as expected. The Ball-Berry scheme links sunlit and shaded leaves stomatal resistance to photosynthesis [49]. Using this physical parameter, Noah-MP takes into account the stomatal resistance per unit LAI of shaded and sunlit leaves [12]. Both resistances are included in the parameterization of LH due to vegetation (vegetation canopy). The values of both stomatal resistances are significantly greater for WRF_MP2 and WRF_MP3 than for the other WRF-MP simulations ( Figure 9). Therefore, WRF_MP2 and WRF_MP3 produce, in fact, more resistance to the release of LH, which is compatible with larger soil moisture contents (Figure 8).  Still, the lowest soil moisture content is simulated by WRF_MP5. However, according to LH and the transpiration fluxes, it should be expected that this simulation produced a soil moisture content in between WRF_MP2 and WRF_MP4. Examining the stomatal resistance to photosynthesis of sunlit and shaded leaves, WRF_MP5 produces greater resistances than WRF_MP1 and WRF_MP4 during the daytime. Figure 9 shows that WRF_MP1 and WRF_MP4 configurations produce a similar pattern in terms of stomatal resistance. The highest resistance simulated by WRF_MP5 in comparison to WRF_MP1 and WRF_MP4 leads to lower LH and transpiration.

Summary and Discussion
Surface exchange coefficients play a key role in controlling the land surface temperature. In fact, the experiments conducted by Yang et al. (2011) [13] indicated that the surface exchange coefficient is a key parameter for modeling LST. In this regard, larger surface exchange coefficients, such as those obtained using the M-O scheme, mean more effective land surface ventilation and larger cooling during the summer season daytime [12]. A previous study [7], conducted over the same area of study and the same simulation period as in the current one, supports the results found here. This previous work showed that the M-O scheme has the potential to additionally improve the simulation when parameterizing a different roughness length for momentum and heat, similarly to the Chen97 scheme. Using a different treatment for the roughness length for momentum and heat corrected the original cold skin temperature bias produced by the Still, the lowest soil moisture content is simulated by WRF_MP5. However, according to LH and the transpiration fluxes, it should be expected that this simulation produced a soil moisture content in between WRF_MP2 and WRF_MP4. Examining the stomatal resistance to photosynthesis of sunlit and shaded leaves, WRF_MP5 produces greater resistances than WRF_MP1 and WRF_MP4 during the daytime. Figure 9 shows that WRF_MP1 and WRF_MP4 configurations produce a similar pattern in terms of stomatal resistance. The highest resistance simulated by WRF_MP5 in comparison to WRF_MP1 and WRF_MP4 leads to lower LH and transpiration.

Summary and Discussion
Surface exchange coefficients play a key role in controlling the land surface temperature. In fact, the experiments conducted by Yang et al. (2011) [13] indicated that the surface exchange coefficient is a key parameter for modeling LST. In this regard, larger surface exchange coefficients, such as those obtained using the M-O scheme, mean more effective land surface ventilation and larger cooling during the summer season daytime [12]. A previous study [7], conducted over the same area of study and the same simulation period as in the current one, supports the results found here. This previous work showed that the M-O scheme has the potential to additionally improve the simulation when parameterizing a different roughness length for momentum and heat, similarly to the Chen97 scheme. Using a different treatment for the roughness length for momentum and heat corrected the original cold skin temperature bias produced by the M-O scheme over the area of study and produced similar temperatures to that provided by the SEVIRI LST product. These results were achieved by means of the significantly reduced surface exchange coefficient compared to that simulated by the M-O scheme using the same roughness length for heat and momentum.
Soil moisture has an effect on the volumetric heat capacity, as well as the surface soil thermal conductivity. It appears that the Noah-MP β factor parameterized with the CLM-type formulation produces a drier environment than that simulated by the Noahtype scheme. This statement agrees with the results presented by Niu et al., (2011) [11], comparing the β factor produced by distinct formulations as a function of soil moisture, where CLM-type β factor presents a sharper and narrower variation with soil moisture than the Noah-type.
The results found in the current study in terms of Noah and Noah-MP suggest that Noah-MP and, therefore, this model coupled to WRF, are very sensitive to distinct physics processes, such as the surface layer drag coefficient (OPT_SFC), the soil moisture factor for stomatal resistance (OPT_BTR), and the calculation of canopy stomatal resistance (OPT_CRS). Considering the results obtained regarding surface energy fluxes, the original Noah produces good skill scores considering the area and period of study. Likewise, Noah-MP configured using the Chen97 surface layer drag coefficient scheme yields similar LH results to those obtained by Noah, but at the expense of a decrease in the SH. Nevertheless, Noah-MP configured using the Chen97 surface layer drag coefficient is a skillful configuration to simulate surface energy fluxes, both LH and SH, producing in general better SH results than the original Noah LSM over BON. However, both the original Noah and Noah-MP configured using the Chen97 surface layer drag coefficient underestimated LH observations, which is translated in an underestimation of the observed surface energy flux. The current study has shown that adjusting model parameters or changing schemes for different processes may correct the mismatch between the model and the observations [12]. Considering the different combinations of the canopy stomatal resistance (OPT_CRS) and the soil moisture factor controlling the stomatal resistance (OPT_BTR) parameterizations, the following results arise. On the one hand, the combination of the Jarvis scheme for the stomatal resistance and the Noah β factor type using soil moisture (WRF_MP1) produces a similar SH to the original Noah scheme (WRF_CON) while producing the largest simulated LH. On the other hand, the combination of the Jarvis scheme for the stomatal resistance and the CLM (WRF_MP2) or the SSiB β (WRF_MP3) factor using matric potential leads to larger SH and to an important decrease in the simulated LH. Moreover, the combination of the Noah β factor and the Ball-Berry scheme for stomatal resistance (WRF_MP4) leads to slightly higher SH and lower LH than the combination of the Jarvis scheme and the Noah β factor (WRF_MP1). However, the combination of the CLM β factor and the Ball-Berry parameterization (WRF_MP5) leads to slightly lower SH and higher LH than the Jarvis stomatal resistance and the CLM β factor type using matric potential (WRF_MP2). Focusing only on the β factor controlling the stomatal resistance, the Noah-type produces lower SH and higher LH than the CLM-type scheme, independently of the stomatal resistance parameterization used. Finally, fixing the β factor to the CLM type, the Ball-Berry (WRF_MP5) scheme for stomatal resistance increases LH in comparison to the results obtained using the Jarvis scheme (WRF_MP2) while slightly reducing SH. In contrast, fixing the β factor to the Noah type, the Ball-Berry scheme for stomatal resistance (WRF_MP4) decreases LH simulated by the Jarvis scheme (WRF_MP1) while slightly raising SH. Therefore, Ball-Berry coupled to the Noah β factor type shows just the opposite trend than Ball-Berry coupled to the CLM β factor type in relation to the results found using Jarvis stomatal resistance. Additionally, WRF_MP2 and WRF_MP3 produce the largest MBE and RMSE statistics in terms of both SH and LH over the area of study and the period of study.
Based on the differences found in the current study among the experiments performed using the WRF modeling framework, we have tried to identify plausible explanations for these differences. Stomatal resistance to photosynthesis of sunlit and shaded leaves has provided some insight in this sense.
To conclude, we must highlight that the results obtained here in relation to landatmosphere coupling have shown that both WRF_CON and WRF_MP6 configurations could be very useful for modeling surface energy fluxes considering the conditions tested. Even though further experiments should be performed considering longer time periods, the outcomes of the current study encourage us to continue testing LSM parameterizations bearing in mind the production of valuable operational meteorological parameters where modeling surface-atmosphere exchanges play an essential role. In this regard, the sensitivity experiments conducted in the present study have shown their significant importance and helpfulness in order to determine optimal configurations of a specific model, in this case, the WRF modeling environment. Therefore, the outcomes found in this work may serve as a reference for future studies and different applications using the WRF model. Funding: This research was funded by the Assistance Programmes of the University of Alicante "Programa de Redes-I3CE de calidad, innovación e investigación en docencia universitaria. Convocatoria 2018-2019. Alicante: Instituto de Ciencias de la Educación (ICE) de la Universidad de Alicante. Ref: [4334]." and "Programa de Redes-I3CE de calidad, innovación e investigación en docencia universitaria. Convocatoria 2020-21. Alicante: Instituto de Ciencias de la Educación (ICE) de la Universidad de Alicante. Ref: [5150]." as well as by Research Group VIGROB-116 (University of Alicante) and by the Spanish Ministerio de Ciencia e Innovación through the project PID2020-118797RB-I00/AEI.