The Role of Parameterizations and Model Coupling on Simulations of Energy and Water Balances – Investigations With the Atmospheric Model WRF and the Hydrologic Model WRF‐Hydro

The distributed hydrologic model WRF‐Hydro can operate in a fully‐coupled mode with the atmospheric Weather, Research and Forecasting (WRF) model. WRF‐Hydro enhances the modeling of terrestrial hydrologic processes in coupled WRF/WRF‐Hydro by simulating lateral surface and subsurface water flows. The objectives of this study are (a) to examine the effect of WRF‐Hydro on the surface energy and water balance in fully‐coupled WRF/WRF‐Hydro simulations and (b) to examine the impact of five WRF physics parameterizations on WRF‐Hydro streamflow. The study area is the Mediterranean island of Cyprus and 31 mountainous watersheds. The domain‐average soil moisture was 20% higher in the coupled WRF/WRF‐Hydro, relative to the standalone WRF model, during a 1‐year simulation. The higher soil moisture could explain the increase in latent heat (36%) and evapotranspiration (33%). The increase in these fluxes was less with a modification in the model transpiration parameterization to represent nocturnal transpiration and the use of remote sensing leaf area index data. The simulated precipitation of the coupled model increased up to 3%, relative to WRF. Two‐year long WRF‐Hydro simulations gave a median Nash‐Sutcliffe Efficiency for daily streamflow of the 31 watersheds of 0.5 for observed precipitation forcing and between −1.9 and 0.2 for the forcing of the five WRF parameterizations. This study showed that the enhancement of the standalone WRF model with lateral water flow processes in the coupled mode with WRF‐Hydro modifies the terrestrial energy and water balance. The improved terrestrial process representation should be considered for future hydrological cycle studies with WRF.


Introduction
Land surface models (LSM) describe land surface processes and regulate hydrological fluxes between the terrestrial and the atmospheric systems in atmospheric and climate models.However, these fluxes are only represented in the vertical direction.The dynamical coupling of coupled atmospheric -hydrological models was presented in order to simulate surface water, soil water and groundwater in the lateral direction as well (Maxwell et al., 2007(Maxwell et al., , 2011)).The simulation of lateral flows in coupled models allows for a horizontal redistribution of soil water, potentially modifying and improving the simulation of components of the terrestrial water balance and of the combined terrestrial -atmospheric water cycle (Arnault et al., 2021;Larsen et al., 2014;Shrestha et al., 2014;Yeh & Eltahir, 2005;York et al., 2002).The accomplishment of the simulation of the entire terrestrial -atmospheric water cycle is accompanied, however, by uncertainties in the two components of the coupled model.
Ensemble modeling techniques are used to assess uncertainties in the simulation of the hydrological and atmospheric models (Arsenault et al., 2018;Rummler et al., 2019).At the same time, the use of ground-based or remotely-sensed observations such as leaf area index, streamflow and evapotranspiration fluxes to parameterize and calibrate the LSM and hydrological model has the potential to improve the representation of terrestrial processes in coupled simulations (Herman et al., 2018;Sofokleous et al., 2023).
The fully-coupled WRF-Hydro model is an atmospheric-hydrological modeling system, which is comprised of the atmospheric WRF model (Weather Research and Forecasting model;Skamarock et al. (2019aSkamarock et al. ( , 2019b))) and its hydrological extension, the WRF-Hydro model (Gochis et al., 2018).It has been applied in different regions at a range of spatial and temporal scales.The coupled model performance was examined, for instance, for mediumsized catchments (800-2210 km 2 ; Givati et al., 2016;Senatore et al., 2015;Wang et al., 2020) and for medium to large catchments (3,280 km 2 ; Kerandi et al., 2017).The interest in the application of the fully-coupled WRF-Hydro model has grown, mainly for its potential to improve the simulation of terrestrial hydrometeorological variables compared to the standalone WRF atmospheric model.A number of studies showed that the coupled WRF-Hydro increases the domain-average soil moisture, evapotranspiration and latent heat fluxes and reduces the surface runoff and sensible heat fluxes in multi-month or multi-year simulation periods (Rummler et al., 2019;Senatore et al., 2015).These differences can be explained by the description of the lateral surface and subsurface flow of water by WRF-Hydro, which can then infiltrate and/or return to the atmosphere as evapotranspiration before reaching the stream.Zhang et al. (2021a) found that the impact of coupled modeling on surface energy fluxes was pronounced in areas of complex terrain in an alpine catchment of the Tibetan plateau.These authors found that the coupled modeling had a minor impact on net radiation and ground heat flux.
The representation of vegetation density and the evapotranspiration equations also regulate the energy and water balance components (Cuntz et al., 2016;Zheng et al., 2019) among other parameterizations of the Noah (Ek et al., 2003) and Noah-MP LSMs (Niu et al., 2011;Yang et al., 2011) in WRF Hydro.For instance, Zheng et al. (2019) showed that the Jarvis parameterization for transpiration (Jarvis, 1976) simulates the evapotranspiration fluxes better than the default Ball-Berry parameterization (Ball et al., 1987) in Noah-MP for semi-arid climate zones of the US.Sofokleous et al. (2023) also found that a modified version of the Jarvis parameterization fitted transpiration observations better and contributed to improved streamflow simulations for the Mediterranean island of Cyprus.However, the combined impact of these vegetation-related parameterizations with that of lateral routing of water in fully-coupled WRF-Hydro simulations has not been studied.
The impact of lateral terrestrial routing on precipitation has gained attention in previous fully-coupled WRF-Hydro studies, because of the possible feedbacks enhanced by the coupling of the atmosphere with the hydrological processes.A reduction of the bias of the total simulated precipitation (2.7%-10%) by coupled WRF-Hydro, compared to WRF, was found for a number of storm events in central Israel (Givati et al., 2016) and in Attica, Greece (Galanaki et al., 2021).Zhang et al. (2019) reported a modification in the spatial patterns of precipitation due to the lateral routing, particularly over mountain areas.Arnault et al. (2016) concluded that the impact of lateral flow and runoff -infiltration partitioning on precipitation in WRF-Hydro, was more evident in a domain with a size at the scale of the catchment area, that is, ∼100 × 100 km 2 , than in a larger domain, from their simulations for the 12,800 km 2 Sissili catchment in West Africa.Arsenault et al. (2018) found that the normalized ensemble spread of daily precipitation was higher in coupled simulations in areas where there is high topographic variability and where the weather conditions are dominated by local processes.Furthermore, a number of combined terrestrial-atmospheric water balance studies revealed that simulated precipitation is mostly dominated by moisture advection from outside the model domains, with a limited contribution of the local soil moisture (e.g., Kerandi et al. (2017); Zhang et al., 2019).
Different studies addressed the uncertainty in the representation of different components of the coupled WRF-Hydro by the ensemble modeling approach.Arsenault et al. (2018) and Zhang et al. (2021b) depicted the uncertainty of the planetary boundary layer (PBL) through alternative PBL parameterizations.Arsenault et al. (2018) found that there is a higher ensemble skill improvement in modeled precipitation, up to 20% with various PBLs, than with variations in lateral flow achieved with different values of the infiltration parameter (REFKDT) of the hydrological model.Rummler et al. (2019) constructed two ensembles, comprised of daily lagged model initializations of the standalone WRF and of the coupled WRF-Hydro.They found enhanced variability of the simulated precipitation and trajectory of a low-pressure system in the coupled model runs.With the standalone version of the WRF model, previous studies found that precipitation is affected by commonly used WRF atmospheric parameterizations, namely, cumulus cloud, PBL and microphysics parameterizations with multi-physics simulation experiments (Avolio & Federico, 2018;Sofokleous et al., 2021;Zittis et al., 2017).The impact of using combinations of the atmospheric parameterizations in a multi-physics atmospheric ensemble to be used as forcing for WRF-Hydro streamflow simulation was not examined in previous WRF-Hydro studies.
This study aims to quantify the impact of atmospheric parameterizations, land surface parameterizations and terrestrial lateral flow processes on hydrometeorological variables simulated by the fully-coupled WRF-Hydro model.The effect of these three parameterizations on the components of the energy and water balance of the fully-coupled WRF-Hydro model has not been investigated together in the literature.Furthermore, previous fullycoupled WRF-Hydro studies have focused on medium and large-sized watersheds, with areas above 1,000 km 2 .In this study, the standalone WRF, fully-coupled WRF-Hydro and one-way coupled WRF-Hydro are employed with different physical parameterizations in simulation experiments for the Eastern Mediterranean island of Cyprus and 31 small watersheds (order of 10-100 km 2 ).The first objective of the study is to examine the role of the terrestrial process parameterizations by comparing surface energy and water balance components simulated by the standalone WRF with those simulated by two configurations of the fully-coupled WRF-Hydro model (one with default WRF-Hydro parameterizations and another with modified vegetation-related parameterizations).The second objective of the study is to examine the role of atmospheric process parameterizations, by evaluating the performance of different one-way coupled WRF-Hydro streamflow simulations driven by five different multiphysics WRF simulations.The simulated variables are evaluated for the entire island, the area of the Troodos mountains and the areas of 31 watersheds for the period 2011-2013.

Study Area
The simulation domain focuses on Cyprus, an island in the Eastern Mediterranean Sea, at latitude 35°N and longitude 33°E (Figure 1).The Troodos massif covers more than 25% of the island and stretches from approximately 300-1,951 m.Its east-west-elongated, dome-like structure is a lithologically complete ophiolite and forms the predominant topographic feature of the island (Christofi et al., 2020).The soils of Troodos are characterized by stony gravelly texture (Camera et al., 2017) and the land use is dominated by natural vegetation, that is, sclerophyllous vegetation and coniferous forests (Zoumides et al., 2017).The plains and coastal areas are characterized by sedimentary soil types and agricultural and urban land use.The 31 watersheds of the study form a radial drainage system around Troodos (Figure 1).The areas of the 31 watersheds range from 5 to 115 km 2 and their average slope ranges from 11% to 53%.
The climate of the island is characterized by the typical for Mediterranean climates annual and interannual variability in precipitation (Hoerling et al., 2012).The mean annual precipitation  ranges from 240 to 400 mm in the plains and coastal areas to nearly 1,100 mm on the highest elevations of Troodos Mountains (Camera et al., 2014).December and January are the wettest months and about 80% of total annual precipitation occurs between November and April.The near surface temperatures, spatially averaged, range between 5°C and 14°C in January and February, and between 21°C and 33°C in July and August.Snowfall occurs in the coldest winter precipitation episodes, which have high interannual variability.Snowcover is spatially limited on the highest elevation area and it affects a small area of 4 of the 31 of the modeled watersheds.Reference evapotranspiration computed with the FAO Penman-Monteith equation (Allen et al., 1998) exceeds 1,000 mm per year for the largest part of the area.The runoff coefficients of the 31 watersheds, namely the ratio of total streamflow to total precipitation, range from 4% to 57%.

Model Description
The state-of-the-art WRF model (version 4.4) is used in the standalone WRF and fully-coupled WRF-Hydro simulation experiments of this study.WRF is a mesoscale numerical weather prediction system widely used by the atmospheric research community for both research and operational forecasting at a broad range of spatial scales (10 m-10 6 m) (e.g., Constantinidou et al., 2019;Liu et al., 2017).Different physics schemes in the model code are used to parameterize physical processes, namely the processes on the land-surface, in the planetary boundary layer, the atmospheric and surface radiation, microphysics-cloud processes and cumulus convection (Skamarock et al. (2019a(Skamarock et al. ( , 2019b))).
WRF-Hydro is a distributed physically-based hydrologic model that simulates terrestrial hydrologic processes, that is, surface, subsurface and channel flow.WRF-Hydro enhances the traditional land surface modeling of processes simulated in the vertical direction, by the component of lateral routing of surface and soil water.Lateral subsurface runoff in WRF-Hydro is calculated with a method that approximates the hydraulic gradient between adjacent grid cells and that includes the effect of topography, saturated hydraulic conductivity and saturated soil layers.Laterally routed surface runoff is subject to reinfiltration.The LSMs used with WRF-Hydro in this study are the Noah LSM (Ek et al., 2003) and Noah-MP LSM (Noah LSM with multiple parameterization options; Niu et al., 2011;Yang et al., 2011).The Noah LSM is used in the standalone WRF and fully-coupled WRF-Hydro simulations and the Noah-MP LSM is used in the one-way coupled WRF/WRF-Hydro simulations.The two LSMs are used with the same parameters and parameterization options.The other WRF-Hydro modules are (a) an overland flow module; (b) a subsurface flow module; (c) a channel flow module; (d) a reservoir storage and flow module; and (e) a groundwater (bucket) module.WRF-Hydro can operate in an offline mode with meteorological forcing provided by gridded observations, by reanalysis data or by the output of atmospheric models, such as the WRF model.In the coupled mode of WRF and WRF-Hydro, the two models interact through the exchange of fluxes at the land surface during the simulation.
Since WRF-Hydro version 5.0, a new capability was added to allow specification of all land surface and lateral routing model parameters on the model grid (Gochis et al., 2018).However, this capability is not supported in the fully-coupled WRF/WRF-Hydro even with the latest WRF-Hydro version 5.2.In the coupled mode, a number of parameters is specified as a single global value (e.g., REFKDT and REFDK), the soil drainage parameter (SLOPE) is specified by land slope classes, whereas other parameters are specified by soil or land cover classes.
The WRF-Hydro model version 5.2 is used in the coupled simulation experiments.For the offline WRF-Hydro simulation experiments, the WRF-Hydro version 5.1 is used.

Model Set-Up
A three-nested domain at 12 km, 4 and 1 km with time steps equal to 72 s, 24 and 6 s, respectively, was defined for the fully-coupled and standalone WRF simulations.The number of model vertical levels was set to 40.Initial and boundary conditions were obtained from the ERA5 reanalysis dataset at 31-km resolution.The WRF configuration used here is described in detail by Sofokleous et al. (2021).WRF-Hydro was activated for the innermost domain of WRF, which has an area of 225 × 145 km 2 (Figure 1).The spatial resolution for the LSM was set to 1km and for the lateral routing processes, it was set to 100 m.
The simulation experiments were initialized at the beginning of each month.Precipitation simulations with a monthly model initialization frequency were found to result in equal model performance as those with more frequent (5-day) initializations (Sofokleous et al., 2021).Additional simulation experiments with less than monthly initialization frequency were found to reduce model performance.A spin-up time of 24 hr was found to be adequate for the atmospheric model (Sofokleous et al., 2021).A spin-up of the land surface model was avoided by manually initializing the states of the model at the beginning of the simulation period.This was possible because the start of the simulation period coincides with the end of the summer period in the study area, when soil moisture is close to the wilting point and streamflow is zero in most watersheds.Thereafter, for each new month of simulations, soil moisture conditions and streamflow were initialized from the final states of the previous month.

Set Up for the Comparison of Standalone WRF With the Fully-Coupled WRF-Hydro
Land surface energy and water balance components simulated by the standalone WRF (WRF.sa) were compared against the simulated components by two configurations of the fully-coupled WRF-Hydro model (WRFH.cpl1and WRFH.cpl2).The WRF atmospheric parameterizations in WRFH.cpl1 and WRFH.cpl2 configurations are the same as WRF.sa.A summary of the model configurations and selected physics parameterizations for the simulation experiments are presented in Table 1.The atmospheric physics parameterizations used in this study were selected based on the results of Sofokleous et al. (2021) for multi-month simulations and Zittis et al. (2017) for extreme events for the study area of Cyprus.The selection of the LSM and WRF-Hydro model parameters in WRFH.cpl1 and WRFH.cpl2 is based on the model calibration study of Sofokleous et al. (2023).
Initial trial-and-error runs of the fully-coupled model (not shown), with parameter values of the land surface model same as the calibrated values in the offline WRF-Hydro presented in Sofokleous et al. (2023), for the hydrological year 2011-2012, produced unexpectedly low streamflow, about 95% less flow than the observed streamflow.A major difference between the two configurations is the time step of the models, namely 6 s in the coupled simulation and 30 min in the offline simulation.Senatore et al. (2015) discussed the differences in the model response due to the different time steps in the coupled and offline models.To reduce this streamflow underestimation, we tested a refkdt value of 0.1 in both WRFH.cpl1 and WRFH.cpl2, instead of the default value of 3.0, used with the WRF.sa.With this parameter change, the median percent bias decreased to 17% for WRFH.cpl1 and to 21% for WRFH.cpl2,relative to the observed streamflow.
The second coupled WRF-Hydro configuration, WRFH.cpl2, was further adjusted with changes in the vegetation parameterizations of Noah-LSM.The default leaf area index (LAI) values were replaced with the Copernicus Global Land Service LAI dataset, which is available at 100-m resolution over the study area (Table S1 in Supporting Information S1).Similar to Sofokleous et al. (2023), the 100-m resolution Copernicus product was averaged for the four dominant land use categories of the study area for each month.The Copernicus LAI values were found to be 50%-80% lower than the default model LAI values.A second adjustment in WRFH.cpl2 is the modification of the radiation component of the Jarvis transpiration parameterization in Noah LSM (Jarvis, 1976),  Kain, 2004), Betts-Miller-Janjć (2; Janjć 1994), Grell-Freitas (3;Grell & Freitas, 2014).c Planetary boundary layer physics schemes: Yonsei University (1; Hong et al., 2006), Mellor-Yamada-Janjć (2; Janjć 1994).d Surface layer physics schemes: Eta similarity (2; Janjć 1994), MM5 similarity scheme (91; Zhang & Anthes, 1982).e refkdt is the infiltration parameter, ksat is the soil hydraulic conductivity, LAI is the leaf area index and Jarvis.noc is the Jarvis transpiration parameterization with a modified radiation component to represent nocturnal transpiration.
to account for nocturnal transpiration.This adjustment caused an average increase in transpiration of 15% relative to the default Jarvis parameterizations (Sofokleous et al., 2023).
The two fully-coupled WRF-Hydro simulations (WRFH.cpl1,WRFH.cpl2) were compared with the standalone WRF simulation (WRF.sa) for the time series of the daily water balance components of soil moisture, evapotranspiration and rainfall, as well as 2-m daily mean temperature and for the diurnal cycles of hourly ground, latent and sensible heat fluxes for the period from September 2011-August 2012.These components were averaged for two areas; the entire island and the land cells of the rectangular domain around Troodos (Figure 1).

Set Up for the Analysis of the Impact of WRF Atmospheric Parameterizations on WRF-Hydro Streamflow
The impact of atmospheric forcing on simulated streamflow was examined with five WRF configurations (WRF1-5).These five WRF configurations, referred to as members, were used as forcing in five one-way coupled WRF and WRF-Hydro simulations, denoted as WRFH1-5.The forcing of WRFH1-5 was at an hourly time step.The physics parameterizations in WRF1-5 (Table 1) were identified as the most skilful over the study area from a set of 18 configurations (Sofokleous et al., 2021).The configuration of WRF-Hydro calibrated for streamflow and evapotranspiration observations, as described in Sofokleous et al. (2023), is used for these simulations.
WRF1-5 were evaluated for their skill in reproducing the observed gridded precipitation dataset CY-OBS over the island and over the areas of 31 watersheds and WRFH1-5 were evaluated for their skill in reproducing the observed streamflow for 31 watersheds (Figure 1) for a 2-year period (September 2011-August 2013).The simulated streamflow from the offline WRF-Hydro simulation, forced by observed precipitation (WRFH.obs),was considered as the baseline for the comparison of the impact of the five atmospheric members on streamflow.
The evaluation measures used in these analyses are presented in Section 4.3.

Evaluation Measures
The evaluation measures for precipitation are four in total.Two measures are categorical scores, that is, the Extreme Event Score (EES; Sofokleous et al., 2021) and the Peirce Skill Score (PSS; Stephenson, 2000), with the 30-mm daily rain as the threshold.The other two measures are the Relative bias of total precipitation (mm/mm) and the Mean Absolute Error of daily precipitation averaged for all grid cells within the boundaries of each watershed (MAE; mm/day).For simulated streamflow, the evaluation measures used are the Nash-Sutcliffe Efficiency (NSE; Nash & Sutcliffe, 1970) and the Kling-Gupta Efficiency (KGE; Kling et al., 2012), which evaluate the goodness-of-fit of simulated to observed streamflow.The other two measures are the MAE of daily flows (mm/day) and the Relative bias of total streamflow in the simulation period (mm/mm).A relative bias larger than unity indicates overestimation by the model.
The relative skill of each member of WRFH1-5 in simulating precipitation and in simulating streamflow, in each watershed, was measured by a Composite Scaled Score (CSS; Sofokleous et al., 2021): where i is the index identifying the member, s is the index of the statistical measure out of a number of N s (four) measures, Y s ,i is the value of measure s obtained by member i and Y s,worst and Y s,best are the worst and the best values for measure s among the five members.The score ranges between 0 and 1.

Comparison of Standalone WRF With the Fully-Coupled WRF-Hydro
The daily volumetric soil water content (SM) and evapotranspiration (ET) in the simulation period in Figures 2 and 3 exhibit a generally similar seasonal pattern, with higher day-to-day variation for ET for all the three model configurations (WRF.sa,WRFH.cpl1,WRFH.cpl2).The peak time of ET, around April, is shifted, however,  about 1 month later than that of the SM.This can be explained by the fact that the vegetation transpiration, which is a component of ET, depends also on other environmental conditions, such as temperature.
The average daily SM for WRF.sa is 0.17 (0.20) m 3 m 3 for the entire island (Troodos mountains), while WRFH.cpl1 and WRFH.cpl2 yield 17% (9%) and 22% (16%) more SM, respectively.The monthly average SM for the island is higher in the coupled model runs than in the standalone model in at least 10 months out of the 12 months of simulations.The analysis for the monthly results for the island and for the Troodos Mountains is shown in Table S2 in Supporting Information S1.The average daily ET is 0.84 mm (0.65 mm) for WRF.sa, 1.12 mm (0.81 mm) for WRFH.cpl1 and 0.99 mm (0.77 mm) for WRFH.cpl2.For the annual total amounts, WRFH.cpl1 and WRFH.cpl2 have 33% (24%) and 18% (18%) higher ET than the standalone WRF model.The monthly ET totals for the island of WRFH.cpl1 and WRFH.cpl2 are higher than those of WRF.sa in 12 and 11 months, out of the 12 months, respectively.
The increased ET in the coupled simulations is a direct result of the higher soil moisture.This can be attributed to the lateral routing in WRFH.cpl1 and WRFH.cpl2, which allows surface and soil water to move across grid-cells, as opposed to the flux of water in the land surface model of WRF.sa, where surface runoff is directly removed from the terrestrial system.This difference can indicatively be seen in the month of November (Figure 2), when the transition from the dry to wet period happens (Figure 4).The domain-average soil moisture in WRF.sa does not respond as directly to rain as in coupled WRF-Hydro.
The comparison of precipitation for the three model runs in Figure 4 shows no patterns of a systematic variation among the three.The comparison on the monthly amounts (Table S2 in Supporting Information S1) indicates that precipitation with the coupled runs WRFH.cpl1 and WRFH.cpl2 is more than that with WRF.sa in 11 and 10 months out of the 12 months of the simulation period, respectively.The increase in total precipitation is 3% (0.7%) for WRFH.cpl1 and 1.5% (0.9%) for WRFH.cpl2relative to the amount of 658 mm (923 mm) simulated by WRF.sa over the island (Troodos area).This minor increase in precipitation with the coupled runs could be the signature of increased vapor flux due to increased ET to the atmosphere from the lateral routing.A small positive change (1.4%-2.7%)for coupled runs was also found by Givati et al. (2016) in two wet months of 2013 for simulations in central Israel.Other studies found that the coupled simulations do not result in significant differences in simulated annual precipitation, compared to standalone WRF (Zhang et al., 2019).A larger effect of the modified surface fluxes to simulated precipitation could be expected in cases where the domain sizes are large enough to resolve mesoscale weather phenomena, as Arnault et al. (2016) suggested.
Compared to observed precipitation, fully-coupled WRF-Hydro and standalone WRF perform similarly, with a total precipitation bias from 10% to 13% for the island and 22%-23% for the Troodos area.These results indicate that lateral routing is not expected to yield any considerable change in simulated precipitation in a relatively small simulation domain and possibly a study area surrounded by sea, such as the one for Cyprus (∼3,300 km 2 ), where most of precipitation originates from large-scale weather events.Similar to simulated precipitation, the simulated temperature (Figure S1 in Supporting Information S1) is not systematically different between the standalone and coupled simulations.The difference in simulated 2-m daily mean temperature is less than 0.5°C between the three configurations.The variation in simulated temperature among the three simulations is smaller than their bias relative to the observed temperature time series, that is, about 1°C and 3°C, on average, over the island and over the mountains.This temperature bias could be attributed to the selected PBL parameterization (García-Díez et al., 2012).
The pattern of daily cycles of latent heat flux (LH) for the three simulations in Figure 5 shows that the simulation with lateral routing, that is, WRFH.cpl1 has the largest value of average daily LH throughout the year.It has 36% higher LH on average than that of WRF.sa (32 W/m 2 ).This difference is maximized in April, when WRFH.cpl1 has almost double the LH of WRF.sa.From these results, it becomes apparent that the lateral routing of water in the coupled model has implications, through the increased soil moisture, on both the water and the energy balance, as seen with the increases in both ET (Figure 3) and LH (Figure 5), respectively.These results are in agreement with previous studies for another Mediterranean area, in southern Italy (Senatore et al., 2015), as well as for the eastern Alps (Rummler et al., 2019).The diurnal cycle of sensible heat flux (SH) (Figure S2 in Supporting Information S1l) exhibits smaller differences among the three simulations, compared to the differences seen with LH.For the entire year, the average daily SH is 66 W/m 2 for WRF.sa.The average daily SH is 7% lower for WRFH.cpl1 and 13% higher for WRFH.cpl2.The lower SH in the simulation WRF.cpl1, which only differs from WRF.sa by the lateral routing, is in agreement with the findings of Senatore et al. (2015) and Zhang et al. (2021a).
The coupled simulation with the modified LAI and modified transpiration parameterization, that is, WRFH.cpl2 has an annual average daily LH about 10% higher than WRF.sa.For WRF.cpl2, the potentially higher transpiration resulting from the modified transpiration parameterization was counter-balanced by the lower, by 50%-80%, LAI values of the Copernicus database, which resulted in a lower ET and LH than WRF.cpl1.The comparatively lower ET and LH imply lower soil water use for soil evaporation and vegetation transpiration, which can also be seen in the up to 5% higher soil moisture of WRFH.cpl2 relative to WRFH.cpl1.This result indicates that the lateral routing process can lead to changes in modeled fluxes that are comparable to the change of the fluxes due to modifications in vegetation parameterizations.This finding confirms and extends the findings of Cuntz et al. (2016) and Sofokleous et al. (2023) who found that runoff and ET depend on both soil and vegetation parameters.Similar to here, change in simulated soil moisture and other basin-scale hydrological losses, due to alternative representations of vegetation descriptors, was found using different NVDI-LAI regression models, by Gigante et al. (2009), and using a biophysically-derived vegetation fraction from LAI, by Fang et al. (2018).
The ground heat flux (GH) (Figure S3 in Supporting Information S1) exhibits also variation among the three runs, although smaller than that for SH and LH.GH has also the smallest magnitude in its diurnal cycle.The average daily GH for all months is 0.9 W/m 2 with WRF.sa, 0.9 W/m 2 with WRFH.cpl1 and 0.4 W/m 2 WRFH.cpl2.The daily average value of the net surface energy balance of SH, LH and GH ranges from 28.8 W/m 2 in December to 163 W/m 2 in June and it has an average value for the year of 98 W/m 2 with WRF.sa.The net energy of the two coupled simulations was about 10% higher than that of the standalone WRF.This increase in net energy could be indicative of a modified net radiation budget in the coupled model.However, the difference in the downward and outgoing shortwave and longwave radiation at and from the land surface in the coupled model runs relative to WRF.sa were found to be minor, that is, between 0.7% and 0.4% (Table S3 and Figures S4-S6 in Supporting Information S1).A slight decrease in incoming shortwave radiation ( 0.2% with WRFH.cpl1 and 0.7% with WRFH.cpl2) could be associated with the also small increase in precipitation in the coupled simulations, up to 3%, through increased cloudiness.Overall, the driver for net energy budget change in the coupled simulations seems to be the soil moisture content, which in turn induces small variations in the radiation budget components.

Impact of WRF Atmospheric Parameterizations on WRF-Hydro Streamflow
The selected WRF physical parameterizations (WRF1-5) capture, on average, the observed precipitation for the area of the island for a 2-year period, as seen in Figure 6.Members WRF1, WRF4 and WRF5 overestimate total precipitation by 4%, 16% and 13%, whereas WRF2 and WRF3 underestimate observed precipitation by 18% and 12%, respectively.These findings agree with the direction of bias of the five members reported in Sofokleous et al. (2021), who used a shorter, 8-month, evaluation period.There are also spatial differences among the five members.For instance, WRF1, with the lowest overall bias, overestimates the total-period precipitation on different parts of the highest elevation areas of the mountains, by 25% up to 50%, and it has a negative, but close to zero bias on the remaining part of the mountains.These spatial differences that exist to a different degree for different members can be indicative of the variable performance of WRF-Hydro, forced by WRF, for different watersheds.
The evaluation of WRF modeled precipitation over the areas of the 31 watersheds by four evaluation metrics is shown in Figure 7.The skill in the reproduction of extreme events, above the 30 mm/day threshold, is generally consistent, with a median EES around 0.3.With the PSS and for the same precipitation threshold, WRF4 and WRF5 achieve a slightly better performance (median PSS of 0.6) than the other three members (median PSS of 0.4).Yet, the result of higher PSS of WRF4 and WRF5 might be inflated by the overestimation of precipitation, including the extreme events, as seen by the generally larger MAE and relative bias values for these two members.The MAE is around 1.7 mm/day and the median relative bias ranges from about 0.2 to about 0.2 for the five members.Overall, there is no member that achieves statistically better performance across all watersheds with either the EES, the PSS or the MAE.The most pronounced variation among the members is for the total period precipitation, as shown by the relative bias.
The ranking of the five members, based on the composite score of the four metrics, that is, the CSS, from highest to lowest performance is WRF2, WRF3, WRF5, WRF4 and WRF1.The simulated precipitation at the convection-resolving resolutions of this study is mainly affected by the microphysics scheme, as shown in the comparison of the CSS of WRF1 (Ferrier) with WRF2-5 (WRF double-moment 6th class) in Figure 7.The best performing WRF members (WRF2, WRF3 and WRF5) have in common the microphysics scheme parameterization (WDM6), in agreement with the findings in Zittis et al. (2017) and Sofokleous et al. (2021) for the area of Cyprus.In contrast to the findings for the area of Cyprus in the Eastern Mediterranean, Cassola et al. (2015) found that the WRF single-moment 6th class (WSM6) outperformed WDM6 for the simulated rainfall in extreme convective events in Liguria region, Italy, in the central Mediterranean.Jeworrek et al. (2019) also found that the less computationally expensive single-moment 5th class (WSM5) yielded better verification scores than other double-moment schemes for forecasts of orographic precipitation in British Columbia, Canada.Jankov et al. (2007) found large differences in precipitation amounts for simulations of orographic precipitation in California with the Ferrier microphysics schemes compared to other schemes, similar to the differences in the relative biases of the Ferrier (WRF1) and WDM6 scheme (WRF2-5) in Figure 6.The best performing WRF2 and WRF3 have also in common the PBL scheme (Yonsei University).The difference in the best performing WRF2 and WRF3 is in the cumulus scheme used, the Kain-Fritch and the Betts-Miller-Janjić, respectively.The good performance of Betts-Miller-Janjić was also reported in Avolio and Federico (2018) and Sofokleous et al. (2021).Figure 8 shows five performance metrics for the streamflow of five WRF-Hydro runs (WRFH1-5) forced by the output of the five different WRF configurations (WRF1-5), based on different atmospheric physics parameterization schemes, shown in Table 1, as well as the WRF-Hydro run forced by the observed precipitation dataset CY-OBS (WRFH.obs).Despite the fact that the same WRF-Hydro configuration is used in the five-member oneway coupled runs (WRFH1-5), there is a large discrepancy in the performance of the WRF-Hydro members (Figure 8).In fact, the ranking of the WRF and WRF-Hydro members, based on the CSS, shown in Figures 7 and  8, respectively, follows the same order.In addition, the sign of the relative bias of precipitation for the areas of the 31 watersheds is consistent with the sign of the bias of streamflow for the same members.WRF2 (WRFH2) and WRF3 (WRFH3) simulate precipitation (streamflow) with negative bias, and WRF4 (WRFH4) and WRF5 (WRFH5) return positive bias for the two variables.
WRFH2 and WRFH3, forced by WRF2 and WRF3, respectively, achieve the highest NSE and KGE for the 31 watersheds, with a median NSE value of 0.2 and 0.02, respectively, and the lowest errors, with a median MAE of about 0.3 mm/day and a median relative bias from 0.3 to 0.1.The higher performance of the two WRF-Hydro members is also seen in the median CSS, which is above 0.8 for both.Streamflow is simulated by WRFH1, WRFH4 and WRFH5 with negative efficiency metrics for most watersheds, a median MAE larger than 0.4 mm/day and median relative bias of total streamflow from 0.5 to 1.The CSS of these three members is below 0.6.
The comparison of WRFH1-5 with the WRF-Hydro simulation forced by observed precipitation (WRFH.obs) in Figure 8 shows that all WRF-Hydro members have degraded performance relative to WRFH.obs.The range of the values of the four metrics for the 31 watersheds indicates that, even with WRFH.obs, streamflow is not always simulated with high performance.These errors are attributed to the hydrological model component, such as low resolution of terrestrial datasets, imperfect model parameterizations and model calibration aspects, as discussed by Sofokleous et al. (2023).The forcing of WRFH.obs, that is, the gridded CY-OBS precipitation dataset, can also be a source of errors, because this dataset was generated from the spatial interpolation of station observations in an area characterized by the high topographic variability of the Troodos mountains.The evaluation metrics for the simulation of streamflow obtained for WRFH3 that achieved the highest CSS, and for WRFH.obs, for each watershed for the two-year simulation are shown in Figure S7 in Supporting Information S1.
The streamflow simulated with the one-way coupled WRF-Hydro for three watersheds with different prevailing hydrological conditions, that is, non-ephemeral (st13) and ephemeral (st17 and st21), shown in Figure 9, follows the patterns of the observed streamflow.The hydrographs of all 31 watersheds are presented in Figure S8 in Supporting Information S1.In general, there is a variable WRF-Hydro model performance with positive and negative NSE and KGE for different members.For st13, with a baseflow index of 0.96, streamflow is almost similarly represented by all five members, except for peak flows.An overestimation of peak flows is seen in all watersheds, especially with the forcings of members WRF1 and WRF4, with generally negative and overall lower values of NSE.The watershed with lowest observed streamflow, st21, is simulated with the lowest overall performance by the five members.The low performance can be attributed, in addition to errors in the simulated precipitation, to the hydrological model parameterization, especially the groundwater parameterization in WRF- Hydro.The average baseflow index is 0.72 for the 31 watersheds, which indicates the significance of baseflow simulation, in addition to precipitation.
These results from Figures 7, Figures 8 and 9 indicate that both the hydrological and the atmospheric model contribute to the streamflow biases.Regarding the atmospheric component, other studies for Mediterranean watersheds, for instance, also found degraded performance in modeled streamflow with the coupled model relative to simulations with observed meteorological forcing.Senatore et al. (2015) reported, for the Crati river basin (1,281 km 2 ) in southern Italy, an NSE above 0.7 for the calibration and validation period with observed meteorological forcing and an NSE of 0.25 for a three-year simulation with the coupled WRF-Hydro.Givati et al. (2016) reported an NSE of 0.8 for two storm events in Ayalon Basin (800 km 2 ) in central Israel with observed precipitation forcing and an NSE equal to 0.6 with the coupled WRF-Hydro.
The lower efficiency of the coupled WRF-Hydro in this study compared to other studies on Mediterranean watersheds, for both the simulations forced by observed precipitation (NSE = 0.49) and the simulations forced by WRF output (NSE from 0.2 to 0.2) can be attributed to the small size of the watersheds (5-115 km 2 ).Camera et al. (2020) reached similar conclusions in the same study area, for the simulation of two extreme events with WRF-Hydro forced by observed precipitation, but found more watersheds with negative NSEs, compared to this study.

Summary and Conclusions
This study investigated the simulation of terrestrial energy and water balance components of the coupled atmospheric -hydrologic WRF-Hydro model in a geographically small domain (225 × 145 km 2 ) and 31 small watersheds (5-115 km 2 ) with different terrestrial and atmospheric parameterizations.Two sets of simulations were conducted for the Mediterranean island of Cyprus for two hydrological years (2011-2013).
In the first set of simulations, the comparison of the standalone WRF model with the fully-coupled WRF-Hydro showed modifications in components of the energy and water balance due to the process of lateral routing of water in the coupled mode.The most important changes of the fully-coupled WRF-Hydro, relative to standalone WRF, resulting from the lateral routing, are a 17% increase in soil moisture, a resulting 36% increase in latent heat flux and 33% increase in ET, and a 3% increase in precipitation.A modification of the transpiration parameterization to include nocturnal transpiration and the use of lower LAI values, as obtained from the Copernicus database, in the coupled WRF-Hydro model, showed a 22% increase in soil moisture, a resulting 10% increase in latent heat flux and 18% increase in ET, and a 1.5% increase in precipitation, relative to standalone WRF.
In the second set of simulations, the effect of WRF simulated precipitation on WRF-Hydro streamflow was studied with five one-way coupled multi-physics WRF and WRF-Hydro simulations (members).The selection of microphysics schemes was found to have a major impact on the simulated precipitation in convection resolving spatial resolutions (1-km used in this study).The WRF Double moment 6th class was found to be the best performing scheme.The relative ranking of the five WRF members for the reproduction of precipitation and the five WRF-Hydro members for the reproduction of streamflow was found to be the same.The bias in total precipitation of the five WRF members varied between 18% and 16% and the bias in the streamflow of the five WRF-Hydro members varied between 37% and 119%, compared to the bias of 6% for WRF-Hydro forced by observed precipitation.
The routing of water in the lateral direction is especially important for areas with high topographic variability.However, in the standalone Noah-MP LSM and the WRF atmospheric model, the vertical flow of water is the only direction of flow simulated, whereas surface runoff is removed from the terrestrial water balance.This study showed that the process of lateral routing of water has a large impact on components of the terrestrial energy and water cycle and extends similar findings from studies in different climate areas and larger domains to a 9,251-km 2 Mediterranean island and small mountain watersheds.Therefore, the process of lateral terrestrial routing should be considered in future applications with distributed, numerical climate, weather or hydrological models for any domain size, to improve the simulation of the terrestrial water cycle, land surface fluxes and atmospheric processes.Although not examined here, the improved simulation of the water cycle in the coupled WRF-Hydro could also enhance the modeling of flood events.
Despite the high-resolution atmospheric model simulations at 1-km, the modeled precipitation contained biases that significantly affected the performance of the calibrated hydrological model in the small watersheds of this study.The simulation of streamflow in small Mediterranean mountain watersheds requires the use of an optimized atmospheric ensemble.Such an ensemble could comprise of a subset of the overall best performing multiphysics atmospheric model configurations, possibly, combined with different initialization times.
For future coupled atmospheric-hydrological studies in small watersheds, where there is a large uncertainty in the precipitation forcing, a skilful hydrological model should be used in order to minimize the uncertainty originating from the hydrological model parameterizations.The use of the WRF-Hydro model with distributed parameters, which is currently not available for use in the fully-coupled WRF/WRF-Hydro, could also improve the simulation of hydrologic processes.

Figure 1 .
Figure 1.Elevation map of Cyprus, with the extent of the map corresponding to the innermost domain of the three-nested WRF model domain, and with the rectangular interior box representing the extent of the Troodos mountains (left).Boundaries and drainage network of 31 watersheds of the Troodos mountains for the WRF-Hydro streamflow simulations (right).

Figure 3 .
Figure 3.As in Figure 2 for the simulated daily evapotranspiration (mm/d).

Figure 6 .
Figure 6.Total observed gridded precipitation (CY-OBS) and percent bias of total simulated precipitation for the five WRF members for the period Sep 2011-Aug 2013.

Figure 7 .
Figure 7. Boxplots of Extreme-Event Score (EES), Peirce Skill Score (PSS), mean absolute error (MAE), relative bias and Composite Scaled Score (CSS) for the five WRF members, with each box representing the variability of the metric for the spatially-average daily precipitation of the 31 watersheds in the period Sep 2011-Aug 2013.

Figure 8 .
Figure 8. Boxplots of the 31 watersheds for the Nash-Sutcliff efficiency (NSE), Kling-Gupta efficiency (KGE), mean absolute error (MAE), relative bias and Composite Scaled Score (CSS) of WRF-Hydro simulated streamflow with precipitation forcing from observations (WRFH.obs)and from the five WRF members (WRFH1-WRFH5).Each box represents the variability of the values of each metric for the 31 watersheds in the period Sep 2011-Aug 2013.The dashed green line represents the best possible value for each metric.

Figure 9 .
Figure 9. Observed streamflow (OBS; red dashed line), WRF-Hydro modeled streamflow with forcing from precipitation observations (WRFH.obs;black line) and range of hydrographs of the five one-way coupled WRF-Hydro members (WRFH1-5; gray shading) in m 3 s 1 of three watersheds.The Nash Sutcliff efficiency (NSE) and Kling-Gupta efficiency (KGE) evaluation metrics for the daily streamflow for the period September 2011-August 2013 are given for the six simulations.