Precipitation is affected by soil moisture spatial variability. However, this variability is not well represented in atmospheric models that do not consider soil moisture transport as a three-dimensional process. This study investigates the sensitivity of precipitation to the uncertainty in the representation of terrestrial water flow. The tools used for this investigation are the Weather Research and Forecasting (WRF) Model and its hydrologically enhanced version, WRF-Hydro, applied over central Europe during April–October 2008. The model grid is convection permitting, with a horizontal spacing of 2.8 km. The WRF-Hydro subgrid employs a 280-m resolution to resolve lateral terrestrial water flow. A WRF/WRF-Hydro ensemble is constructed by modifying the parameter controlling the partitioning between surface runoff and infiltration and by varying the planetary boundary layer (PBL) scheme. This ensemble represents terrestrial water flow uncertainty originating from the consideration of resolved lateral flow, terrestrial water flow uncertainty in the vertical direction, and turbulence parameterization uncertainty. The uncertainty of terrestrial water flow noticeably increases the normalized ensemble spread of daily precipitation where topography is moderate, surface flux spatial variability is high, and the weather regime is dominated by local processes. The adjusted continuous ranked probability score shows that the PBL uncertainty improves the skill of an ensemble subset in reproducing daily precipitation from the E-OBS observational product by 16%–20%. In comparison to WRF, WRF-Hydro improves this skill by 0.4%–0.7%. The reproduction of observed daily discharge with Nash–Sutcliffe model efficiency coefficients generally above 0.3 demonstrates the potential of WRF-Hydro in hydrological science.
Numerical atmospheric models generally consider terrestrial hydrological processes as only being vertical, in order to estimate the surface heat fluxes for constraining the atmospheric lower boundary condition. This is, for example, the case for the Weather Research and Forecasting (WRF) Model (Skamarock and Klemp 2008) coupled with the Noah land surface model (LSM; Chen and Dudhia 2001). In this approach, the lateral redistribution of soil moisture according to the topography and groundwater depth is, however, neglected. To relax this constraint and better represent soil moisture spatial variability, coupled atmospheric–hydrological models have been developed in recent years (e.g., Maxwell et al. 2007, 2011; Anyah et al. 2008; Gochis et al. 2015; Rahman et al. 2015; Wagner et al. 2016; Larsen et al. 2016).
The initiation and development of moist convection is sensitive to thermally induced wind circulations originating from the spatial variability in surface heat fluxes (e.g., Pielke 2001). Soil moisture heterogeneities can generate such circulations (e.g., Chen and Avissar 1994; Cheng and Cotton 2004; Taylor et al. 2007; Rieck et al. 2014). Coupling the Advanced Regional Prediction System (ARPS; Xue et al. 2000) with a three-dimensional and variably saturated groundwater flow model (ParFlow; Jones and Woodward 2001) for an idealized case study, Maxwell et al. (2007) found a shallow water table–induced circulation that impacted the location of convective cells after a 36-h run. Coupling WRF with the Hydrological Modeling System (HMS; Yu et al. 2006), Wagner et al. (2016) investigated groundwater effects on surface and atmospheric variables in a catchment of southeast China for an 8-yr period. Comparing WRF and WRF-HMS precipitation results, basin-averaged differences were minor, although spatial redistribution on the order of ±5% occurred. Rahman et al. (2015) simulated two convective events in western Germany with the Terrestrial System Modeling Platform (TerrSysMP; Shrestha et al. 2014), a version of the COSMO model coupled with ParFlow. They performed ensemble simulations based on perturbed initial conditions, with and without groundwater coupling. Their ensemble mean results supported the fact that groundwater dynamics noticeably affects soil moisture, surface fluxes, convective initiation, and the precipitation amounts. Recently, Larsen et al. (2016) brought evidence that including groundwater feedbacks in an atmospheric model can reduce the difference between simulated and observed seasonal precipitation, at least in the case of a river basin in Denmark.
Senatore et al. (2015) applied WRF and its hydrologically enhanced version, that is, WRF-Hydro (Gochis et al. 2015), to a catchment in southern Italy for a 3-yr period. Senatore et al. (2015) concluded that the lateral redistribution of soil moisture additionally resolved in WRF-Hydro reduced surface runoff and increased soil moisture amounts and drainage. However, the change in precipitation between WRF and WRF-Hydro was modest due to strong oceanic and orographic forcing in their study region. Differences between WRF and WRF-Hydro seasonal precipitation were also small in the case of a steep catchment at the foothills of Mount Kenya, East Africa (Kerandi et al. 2018). In West Africa, Arnault et al. (2016) found that the impact of overland flow and runoff–infiltration partitioning on precipitation was scale dependent, that is, much more noticeable in a 100 × 100 km2 domain, but not in a 500 × 2500 km2 domain. WRF-Hydro also produced daily discharge moderately close to observations according to the Nash–Sutcliffe model efficiency coefficient (NSE; Nash and Sutcliffe 1970): 0.27 in the case of Senatore et al. (2015), 0.43 in the case of Arnault et al. (2016), and 0.02 in the case of Kerandi et al. (2018).
The above studies show that the representation of terrestrial water flow indeed impacts the surface fluxes and planetary boundary layer (PBL) dynamics, thus potentially influencing precipitation. The representation of PBL dynamics in atmospheric models, through turbulence parameterization, also influences precipitation. However, it is unknown if the uncertainty in the representation of terrestrial water flow increases the precipitation spread originating from PBL scheme uncertainty.
This study addresses the precipitation sensitivity to the uncertainty in the representation of terrestrial water flow for central Europe during the warm season April–October 2008. An ensemble of WRF and WRF-Hydro simulations is performed based on various turbulence parameterization schemes and a varied runoff–infiltration partitioning parameter. The uncertainty in the representation of terrestrial water flow is accounted for by comparing simulations with a varied runoff–infiltration partitioning parameter and by comparing WRF and WRF-Hydro simulations. The first objective of this study is to assess the respective impacts of the uncertainties in turbulence parameterization and terrestrial water flow representation on precipitation. The second objective is to assess if these impacts depend on the weather regime (e.g., Barthlott et al. 2011; Keil et al. 2014) or the spatial variability in surface fluxes (e.g., Pielke 2001). The third objective is to assess the skill of the WRF and WRF-Hydro simulations in reproducing observed daily precipitation, as well as observed daily discharge for WRF-Hydro. The method to address these objectives, including the models used, the ensemble strategy, the observational validation dataset, and quantitative metrics, is detailed in section 2. Results are given in section 3, and a summary and perspective are finally provided in section 4.
a. Modeling approach: WRF and WRF-Hydro setups
The WRF Model (Skamarock and Klemp 2008) and the hydrologically enhanced version of WRF, that is, WRF-Hydro (Gochis et al. 2015), are used to simulate the regional land–atmosphere system over central Europe and investigate the impact of terrestrial water flow representation uncertainty on precipitation. The simulation period is from 1 January to 31 October 2008, with the first three months being considered for spinup time. Initial and lateral boundary conditions of the regional model are provided by the 6-hourly operational analyses (OP) at 0.125° resolution from the European Centre for Medium-Range Weather Forecasts (ECMWF).
In both WRF and WRF-Hydro setups the equations of atmospheric motion are solved at a time step of 10 s on a rotated grid of 0.025° (2.8 km) horizontal resolution centered over central Europe and slightly larger than that used in the Consortium for Small Scale Modeling (COSMO-DE) setup (Baldauf et al. 2011; Gebhardt et al. 2011) (see Fig. 1a). The topography is derived from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM; NASA 2015). The vertical coordinate is a terrain-following hybrid pressure coordinate (Skamarock and Klemp 2008), with 50 vertical levels and a pressure top at 10 hPa. Subgrid processes additionally parameterized are the longwave and shortwave radiative fluxes [Mlawer et al. (1997) and Dudhia (1989), respectively], cloud microphysics (Hong and Lim 2006), atmospheric turbulence, and surface heat and moisture fluxes as follows.
Three different PBL schemes are considered for the parameterization of turbulence: the Asymmetrical Convective Model version 2 (ACM2) scheme of Pleim (2007), the Mellor–Yamada–Janjić (MYJ) scheme of Janjić (1994), and the Yonsei University (YSU) scheme of Hong et al. (2006). This is to generate the model ensemble detailed in section 2b. These three particular PBL schemes have already been considered for evaluating turbulence parameterization uncertainty in WRF ensembles (e.g., García-Díez et al. 2013).
Surface fluxes are calculated with the Noah LSM predicting soil temperature and soil moisture in a 2-m-depth, four-layer column and taking into account vegetation effects (Chen and Dudhia 2001). Albedo, vegetation fraction, and leaf area index (LAI) are taken from satellite-derived climatology (Csiszar and Gutman 1999; Gutman and Ignatov 1998; Kumar et al. 2014), whereas other land surface and soil parameters are assigned for each land category and soil texture from the GlobCover2009 land cover map (Arino et al. 2012) and the Harmonized World Soil Database (HWSD; FAO/IIASA/ISRIC/ISS-CAS/JRC 2012), respectively. It is noted that the conversion table of Smiatek (2014) is used to translate GlobCover2009 and HWSD original classes into WRF Preprocessing System indexes.
The ratio between surface water WS and surface infiltration IS in the Noah LSM is computed at each model time step as
where Pd (m) is the precipitation falling on the bare soil, ΔZi (m) is the depth of soil layer i, θi (m3 m−3) is the volumetric water content (soil moisture) in soil layer i, θs (m3 m−3) is the saturated soil moisture (porosity) that depends on soil texture, Ks (m s−1) is the saturated hydraulic conductivity, Kref (m s−1) is the saturated hydraulic conductivity of the silty–clay–loam soil texture chosen as a reference, δt (s) is the model time step divided here by 86 400 s for a conversion in days, and k is the calibration parameter [k stands for kdtref in Chen and Dudhia (2001)]. In Eq. (1), the ratio WS/IS is decreased by increasing k (Schaake et al. 1996). The strength of the vertical water flow therefore depends on k. In the WRF case, WS stands for the surface runoff RS, which is the reason why k is referred as the runoff–infiltration partitioning parameter (Arnault et al. 2016).
Furthermore, the terrestrial water budget equation resolved in WRF can be written as
with precipitation P being equal to the sum of surface evaporation/sublimation E; surface and underground runoff being RS and RG, respectively; and a terrestrial water storage term ΔS that includes soil moisture, canopy water, and snow cover change. In this study, the terms of Eq. (2) are computed as water flux rates in millimeters per day.
In the WRF-Hydro setup, which also includes the Noah LSM as in the WRF setup, the atmospheric grid at 2.8-km resolution (Fig. 1a) is coupled with a terrestrial subgrid at 280-m resolution including a river network (Fig. 1b) in order to route surface water WS overland, soil moisture in the subsurface, and stream water in the river channels (Gochis et al. 2015). The time step to resolve these terrestrial processes is set to 10 s, as for the atmospheric processes. Technically, WS and soil moisture are conservatively disaggregated on the terrestrial subgrid, routed, and aggregated back to the atmospheric grid each time step. The topography of the terrestrial subgrid is derived from the ASTER dataset. The river network is obtained with ArcGIS software using the Catchment Characterization and Modeling version 2.1 (CCM2) database (de Jager and Vogt 2010). It is emphasized that this WRF-Hydro setup allows for reinfiltration, that is, WS infiltrating at a later time step eventually at a different grid point, and exfiltration, that is, WS originating from water excess in a fully saturated soil column. The WS reaching a river channel grid point in the terrestrial subgrid is finally removed from the land and added to the water in the river channel. The part of WS contributing to the river flow stands for the surface runoff RS [Eq. (2)] in the WRF-Hydro setup and is added to the model outputs (as in Arnault et al. 2016). This allows us to define a WRF-Hydro-derived terrestrial water budget equation as in Eq. (2), in which the storage term ΔS also includes the surface water change ΔWS.
In the WRF-Hydro setup, the river water volume is routed on a pixel-by-pixel basis using a diffusive wave formulation allowing for backwater effects. Channel parameters, including the initial river head, bottom width, and side slope of the river channel, and Manning’s channel roughness coefficient are prescribed as functions of Strahler stream order (Strahler 1957). In comparison to other WRF-Hydro studies (e.g., Yucel et al. 2015; Arnault et al. 2016; Kerandi et al. 2018), the Manning coefficients used here have been decreased for calibration purposes, ranging from 0.35 at the order 1 to 0.03 at the order 8. Otherwise, the further channel parameters are default. The WRF-Hydro model could also be further calibrated with respect to discharge, by tuning land surface and soil parameters such as, among others, the stomatal resistance, surface roughness, and the soil’s hydraulic conductivity, which are currently assigned as a function of land category and soil texture. However, this study focuses only on the runoff–infiltration partitioning parameter as it plays the key role for surface water retention and subsequent redistribution by the lateral routing components of WRF-Hydro. Investigating the impact of further parameters would be beyond the scope of this article.
Moreover, the baseflow contribution to the river flow in each basin is evaluated with a basin-attributed linear drainage bucket model, using the pass-through option (Gochis et al. 2015). This means that the underground runoff RG generated in a basin area is collected and directly redistributed to all the channel grid cells of the basin. The following equation of specific discharge Q generation at basin scale is derived as
where RS and RG have already been defined above, and ΔWR is the river water storage term in the basin. It is noted that the specific discharge is computed as the ratio between discharge (m3 s−1) and basin area (m2). As for Eq. (2), terms of Eq. (3) are also computed in millimeters per day.
Surface variables, including those in the budget equation [Eq. (2)], are stored at an hourly interval, as well as the geopotential heights at 500 hPa and the convective available potential energy (CAPE). For WRF-Hydro, terms of Eq. (3) are additionally saved on the terrestrial subgrid at a daily interval.
b. Ensemble strategy
There is a potential link between terrestrial water flow, soil moisture, surface fluxes, PBL dynamics, and precipitation (e.g., Chen and Avissar 1994; Cheng and Cotton 2004; Taylor et al. 2007; Maxwell et al. 2007; Rieck et al. 2014; Senatore et al. 2015; Rahman et al. 2015; Arnault et al. 2016; Wagner et al. 2016; Larsen et al. 2016; Kerandi et al. 2018). A comparison between the respective effects of terrestrial water flow and PBL dynamics uncertainty could provide further insight on the precipitation sensitivity to terrestrial water flow.
Accordingly, the following WRF/WRF-Hydro ensemble is constructed in order to allow for a comparison between turbulence and terrestrial water flow uncertainty effects. The WRF and WRF-Hydro setups of section 2a are run for the three PBL schemes, ACM2, MYJ, and YSU, and for two values of k [Eq. (1)], that is, k = 1 and k = 3, with 3 being the default value. The two values for k represent the terrestrial water flow uncertainty in the vertical, the two models, WRF and WRF-Hydro, represent the contribution of lateral flow to terrestrial water flow uncertainty, and the three PBL schemes represent the turbulence parameterization uncertainty. This makes an ensemble (ENS) of 12 members: six WRF members and six WRF-Hydro members. Figure 2 provides a conceptual view of ENS. The ensemble subsets of the members using the PBL scheme X are called ENS(X), where X stands for ACM2, MYJ, or YSU. The ensemble subsets of the members using k = 1 and k = 3 are called ENS(k = 1) and ENS(k = 3), respectively, whereas those of the WRF and WRF-Hydro members are called ENS(WRF) and ENS(Hydro), respectively.
A so-called control ensemble, similar to ENS except for the initial time, is generated, in order to evaluate the magnitude of random noise in the model results. Here, we arbitrarily choose to initialize this control ensemble on 2 January 2008.
c. Validation datasets
Precipitation ensemble results are validated with the E-OBS gridded precipitation product PEOBS from the European Climate Assessment and Dataset project (Haylock et al. 2008). The PEOBS product is available daily on a grid at 0.25° resolution. As remarked by Haylock et al. (2008), PEOBS is the product of an interpolation method that has been designed to facilitate the comparison with regional climate models at the same spatial scale.
Daily discharge QGRDC from the Global Runoff Data Center (GRDC 2013) is also considered in this study for validating the results from the WRF-Hydro hydrological extension. The selected stations are Koeln, Liebenau, Neu Darchau, Kienstock, and Passau, located at the rivers Rhine, Weser, Elbe, Danube, and Inn, respectively (see Fig. 1b). The covered basin areas are referred to as Rhine–Koeln (144 × 103 km2), Weser–Liebenau (20 × 103 km2), Elbe–Neu Darchau (131 × 103 km2), Danube–Kienstock (96 × 103 km2), and Inn–Passau (28 × 103 km2), respectively. The total area covered by these five river basins has a size of about 625 km × 625 km and is named area A. This set of river basins has been chosen as it covers a large part of central Europe, that is, area A (Fig. 1b), and allows us to test WRF-Hydro for different catchment sizes in various hydrological environments: a mixed high–low mountainous region (Rhine–Koeln, Danube–Kienstock, Inn–Passau) and low mountainous region (Weser–Liebenau, Elbe–Neu Darchau).
d. Continuous ranked probability score
The skill of an ensemble subset ENS(X) in reproducing daily precipitation PEOBS is assessed with the continuous ranked probability score (CRPS; e.g., Matheson and Winkler 1976; Hersbach 2000; Gneiting and Raftery 2007). The CRPS has been selected for this skill assessment as it is a strictly proper scoring rule, following decision-theoretically justified principles (e.g., Gneiting and Raftery 2007). A high or low skill of the ensemble subset ENS(X) is indicated by a low or high CRPSENS(X), respectively:
where FENS(X) and FEOBS stand for the cumulative distribution functions of modeled and observed daily precipitation P, respectively.
It is noted that the unit of the CRPS is identical with that of daily precipitation, which in our case is millimeters per day. One could decide to normalize the CRPS [Eq. (4)] by the observed daily precipitation in order to have a dimensionless score. However, such a normalized CRPS would be improper as it would favor forecast distributions that unduly emphasize low precipitation events, and therefore become misleading (Lerch et al. 2017).
The integral in Eq. (4) can be evaluated via
where NENS(X) is the number of members in ensemble subset ENS(X), and PL and PM are the precipitation from members L and M, respectively.
However, as noted by Ferro et al. (2008), the CRPS computed as such favors ensembles with a larger number of members. Since our interest focuses on the relative comparison of the skill of ensemble subsets of different size, we follow Ferro et al. (2008) and Fricker et al. (2013), who proposed the following “unbiased” expression for the CRPS, the so-called fair CRPS:
Equation (6) is asymptotically equivalent to Eq. (5) for an infinite number of members NENS(X) and allows for a fair comparison of differently sized ensembles (Ferro et al. 2008). In the following, the CRPS of each ensemble subset of section 2b is evaluated with Eq. (6) as horizontal maps and spatially averaged in area A (see Fig. 1b).
e. Normalized ensemble spread
The respective effects of terrestrial water flow and turbulence uncertainty are evaluated with the normalized ensemble spread S (Hohenegger et al. 2006; Keil et al. 2014). Parameter S is computed for the ensemble ENS as
where PENS is the ensemble mean daily precipitation, and SENS [Eq. (7)] is displayed as horizontal maps and computed as area averages on grid points receiving more than 1 mm day−1 in area A.
Parameter S is also evaluated for groups of ensemble subsets as
where G is a group of NG ensemble subsets (SUB) of size NSUB, PSUB is the ensemble mean daily precipitation from SUB, and SG is the normalized ensemble spread relative to G. Three groups of ensemble subsets are considered in this study, as illustrated in Fig. 2: 1) the four ensemble subsets in which only the PBL scheme is varied, 2) the six ensemble subsets in which only the value of k is varied, and 3) the six ensemble subsets in which only the model, that is, WRF or WRF-Hydro, is varied. Note that each group comprises 12 individual simulations. These three groups aim at quantifying the precipitation sensitivity to 1) turbulence parameterization uncertainty, 2) terrestrial water flow uncertainty in the vertical direction, and 3) terrestrial water flow uncertainty originating from the consideration of resolved lateral flow. Their associated normalized precipitation spreads are named SPBL, Sk, and SHydro. The impact of these respective effects can finally be assessed by comparison to SENS.
f. Weather regime dependence
In regional atmospheric modeling, internal processes uncertainty preferentially affects precipitation during weak synoptic forcing episodes (e.g., Stensrud et al. 2000; Keil et al. 2014). The question here is, if the precipitation sensitivity to the uncertainty in the representation of terrestrial water flow also depends on the level of synoptic forcing.
Keil et al. (2014) used the convective adjustment time scale originally proposed by Done et al. (2006) and modified by Keil and Craig (2011) and Zimmer et al. (2011), in order to objectively determine the level of synoptic forcing and associated weather regime:
where is proportional to the ratio between the CAPE (J kg−1) and precipitation P. Other parameters in Eq. (9) are the density reference ρ0 (kg m−3), the specific heat of air at constant pressure cp (J kg−1 K−1), the temperature reference T0 (K), the latent heat of vaporization Lυ (J kg−1), and the terrestrial gravitational acceleration g (m s−2). Choosing millimeters per hour for the precipitation unit, which is equivalent to kilograms per square meter per hour as water density is 1000 kg m−3, the resulting [Eq. (9)] comes in hours.
Parameter provides a measure of how fast conditional instability is removed in an atmospheric column through the release of moist convection. If is much smaller than the time scale characterizing the development of the synoptic environment, convection is in equilibrium with the synoptic-scale forcing and hence controlled by it. In the case of weak synoptically forced situations, CAPE can build up as local forcing is generally not as efficient as synoptic forcing in triggering precipitation. In this case, area-averaged is larger than 3–6 h (Keil et al. 2014; Kühnlein et al. 2014).
In this study the weather regime dependence is tested for the contribution of the terrestrial water flow representation uncertainty to the normalized ensemble spread of daily precipitation. As in Keil et al. (2014), hourly CAPE and P are convolved with a Gaussian kernel of half-width size of 56 km before computing an hourly value for [Eq. (9)], which ensures that is characteristic of an enlarged environment around the precipitating systems. Parameter is then spatially averaged in area A for pixels receiving more than 1 mm h−1. Following Kühnlein et al. (2014), daily values of are computed as the hourly maxima reached in a day. If is larger than 6 h, the weather regime is governed by synoptic processes. Finally, the ensemble mean of these daily values, called , is compared to SPBL, Sk, and SHydro (see definition in section 2e).
g. Surface flux spatial variability dependence
Convective precipitation is expected to be sensitive to the wind circulation induced by surface flux spatial variability or heterogeneity H (e.g., Chen and Avissar 1994; Pielke 2001; Cheng and Cotton 2004; Taylor et al. 2007; Rieck et al. 2014; Rahman et al. 2015). It is therefore desirable to objectively define H and see if it relates to the normalized precipitation spreads, and if the precipitation sensitivity to the uncertainty in the representation of terrestrial water flow depends on H. It is chosen here to compute H as the norm of the gradient of surface evaporation/sublimation E:
where symbols and ∥ ∥ stand for the gradient and norm operator. Parameter E is convolved with a 5-pixel/14-km diameter circular mean filter kernel, before computing the gradient of E, as boundary layer motions are known to be more sensitive to surface heterogeneities of the neighborhood size 10–20 km (e.g., Clark et al. 2004). Parameter is then convolved with a Gaussian kernel of half-width size of 56 km, as for CAPE and P in the computation of (section 2f).
As for , H is evaluated at the hourly time scale and spatially averaged for the areas with hourly precipitation exceeding a threshold of 1 mm h−1 within area A. Surface flux spatial variability has a potential impact on daily precipitation uncertainty through its effect on PBL processes during the initiation of moist convection (e.g., Pielke 2001). Consequently, the value of H at the hour of the day when reaches its maximum, that is, the initial stage of convection, is selected for computing a daily ensemble mean HENS, which is then compared to SPBL, Sk, and SHydro.
The fact that HENS is computed in a similar manner as (section 2f) ensures that both quantities are representative of the same environment. Parameter HENS is finally evaluated in millimeters per day per kilometer, as E is computed as a water flux rate in millimeters per day [see section 2a, Eq. (2)] and the grid spacing for the spatial derivate is provided in kilometers.
a. Model validation
1) Synoptic-scale dynamics
Synoptic dynamics is assessed here with the geopotential heights at 500 hPa called Z500, as usually done in midlatitude meteorology (e.g., Holton 2004). Figure 3a displays the 6-hourly evolution of the area A–averaged Z500 from the ECMWF OP, that is, Z500ECMWF, and from the ENS members, that is, Z500ENS. It shows that Z500ENS remains close to Z500ECMWF, and that there is not much variation between ENS members. This is confirmed by the 6-hourly evolution of the spatial root-mean-square error (RMSE) between these two quantities, with values between 5 and 30 m (Fig. 3b). The small spread of Z500ENS at the 6-h time scale suggests that the simulation domain (Fig. 1a) is small enough that the boundary conditions effectively control the synoptic-scale dynamics of the solution in the interior of the domain.
2) Mean precipitation
The mean precipitation of the ensemble PENS, as well as that of the ensemble subsets PENS(ACM2), PENS(MYJ), PENS(YSU), PENS(k=1), and PENS(Hydro), are averaged for the simulation period April–October 2008 and are compared with the observation PEOBS in Fig. 4. Averaged PEOBS is above 3 mm day−1 in the southern high mountainous region of area A, and below 3 mm day−1 in the lower regions to the north. This spatial distribution is qualitatively reproduced by PENS (cf. Figs. 4a and 4b), although with a mean overestimation of 23%. This overestimation is mainly distributed in the southeastern half of area A, and particularly in the upper Elbe River basin, where the difference-to-observation value reaches 100% (Fig. 4c).
YSU mainly increases the average precipitation, while ACM2 mainly decreases it and MYJ has a mixed effect (Figs. 4d–f). However, these PBL effects are much smaller than the difference to PEOBS (cf. color scales in Figs. 4c and Figs. 4d–f). The change in average precipitation for the period April–October 2008 induced by modifications in the representation of terrestrial water flow ranges from −5% to +5% (Figs. 4g,h), which is comparable to that obtained by Wagner et al. (2016). However, this is smaller than the change induced by the PBL scheme, that is, ranging from −20% to +20%.
3) Daily precipitation
The skill of the ENS members in reproducing daily PEOBS is assessed with the continuous ranked probability score CRPSENS defined in section 2d [Eq. (6)] and displayed in Fig. 5a. In average for the period April–October 2008, CRPSENS is above 2 mm day−1 in the southern high mountainous region of area A and below 2 mm day−1 in the lower regions to the north. This pattern is similar to that of the ensemble mean precipitation PENS (cf. Figs. 4b and 5a), so that the lowest ensemble skill is associated with the highest precipitation, mainly in the high mountainous region. On average over the entire area A, CRPSENS is about 1.5 mm day−1 (Table 1).
The ability of each PBL scheme to improve the ensemble skill is investigated and is shown in Figs. 5b–d, displaying the respective relative differences between CRPSENS(ACM2), CRPSENS(MYJ), CRPSENS(YSU), and CRPSENS. All three of the PBL ensemble subsets have a CRPS that is substantially higher than CRPSENS. This indicates that even the best PBL subensemble is inferior to the ensemble containing all three PBL schemes in terms of probabilistic skill. Overconfidence is known to penalize the probabilistic skill of an ensemble (e.g., Weigel et al. 2008). Accordingly, each PBL ensemble subset appears to be overconfident, and the full ensemble reduces this overconfidence.
The impact of the other ensemble subsets on the CRPS, that is, ENS(WRF, k = 1), ENS(WRF, k = 3), ENS(Hydro, k = 1), ENS(Hydro, k = 3), is investigated in Figs. 5e–h. These panels show that each of these ensemble subsets has a CRPS slightly smaller than that of ENS. This means that an ensemble subset considering the three PBL schemes, but only one value for k, and only one model, that is, WRF or WRF-Hydro, has a better skill than the full ensemble itself. Accordingly, following Weigel et al. (2008), considering only one option for the representation of terrestrial water flow reduces the overconfidence of the ensemble. Table 1 further shows that the ensemble subset using WRF-Hydro and k = 1, that is, ENS(Hydro, k = 1), exhibits the lowest CRPS, with an induced area-average diminution of −5.5%.
The robustness of this result is assessed in a bootstrap approach by computing the CRPS based on 100 randomly selected days instead of the full period April–October 2008. Reconducting this random evaluation one thousand times, it is found that ENS(Hydro, k = 1) exhibits the lowest mean CRPS in 90% of the cases. This suggests that the skill of a WRF ensemble can be improved by considering several atmospheric turbulence parameterization options, but only one best option for the terrestrial water flow representation, which in this case is WRF-Hydro with k = 1.
4) Basin-averaged results
Basin-averaged daily time series of PENS in the river basins Rhine–Koeln, Weser–Liebenau, Elbe–Neu Darchau, Danube–Kienstock, and Inn–Passau are close to those of PEOBS (Fig. 6), with a correlation coefficient r between 0.67 and 0.91 and a mean difference-to-observation value between +12% and +39% (Table 2). This is in agreement with the precipitation overestimation discussed in section 3a(2) (see Fig. 4c). The timing of the daily peaks is generally well captured, and there is not much difference in terms of daily variation among ENS members’ basin-averaged precipitation time series (see range of the red lines in Fig. 6).
WRF-Hydro produces daily time series of discharge QENS(Hydro) moderately close to that observed at the outlets of the five selected river basins (Fig. 7), with a mean difference to observations between −39% and +17% and an NSE value ( Nash and Sutcliffe 1970) generally above 0.3 (Table 2). This is comparable to what has previously been obtained with WRF-Hydro (e.g., Yucel et al. 2015; Senatore et al. 2015; Arnault et al. 2016; Kerandi et al. 2018), but noticeably lower to what is usually obtained with traditional “uncoupled” hydrological models in part because the meteorological forcing is prescribed in this latter case (e.g., Newman et al. 2015; Zink et al. 2017). Indeed, the quality of the WRF-Hydro discharge not only depends on the simulated amount of basin-averaged precipitation (Fig. 6), but also on the spatial distribution of simulated precipitation within the basin, which is challenging for an atmospheric model. It is also acknowledged that the discrepancies between simulated and observed discharge in Fig. 7 may also be related to 1) distributed parameters not properly calibrated for this particular application, such as the unsaturated hydraulic conductivity, surface roughness or stomatal resistance, and 2) the relatively short 3-month spinup time considered in this study for snow and soil moisture in deeper soil layers to reach equilibrium. Still, it is of high interest to assess how well the coupled WRF-Hydro modeling system is able to reproduce observed streamflow at respective gauge locations.
The best discharge result is obtained for Elbe–Neu Darchau, that is, NSE > 0.8, with a peak in late April 2008 probably related to snowmelt, and a slow descent afterward (Fig. 7c). The Weser–Liebenau’s discharge has a similar temporal variation, although in this case QENS(Hydro) largely underestimates the baseflow contribution from June to October 2008 (Fig. 7b). NSE is accordingly smaller, with values above 0.39. This could be related to a relatively large contribution of soil water older than 3 months in the case of the Weser–Liebenau discharge, which would require a longer spinup time in order to be properly taken into account. Danube–Kienstock and Inn–Passau’s discharge displays many more peaks, which makes hydrological modeling in these basins more challenging. Parts of these two river basins are located in a high mountainous region where fast surface runoff generation during a storm event is more likely to occur. Peaks simulated in QENS(Hydro) are still relatively close to those in QGRDC (see Figs. 7d,e), as deduced by a NSE above 0.3. Finally, concerning the discharge of Rhine–Koeln, QENS(Hydro) displays such peaks as well, in relation to an upper catchment area also in the high mountainous region. However, these peaks are not in QGRDC, with this last one also displaying a higher base flow (see Fig. 7a). These discrepancies are confirmed by relatively low NSEs, between 0.0 and 0.39 (Table 2), and could be related to the fact that the storing influence of lakes is not taken into account in WRF-Hydro, as it is known that Lake Constance and the lakes of the Swiss Alps in the Upper Rhine attenuate the flood peaks in the Lower Rhine (e.g., Lohre et al. 2003; Bronstert et al. 2007). Such a smoothing is not seen in the case of the Danube–Kienstock and Inn–Passau discharge, which is potentially explained by the fact that these two basins do not have a lake as big as Lake Constance.
b. Relative impact of terrestrial water flow representation on land surface variables
In coupled hydrological models, the consideration of terrestrial water flow primarily modifies the distribution of soil moisture and surface fluxes (Maxwell et al. 2007; Senatore et al. 2015; Rahman et al. 2015; Arnault et al. 2016; Wagner et al. 2016; Larsen et al. 2016; Kerandi et al. 2018). In the following, this impact is discussed with the terrestrial water budget [Eq. (2)] and river water budget [Eq. (3)]. Figure 8a displays daily time series of the terms of Eq. (2) spatially averaged in area A (see Fig. 1b) and derived from the mean of ENS members. It shows that ΔS and RS mainly follow the daily variations of P, whereas E and RG have a much smoother temporal variation. The differential daily time series between the mean of ENS(k = 1) [ENS(Hydro)] and ENS members are displayed in Fig. 8b (Fig. 8c) and summed in Table 3. Figures 8b and 8c show that daily changes in ΔS mainly counterbalance those in RS. Parameters E and RG are not much affected on the daily scale, but noticeable changes in the total sums are found (see Table 3). The increase of ΔS is associated with an increase of E, RG, and P, whereas the decrease of ΔS is associated with a decrease of E, RG, and P. However, the impact of modified terrestrial water flow representation on the average precipitation in area A is small: ±2 mm over 631 mm for the considered time period, which is about ±0.3% (Table 3).
The robustness of the above result is assessed by comparing with the control ensemble result. In the control ensemble the total sum of precipitation is increased by about 0.25%, which gives the magnitude of random noise. The change in precipitation induced by modifying the representation of terrestrial water flow is slightly larger than this random noise. According to a standard t test, the change is significant only at the p > 0.1 level. This confirms that the terrestrial water flow uncertainty has a reduced and barely significant impact on the total sum of precipitation in area A. In the control ensemble, we find that modifying the representation of terrestrial water flows induces a change in ΔS, E, RG, and P of the same sign and with a similar amplitude as in Table 3 (not shown). In terms of physical processes, this suggests that some of the water evaporating from the surface in area A falls back in area A as precipitation, so that changing the amount of surface evaporation by decreasing surface infiltration or enabling lateral terrestrial water flow would have a direct effect on the amount of precipitation.
Figure 9a displays daily time series of the terms of Eq. (3) spatially averaged in area A (see Fig. 1b) and derived from the mean of ENS(Hydro) members. It shows the contributions of RS and RG to daily discharge generation in area A. Isolated peaks in simulated Q are mainly related to RS. The river water storage term ΔWR acts as a buffer and smooths the Q isolated peaks with respect to those of RS. This smoothing can be further calibrated with the Manning roughness coefficients, which are currently specified as a function of stream order.
Decreasing k to 1 slightly increases RS, although the associated increase in Q is partly counterbalanced by a smaller decrease in RG (Fig. 9b, Table 4). However, the impact of decreasing k on the area-average discharge is small: +2 mm over 182 mm for the considered time period, which is about +1.1% (Table 4). As a side note, this is a noticeable difference with West Africa, where Arnault et al. (2016) found that the amount of simulated discharge was tightly related to the value of k.
c. Relative impact of terrestrial water flow representation on precipitation
The impact of terrestrial water flow representation is further discussed with the normalized ensemble spreads of daily precipitation SENS, SPBL, Sk, and SHydro (see section 2e), displayed in Fig. 10 as averaged maps for the period April–October 2008. Parameter SENS is generally between 0.5 and 0.8 in the central and northern parts of area A but is below 0.5 in the southern part in association with a strong orographic forcing on precipitation there.
Parameter SPBL is generally slightly lower than SENS (Fig. 10b), with an area-average difference of about 5%. In comparison, Sk and SHydro are much lower (see Figs. 10c,d), with an area-average difference of about 40%–50%. This shows that the uncertainty in the modeled atmospheric turbulence is largely dominating the full ensemble spread SENS and that the averaged contribution of the terrestrial water flows representation uncertainty to this variability is 5%. The reduced spread between the ensemble members using the same PBL scheme is in accordance with the fact each of the PBL ensemble subsets ENS(ACM2), ENS(MYJ), and ENS(YSU) is overconfident in comparison to the full ensemble ENS and that ENS is overconfident in comparison to each of the ensemble subsets that consider only one terrestrial water flow option, as discussed in section 3a(3).
Parameter SENS is much closer to SPBL in the southern part of area A, suggesting that a strong orographic forcing inhibits the impact of terrestrial water flow representation uncertainty on precipitation. The largest impact occurs in the central and northern parts, where SENS locally exceeds SPBL by 5%–20% (Fig. 10b). Accordingly, the uncertainty in the representation of terrestrial water flow preferentially affects precipitation in regions with moderate topography. The next sections investigate if the sensitivity of precipitation to terrestrial water flow uncertainty varies from day to day, potentially in relation with the weather regime and surface flux spatial variability.
d. Role of weather regime
The weather regime dependence of SPBL, Sk, and SHydro is investigated with the daily convective adjustment time scale averaged for the ENS members in area A, referred as (section 2f). Figure 11a displays the scatterplot between daily values of and SPBL. There is a positive correlation between SPBL and , with a coefficient of determination R2 of 0.53. A similar correlation is obtained for Sk and SHydro (not shown). This is much higher than that obtained in Keil et al. (2014) for forecast ensembles of hourly precipitation, possibly because the analysis here is conducted with daily values. This confirms that internal processes uncertainty preferentially affects precipitation during weak synoptic forcing episodes (e.g., Stensrud et al. 2000; Keil et al. 2014).
The red bold plus signs in Fig. 11 indicate the days when area-average SENS exceeds area-average SPBL by more than 20%. These are the days when the precipitation sensitivity to the uncertainty in the representation of terrestrial water flow is most noticeable. According to the red bold plus signs in these scatterplots, these days are mostly associated with a weak synoptically forced weather regime ( higher than 6 h).
e. Role of surface flux spatial variability
The surface flux spatial variability dependence of SPBL, Sk, and SHydro is investigated with the surface flux spatial heterogeneity averaged for the ENS members, referred as HENS (section 2g). Figure 11b displays the scatterplot between daily values of HENS and SPBL. The correlation is close to that obtained in previous section, with a R2 of 0.40. A similar correlation is obtained for Sk and SHydro (not shown). This means that internal processes uncertainty preferentially affects precipitation when surface flux spatial variability is high. The fact that the coefficient of determination in Fig. 11a is close to that in Fig. 11b further suggests that surface flux spatial variability is as important as the weather regime to explain the precipitation ensemble spread driven by model physics.
As in Fig. 11a, the red bold plus signs in Fig. 11b indicate the days when area-average SENS exceeds area-average SPBL by more than 20%. It is noted that these days are associated with relatively large values of HENS, that is, above 0.15 mm day−1 km−1, whereas the total range of HENS values spreads from 0 to 0.4 mm day−1 km−1. This result suggests that, in central Europe, the uncertainty in the representation of terrestrial water flow most noticeably affects precipitation preferentially when surface flux spatial variability is high.
4. Summary and perspective
This study addressed the precipitation sensitivity to the uncertainty in the representation of terrestrial water flow in central Europe. For this purpose, an ensemble of WRF and WRF-Hydro simulations was generated for the period April–October 2008, using the COSMO-DE grid at 2.8-km resolution and ECMWF operational analyses as forcing data, taking into account turbulence uncertainty with three PBL parameterization (ACM2, MYJ, or YSU schemes), and vertical water flow representation uncertainty with two values for the runoff–infiltration partitioning parameter k, that is, 1 and 3. The contribution of lateral flow to the uncertainty in the representation of terrestrial water flow was considered by comparing the WRF and WRF-Hydro simulations.
Enabling lateral terrestrial water flow in WRF-Hydro generally increased soil moisture, surface evaporation, and precipitation in the study region, whereas these variables were generally decreased by reducing the runoff–infiltration partitioning parameter. However, a similar precipitation change was found by changing the initial time of the ensemble, which confirmed that the total sum of precipitation in the study region was not much affected by terrestrial water flow uncertainty. The fact that the change in surface evaporation and precipitation were in phase, both in the ensemble and in the control ensemble, is thought to be related to regional precipitation recycling.
The impact of these terrestrial water flow effects on precipitation was evaluated with the normalized ensemble spread. The full ensemble spread was largely dominated by the uncertainty in the modeled atmospheric turbulence. On average for April–October 2008, the difference between the full and turbulence-uncertainty-driven ensemble spread was about 5%, which was considered to be the averaged effect of terrestrial water flow representation uncertainty. This averaged effect was found to be inhibited over steep terrain, but enhanced in regions where topography is moderate. Special attention was drawn to the days when the normalized ensemble spread of the full ensemble was about 20% higher than that originating from turbulence parameterization uncertainty. These particular days were associated with relatively high values of the convective adjustment time scale and surface flux spatial variability. It was therefore concluded that the uncertainty in the representation of terrestrial water flow is more likely to noticeably affect precipitation when the weather regime is weakly synoptically forced and surface flux spatial variability is high. As 1) the additional description of lateral terrestrial water flow noticeably increases the normalized ensemble spread in certain conditions and 2) the WRF-Hydro coupled modeling system allows us to describe this lateral flow at a moderate computational cost, we argue that WRF-Hydro is suitable for ensemble forecasting.
Further focus of the study was on assessing the skill of ensemble members in reproducing the E-OBS daily precipitation product. Spatially averaged daily time series of E-OBS precipitation in five large river basins in central Europe were approximately well reproduced, with a correlation coefficient above 0.67. Precipitation was on average overestimated between +12% and +39%. As for the normalized ensemble spread, modifying the representation of terrestrial water flow had much less impact on the average precipitation amount than the choice of the PBL scheme. The ensemble skill in reproducing E-OBS daily precipitation was further evaluated with the CRPS. The ensemble subset considering the three PBL schemes, but only the WRF-Hydro model and the value of 1 for k provided the lowest CRPS (highest ensemble skill). The decrease of probabilistic skill, or increase of confidence, for the ensemble subsets considering several terrestrial water flow options is coherent with the reduced ensemble spread induced by the terrestrial water flow uncertainty. Future studies should investigate if for other years and/or other regions the skill of a WRF ensemble is also improved by considering several atmospheric turbulence parameterization options, but only one best option for the representation of terrestrial water flow.
Finally, the WRF-Hydro ensemble driven by ECMWF operational analyses showed diverse skills in reproducing streamflow in large European river basins, with NSE ranging from 0 to 0.91. Future research with WRF-Hydro could also focus on an enhanced calibration of the various distributed land surface and soil parameters and the use of multiyear spinup to investigate the potential of further improvement in reproducing streamflow.
This research is part of the subproject “A5-The role of soil moisture and surface- and subsurface water flows on predictability of convection” of the Transregional Collaborative Research Center SFB/TRR 165 “Waves To Weather” funded by the German Science Foundation (DFG). ECMWF operational analyses were obtained within the framework of a special project at the ECMWF. WRF and WRF-Hydro simulations were performed on the computational resource ForHLR I funded by the Ministry of Science, Research and the Arts Baden-Württemberg and DFG (Deutsche Forschungsgemeinschaft). We acknowledge the P-EOBS dataset from the ECA&D project (http://www.ecad.eu) and the Q-GRDC dataset from the Global Runoff Data Center (http://www.bafg.de/GRDC/EN/Home/homepage_node.html). Special thanks go to Christian Barthlott for providing COSMO static input data; George Craig, Stephan Rasp, Michael Riemer, and Lotte Bierdel for fruitful discussions about skill scores; Gerhard Smiatek, Dominikus Heinzeller, Christof Lorenz, Patrick Laux, Thomas de Couet, Carsten Maas, Hartmut Häfner, and Christoph Sörgel for the logistic and computer support; and two reviewers for helping to clarify the manuscript.
Current affiliation: Fraunhofer Institute for Industrial Engineering IAO, Stuttgart, Germany.
This article is included in the Waves to Weather (W2W) Special Collection.