The spatial variability of topsoil water content (SWC) is often expressed through the relationship between its spatial mean 〈θ〉 and standard deviation σθ. The present study tests the concept that a reasonably performing land surface model (LSM) should be able to produce σθ–〈θ〉 data pairs that fall into a polygon, spanned by the cloud of observed data and two anchor points: σθ at the permanent wilting point σθ–〈θwp〉 and σθ at saturation σθ–〈θs〉. A state-of-the-art LSM, Noah-MP, was driven by atmospheric forcing data obtained from eddy covariance field measurements in two regions of southwestern Germany, Kraichgau (KR) and Swabian Alb (SA). KR is characterized with deep loess soils, whereas the soils in SA are shallow, clayey, and stony. The simulations series were compared with SWC data from soil moisture networks operating in the two study regions. The results demonstrate that Noah-MP matches temporal 〈θ〉 dynamics fairly well in KR, but performs poorly in SA. The best match is achieved with the van Genuchten–Mualem representation of soil hydraulic functions and site-specific rainfall, soil texture, green vegetation fraction (GVF) and leaf area index (LAI) input data. Nevertheless, most of the simulated σθ–〈θ〉 pairs are located outside the envelope of measurements and below the lower bound, which shows that the model smooths spatial SWC variability. This can be mainly attributed to missing topography and terrain information and inadequate representation of spatial variability of soil texture and hydraulic parameters, as well as the model assumption of a uniform root distribution.
Soil moisture, next to sea surface temperature, is the second-most important factor determining the predictability of the atmosphere’s state (Dirmeyer and Shukla 1994; Dirmeyer 1995). Soil water content (SWC), in particular that of the topsoil, interacts with atmospheric temperature, clouds, and precipitation formation over land. It also influences surface albedo, boundary layer evolution, and energy partitioning into sensible and latent heat fluxes (Yeh et al. 1984; Goodrich et al. 1994; Heathman et al. 2009; Seneviratne et al. 2010; Heathman et al. 2012). SWC is an important variable at the lower boundary of large-eddy simulation (LES) and numerical atmospheric models. Its realistic representation is essential, because it is a key biophysical and hydrologic state variable, used in a range of practical applications including irrigation scheduling, agriculture, turbulent flux simulations, quantitative rainfall forecasting, climate simulation, and weather prediction (e.g., Morari and Giardini 2002; Leib et al. 2003; Wraith et al. 2005).
Currently, individual atmospheric models provide uncertainties and biases in simulated SWC (Koster et al. 2004; Dirmeyer et al. 2006). Hauck et al. (2011) showed that biases in soil moisture significantly influence simulated precipitation. Gantner and Kalthoff (2010) found that spatial variability of the top 0.1-m SWC impacted the formation of precipitation in West African conditions—an area of strong land surface–atmosphere coupling (Koster et al. 2004). Simulations showed that inhomogeneous soil moisture conditions induced more intense rainfalls than homogeneous conditions. When a mature convective system moved over the study area, the precipitation pattern was significantly reduced over drier areas. This calls for further model improvements and validation against observations.
During the past five decades, land surface models (LSMs) have been continuously evolving according to the requirements of atmospheric and hydrological disciplines (Chen and Dudhia 2001; Seneviratne et al. 2010; Niu et al. 2011). Noah is one of the commonly used land surface models for hydrological processes (Chen and Dudhia 2001). Using the dataset of the First International Satellite Land Surface Climatology Project (ISLSCP) Field Experiment (FIFE; Betts and Ball 1998), the test of the Noah LSM with multiparameterization options (Noah-MP) showed that simulated daily soil water storage over 1 m matched observations well (Niu et al. 2011). The study period lasted 2 months during summer 1987. The data were averaged over 20 grassland sites within a domain area of 15 km × 15 km. Cai et al. (2014) reported fairly low root-mean-square errors (RMSEs) against monthly averaged observations of 60 Soil Climate Analysis Network (SCAN; Schaefer et al. 2007) stations during a 5-yr period. In the top 0.1-m soil layer, the RMSE was 0.016 m3 m−3. The research sites were located on cropland, grassland, and shrubland. In contrast, Yang et al. (2011) found, during a 10-yr simulation run, that Noah-MP either significantly overestimated monthly SWC or excessively depleted the top 1-m soil profile in summertime, particularly in dry seasons in Illinois (Illinois Climate Network; Hollinger and Isard 1994). On a winter wheat stand in southern Germany, Ingwersen et al. (2011) showed that the Noah LSM tends to uniformly deplete the soil profile; this overestimates SWC in the top (0.15 m) soil layer and underestimates it in deeper layers.
The Noah-MP LSM was built on the original Noah LSM (Niu et al. 2011). It has improved physics, including dynamic vegetation, introduction of a simple groundwater model, and a three-layer snowpack. The Noah-MP has a choice of parameterization options in leaf dynamics, canopy stomatal resistance, a soil moisture factor for stomatal resistance, runoff, and groundwater. Another reason to test the Noah-MP LSM in our study is its coupling with the Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008), where Noah-MP LSM represents a land component.
Spatial variability of SWC is often evaluated by exploring the relationship between the spatial mean of SWC 〈θ〉 and the corresponding standard deviation σθ (e.g., Famiglietti et al. 2008). In the following, this relationship will be denoted σθ–〈θ〉. The range of all possible σθ–〈θ〉 states forms a phase space. Poltoradnev et al. (2016) suggested creating a σθ–〈θ〉 phase space from a polygon that is stretched upon the cloud of observed data and two anchor points: σθ at the permanent wilting point σθ–〈θwp〉 and σθ at saturation σθ–〈θs〉. The lower edge of the σθ–〈θ〉 envelope reflects the regional variability associated with the variability of soil water retention curves.
The objective of the present study was to evaluate the performance of Noah-MP in simulating the soil water dynamics in plowing horizons of two agricultural landscapes in southwest Germany. In particular, we focused on the question of how well the model reproduces the spatial variability of SWC. In this context, we followed the concept that a reasonably performing model should be able to produce σθ–〈θ〉 data pairs falling into the observed σθ–〈θ〉 envelope. For this purpose, we performed Noah-MP simulations and tested the model performance against a dataset of regional SWC observations (spatial extent: 27 km × 27 km) from our previously published work (Poltoradnev et al. 2016).
2. Materials and methods
a. Study sites
The study area has already been described in detail in Poltoradnev et al. (2015) and Wizemann et al. (2015). The Kraichgau (KR) region is covered with loess soils (predominantly luvisols; WRB 2006). The mean temperature is about 9°C. Annual precipitation ranges between 720 and 830 mm. The dominant soil type of the Swabian Alb (SA) region is leptosol (WRB 2006) with shallow stony solum and prevailing clay fraction. The mean temperature is 6°–7°C, and annual precipitation ranges between 800 and 1000 mm.
b. Soil moisture networks and eddy covariance stations
In total, 42 soil moisture stations were installed in KR and SA during the summer of 2009. Per region, 21 stations were distributed across the three spatial domains (inner, middle, and outer, with spatial extents of 3 km × 3 km, 9 km × 9 km, and 27 km × 27 km, respectively), as shown in Fig. 1. Soil moisture was measured on cropland with time domain transmission (TDT) soil moisture sensors (SI.99 Aquaflex Soil Moisture Sensor, Streat Moisture Solutions, United States) buried 0.15 m deep. Rainfall was recorded with a tipping-bucket rain gauge (resolution of 0.2 mm per tip) installed 1.8 m above ground.
Six eddy covariance (EC) stations were erected on arable land in 2009 (Ingwersen et al. 2011; Wizemann et al. 2015). In KR, three EC stations (hereafter EC1, EC2, and EC3) were set up near the city of Pforzheim close to the “Katharinentaler Hof” (48.9°N, 8.7°E), about 30 km south from the center of the KR soil moisture network. Another three EC stations (EC4, EC5, and EC6) were set up close to the village of Nellingen (48.5°N, 9.8°E) in SA, about 28 km from the SA soil moisture network center. More information on the soil moisture networks and EC measurements is given in Ingwersen et al. (2011), Poltoradnev et al. (2015), and Wizemann et al. (2015).
c. Tillage and crops within the study sites
Local soil properties determine different soil tillage strategies in KR and SA. Because of water erosion concerns, KR farmers grub instead of plow. Grubbing depth ranges between 0.15 and 0.20 m. In contrast, the common practice for SA soils is plowing. However, the typical depth is only about 0.10 m, because the soils are shallow. In KR, winter wheat is the most commonly cultivated crop, whereas in SA, the crop rotation is more diverse (Table 1).
d. Field measurements
The continuous SWC and rainfall observations were obtained from soil moisture networks over a 3-yr period from 2010 to 2012. Stations that contained gaps in the SWC series were excluded from the analysis. Thus, the KR data pool contained SWC measurements from 17 stations in 2010 and 2011 and 16 stations in 2012. The data pool of SA consisted of SWC records from 15 stations in 2010, 14 stations in 2011, and 10 stations in 2012. Further, seasonal datasets were extracted from the 3-yr data pools. These datasets covered the period from 1 April to 31 October, that is, 214 days, in 2010, 2011, and 2012. Then, the original 15-min resolution data were aggregated to 30-min, hourly, and daily averages. The spatial mean 〈θ〉 and standard deviation of SWC σθ were estimated using the standard formulas.
On the day of installation at each station, soil samples were collected along the sensor area from the same depth. Particle size distribution and SWC at the wilting point θwp (pF = 4.2) were determined in the laboratory. Fractions of sand, silt, and clay were measured with the standard pipette method (Scheffer and Schachtschabel 2008). Soil texture was classified using the U.S. Department of Agriculture (USDA) scheme (Soil Survey Division Staff 1993). To determine θwp, soil samples were placed on a porous ceramic plate at −1.5 MPa pressure. After equilibration, θwp was determined gravimetrically by drying at 105°C. SWC at saturation θs was estimated from bulk density ρs (g cm−3) as follows:
where the coefficient ρf stands for particle density (2.65 g cm−3). Soil bulk density was determined on five 100-cm3 soil cores taken close to the sensor from 0.15-m soil depth on the day of sensor installation.
Crop evapotranspiration (ETc; mm day−1) was calculated using the FAO Penman–Monteith approach (Allen et al. 1998). More details on calculation steps are given elsewhere (Poltoradnev et al. 2016). From the difference between total rainfall R and ETc, we computed the seasonal water balance (SWB).
The weather data (wind speed, wind direction, air temperature, relative humidity, atmospheric pressure, radiation flux) as well as soil temperature and SWC that were used for model initialization were measured in 30-min time lags at EC stations under winter wheat cultivation.
Within the footprint of every EC station, the leaf area index (LAI) was measured biweekly. For that purpose, at the beginning of each season five subplots of 4 m2 were randomly selected within the footprint. LAI was tracked at the central square meter of every subplot with an LAI-2000 Plant Canopy Analyzer (LI-COR Biosciences Inc., United States) throughout each season. The phenological development of 10 labeled plants (wheat, rape, or maize) per subplot was assessed using the Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie (BBCH) growth stage code (Meier 2001). The green vegetation fraction (GVF) was measured at KR EC stations in 2012 and 2013. Color photos were taken weekly at 1-m height above the canopy with a digital camera at five permanently marked 1-m2 plots within each field (Imukova et al. 2015). The GVF was not measured in SA, but was adapted by shifting corresponding KR GVF observation dates by 2 weeks. This is the typical time lag between vegetation periods and farmers’ activity between SA and KR.
e. The σθ–〈θ〉 phase space
To indicate possible σθ–〈θ〉 states, we used the PlotRegionHighlighter package (Noma 2013) available as an R software routine (version 3.1.1, Lucent Technologies Bell Laboratories Innovations, France). To compute the σθ–〈θ〉 envelope for each region, we combined observations of three seasons and two anchor points (σθ–〈θwp〉, σθ–〈θs〉) to form a dataset. Then, we generated a polygon, stretched over the outermost data points, for each dataset using the PlotRegionHighlighter routine.
As discussed in our previous study (Poltoradnev et al. 2016), the lower boundary of a σθ–〈θ〉 phase space reflects the spatial variability of soil water retention (pF) curves. To build an ensemble of retention curves, we computed SWC at given matric potentials Ψm, denoted as SWC(Ψm), starting at 0 hPa and ending at −16 000 hPa with a 10-hPa lag, for the topsoil of each station. For this, we used two hydraulic functions: van Genuchten–Mualem (VG; van Genuchten 1980) and Clapp and Hornberger (CH; Cosby et al. 1984). The VG parameters were estimated using Rosetta Lite version 1.1 (USDA; Schaap et al. 2001) based on measured soil texture and bulk density. CH parameters were derived from fitting the CH function to the VG retention curve. Initial parameters of CH were taken from Cosby et al. (1984). To judge CH versus VG curves, the index of disagreement F (Nash and Sutcliffe 1970) was computed as
wherein Pi and Oi denote SWC(Ψm), computed with CH and VG, respectively. The F value was minimized by changing the CH parameters using the Solver Excel add-in (version 2010, Microsoft Ltd.). The outputted CH retention curves are hence as close as possible to the VG curves. In this way, we computed 〈θ〉 and σθ as a function of Ψm for every soil moisture network (KR and SA).
To characterize the upper bound of the σθ–〈θ〉 envelope, we computed the greatest σθ–〈θ〉 value theoretically possible. To achieve this we assumed that one of the two utmost SWC states (θs or θwp) per station materializes at a given point in time. This results in 221 (2 097 152) combinations of binary SWC states per soil moisture network. To compute σθ and 〈θ〉 of all possible combinations, a Fortran computer code was written (Annex 1; see online supplemental material). After this, the entire dataset was grouped into 〈θ〉 classes at 1 Vol.% intervals (i.e., 15 ≤ 〈θ〉 < 16 Vol.%, 16 ≤ 〈θ〉 < 17 Vol.%, etc.; Vol.% is the volumetric percentage). The maximum σθ and the corresponding 〈θ〉 of each class were determined and assigned to an array giving the upper bound curve of σθ–〈θ〉 envelope.
f. The Noah-MP LSM (version 1.1)
1) Soil water regime
In the Noah-MP LSM (version 1.1), SWC dynamics are simulated based on the Richards equation:
Here, θ (m3 m−3) stands for the volumetric soil water content, t (s) is time, z (m) denotes the vertical coordinate (positive upward), D(θ) (m2 s−1) is the soil water diffusivity, K(θ) (m s−1) is the hydraulic conductivity, and F(θ) (s−1) is a sink term for root water uptake. In the present study, we tested two functions to describe K(θ) and D(θ). The first is the CH approach by Cosby et al. (1984) used in Noah-MP by default:
where Ks denotes the saturated hydraulic conductivity; Ψm and Ψs stand for the matric potential and matric potential at the air entry point, respectively; and b determines the slope of the retention curve. The second function, the VG approach according to van Genuchten (1980), reads
where Se is the effective water content and θr represents the residual water content, while α and m = (1 − n)/n are parameters.
2) Model parameterization
Simulations were run with a time step of 1800 s (30 min) point-by-point for each station of the two soil moisture networks left in the analysis (see section 2d) in order to get the same dataset of SWC as observed. The spinup time was set to 3 months prior to April of each year, that is, the first day of every simulation run was 1 January. The soil profile was divided into five layers (0.1, 0.1, 0.3, 0.5, and 1.0 m thick). Accordingly, the total thickness of the soil profile was set to 2 m. The initial temperatures and water contents of the soil layers were derived from EC field measurements. The soil temperature at 2-m depth was set to the mean temperature. A free drainage boundary condition was set at the bottom of the soil column. The vegetation parameters were based on the U.S. Geological Survey (USGS) “dryland cropland and pasture” land cover class. The main rooted zone was assumed to extend to the uppermost 0.3 m.
3) Simulation setup and output
To determine the contribution of different factors to the spatial variability of SWC, rainfall, soil texture, and vegetation were stepwise included as spatial variables. The modification process included a stepwise introduction of station-by-station observed rainfall, soil texture, LAI, and GVF data into the model input files and consequent replacement of the default data layout. These model runs were performed only for 2010. The 2011 and 2012 observations were investigated only with the full set of influencing factors. Each simulation produced a set of SWC data for the second soil layer (0.15-m soil depth). From these, we estimated 〈θ〉 and σθ and compared them with the corresponding statistical parameters of the observed data.
(i) S1: Simulation with spatially variable rainfall
In the first set of simulations (S1), rainfall was introduced as the only factor varying among stations. Rainfall data were derived from the observations at the stations and aggregated to 30-min time spans. Soil texture was set to be uniform among all stations, defined as silty loam for KR and silty clay for SA. The soil parameters used in the CH equation were taken from Cosby et al. (1984). Class average VG parameters were estimated with Rosetta Lite version 1.1 (Schaap et al. 2001). The default monthly mean LAI and GVF were taken from the USGS dryland cropland and pasture category. Monthly LAI and GVF values from January to December were set to 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 1.5, 0.0, 0.0, 0.0 and 0.05, 0.05, 0.25, 0.32, 0.93, 0.96, 0.25, 0.05, 0.05, 0.05, 0.05, 0.05, respectively.
(ii) S2: Simulation with spatially variable rainfall and soil texture
The soil texture of measuring locations varied among the USDA texture classes as follows: in KR, silt loam (N = 17) and silty clay loam (N = 4); in SA clay (N = 8), silty clay (N = 8), silty clay loam (N = 4), and silt loam (N = 1) (Fig. 2). The S2 simulation was divided into two parts. In a first variation (S2a), hydraulic parameters specific for the respective soil texture class were used. Parameters of the CH function are provided in the Noah-MP lookup table by default, while VG parameters for the modified Noah-MP were taken from Table 3 of Cosby et al. (1984). In a second variation (S2b), VG parameters were estimated with Rosetta Lite version 1.1 based on soil texture and bulk density, measured at each station. Site-specific CH functions were fitted to the VG retention curve as described in section 2e. Monthly mean LAI and GVF were kept at the default values.
(iii) S3: Simulation with spatially variable rainfall, soil texture, and vegetation dynamics
In the S3 simulations, we additionally considered the heterogeneity of vegetation. Similar to Imukova et al. (2015), the croplands were divided into four groups: “early covering crops winter wheat like” (ECC1; winter wheat, summer barley, winter barley, triticale, spelt, pea, and oats), “early covering crops winter rape like” (ECC2; here only winter rape), “late covering crops silage maize like” (LCC; silage maize, corn, sugar beet, and sunflower) and “grassland.” LAI and GVF data were compiled for crops and regions separately. LAI and GVF data for the ECC1, ECC2, and LCC groups were derived from winter wheat, winter rape, and silage maize data measured at the EC stations (Wizemann et al. 2015). The default LAI and GVF were used for “grassland.” Values for the twenty-fourth day of every month of a year were derived by linear interpolation. LAI and GVF dynamics of the different crop groups are shown in Fig. 3.
(iv) S4: Monte Carlo–based sensitivity analysis
To determine the contribution of the spatial distribution of soil hydraulic properties and atmospheric forcing into the σθ–〈θ〉 relationship, two additional simulations were run for KR. The S4a set of simulations was performed in the same way as S2b (VG approach). Yet, instead of deriving Ks from Rosetta, Ks values were drawn from a lognormal distribution with the mean and standard deviation taken from Wang et al. (2013). These authors measured Ks in the surface layer of a loess soil. The S4b simulation was also set up like S2b (VG approach), but with variable atmospheric forcing. To model the variability of weather conditions, weather data were collected from the EC stations and five German Weather Service (DWD; Offenbach am Main, Germany) weather stations, located close to the soil moisture networks. Wind speed, wind direction, air temperature, relative humidity, atmospheric pressure, and radiation fluxes were randomly generated for each soil moisture station within the range of the mean standard deviation of each weather variable.
4) Model evaluation
where Pi stands for the predicted value, Oi is the corresponding observed value, and n is the number of observations.
Bias and model efficiency were computed as
where denotes the mean of Oi.
To verify the performance of a simulation in regard to the spatial variability of SWC, we computed the fraction of data points within the σθ–〈θ〉 envelope, the mean distance from simulated data points to the lower edge of the envelope, and the distance between measured and simulated cluster centers. The R “sp” package (Pebesma and Bivand 2005; Roger et al. 2013) was used to calculate the number of data points falling into the σθ–〈θ〉 envelope.
a. Comparison of Noah-MP model runs with observations of the 2010 season
1) 〈θ〉 dynamics
Figure 4 shows the observed and simulated 〈θ〉 dynamics in KR and in SA during the 2010 season. For both networks, S3 performed best (Fig. 4, Table 2). In general, 〈θ〉 dynamics was reproduced better for the loess soils of the KR region than for the shallow clayey soils of SA. In KR, when only spatially variable rainfall (S1) or both the rainfall and texture class (S2a) were considered, 〈θ〉 was overestimated during April, early May, and the June–July time periods using the CH approach. Modeled 〈θ〉 immediately responded to rainfall, but depleted rather slowly. In contrast, simulated 〈θ〉 drying processes tended to be faster than observed during August. VG performed better at the beginning of the vegetation season but delivered systematically lower 〈θ〉 from the end of June onward. In SA, for S1 and S2a the CH approach underestimated 〈θ〉 in spring and fall and significantly overestimated 〈θ〉 during the summer drought. VG strongly underestimated the observed 〈θ〉 during the entire observation period. In both regions, adding a site-specific parameterization of the hydraulic functions (S2b) significantly improved the model performance. Nevertheless, the topsoils of SA were always wetter than observed, except for certain time periods in spring and fall. Parameter 〈θ〉 was overestimated during April drying periods and during summer drought (June–July) in KR. The underestimation of soil depletion during drying phases was greater when using CH than with VG. A further improvement was achieved by introducing heterogeneous vegetation (S3). This brought both simulated 〈θ〉 dynamics closer to the measurements: EF and R2 significantly increased, and RMSE and bias decreased (Table 2). In SA, however, Noah-MP continued to deliver systematically higher 〈θ〉 than observed, except for one month in spring. The highest R2 and EF along with the lowest RMSE and bias were delivered by VG S3 in KR.
2) Spatial SWC variability
Simulated and observed σθ–〈θ〉 data pairs are plotted in Fig. 5 against the constructed envelopes. The model performance statistics are given in Table 2. In both regions, S1 resulted in the poorest Noah-MP simulation. In KR, a great improvement was achieved by introducing spatially variable soil texture classes. The VG parameterization performed better than the CH approach. Moreover, VG S2a simulations delivered the highest fraction of modeled data points falling into the envelope, among all runs. It brought σθ–〈θ〉 on average slightly above the lower edge and greatly reduced the distance between two cluster centers (Table 2, Fig. 5; S2a). For KR, simulated σθ converged to the left anchor point during the long drying-out phase and showed the same pattern as the measurements. Parameter σθ increased as 〈θ〉 decreased in the beginning of the drying period, reached the minimum point, and increased with further 〈θ〉 decrease thereafter. Introducing station-specific soil texture (S2b) and spatially heterogeneous vegetation dynamics (S3) brought less data into the envelope than S2a, though the simulated and observed data clusters approached each other.
In SA, the simulated data did not fall into the σθ–〈θ〉 envelope when the spatially heterogeneous soil texture was considered (Fig. 5, S2a). Considering station-specific soil texture (S2b) and, further, spatially variable dynamic vegetation (S3) positively influenced model performance. The maximum fraction of data points in the envelope was achieved with the CH S3 setup. S2b and S3 simulations formed a cloud of σθ–〈θ〉 data.
Both approaches could not adequately represent the SWC variability at the lower bound of the σθ–〈θ〉 envelope (black solid and dashed lines, Fig. 5; S2a). When a station-specific soil texture was introduced (S2b), the lines that reflect the variability within the water retention curve ensemble were located closer to the lower boundary of the σθ–〈θ〉 envelope.
The upper boundary of the σθ–〈θ〉 envelope had a top-closed convex shape, and the observed σθ values were positioned inside the polygon. Theoretically, however, measured values could be greater. The dashed curve in Fig. 6 indicates the maximum possible variability of SWC as estimated in the binary soil water scenario.
The results of the sensitivity analyses are shown in Fig. 7 and are summarized in Table 2 (VG S4a,b). Each of VG S4 simulations brought more data points into the envelope than VG S2b (KR, Table 2). Among all runs, VG S4b yielded the highest positive mean distance between the data points and the lower envelope boundary and also the lowest distance between cluster centers of modeled and observed data.
The temporal dynamics of observed σθ versus σθ modeled with different Noah-MP setups are presented in Fig. 8. The statistical moments (temporal mean and σ) are given in Table 3. Based on the σθ–〈θ〉 relationships, the model generally failed to reproduce the observations. However, in KR the S2a and S3 variations yielded σθ values, the temporal variation of which was close to the σ of observed σθ (Table 3).
b. S3 simulations from 2010 to 2012
Figure 9 shows observed and simulated (S3 simulations) SWC dynamics in the 2010, 2011, and 2012 seasons. As expected, the response of the simulated 〈θ〉 to rainfall was well pronounced in all three years. This led to frequent overestimation of wetting and 〈θ〉 magnitude in both regions. In KR, VG tended to underestimate SWC in early April 2010 and April–May 2011 and 2012. Both CH and VG approaches repeatedly depleted the topsoil during the second half of July and August 2011 and 2012. The statistics show that the highest R2 and EF values (Table 4) were observed in the season with positive SWB (2010, Table 5), while the lowest was found in the season with negative SWB (2012, Table 5). In SA, the best model fit was achieved in 2011, when the lowest SWB among the seasons occurred. In 2012, both approaches underestimated 〈θ〉, particularly during the second half of the season.
The σθ–〈θ〉 pattern of the simulated SWC dynamics poorly reflected the observed SWC variability in all three years. This holds true for both regions (Fig. 10). In KR, VG performed better than CH in all three seasons (Table 6). The year 2011 was the year in which the largest fraction of modeled data points was within the envelope. In SA, the performance of VG and CH was different among the years. In 2010 and 2011, CH showed the best match and located more data points into the σθ–〈θ〉 envelope. In 2012, VG yielded a better match.
In the present study, the best match between the mean spatial average of simulated and measured SWC values was achieved when the model used site-specific rainfall, soil texture, GVF, and LAI. Additionally, the choice of the hydraulic functions affected the performance of Noah-MP. The latter was best with VG. Gayler et al. (2013) suggested that one of the explanations for the better performance of VG parameterization in their work is that the VG model is more flexible than CH. It represents the inflection point of the retention curve better (data not shown). Moreover, with VG the unsaturated hydraulic function sloping was smoother than with CH (data not shown), which may have had a strong effect on SWC modeling. The same parameterization, however, performed less accurately for the seasons with neutral and negative SWB. During these seasons the observed 〈θ〉 at 0.15-m depth often did not respond to a rainfall event, while the response of modeled 〈θ〉 was well pronounced in all three years. One possible reason for this mismatch is that the modeled infiltration rate was too high, resulting in a sharp increase of SWC at 0.15 m. Noah-MP uses a simple water balance scheme (Schaake et al. 1996) to partition incoming rainfall into infiltration and runoff at the soil surface. Chen and Dudhia (2001) have pointed out that some of the parameters used in the simple water balance scheme are purely empirical. They depend on the precipitation climatology of a study region and have to be calibrated. Another reason for the mismatch may be a bias in simulated evapotranspiration, which controls the pace of dissipation of SWC. Ingwersen et al. (2011) compared latent heat flux, simulated by the Noah LSM, against the eddy covariance data measured in 2009 at a winter wheat stand in KR. The model underestimated latent heat flux in April–June and considerably overestimated it in July and August. The researchers could markedly improve the simulations by introducing a monthly varying minimum stomatal resistance into the Noah LSM.
The reasons for the better match of simulated and observed 〈θ〉 in KR are twofold. First, the SA soils are usually shallow. The plowing horizon often directly overlies weathered lime rock. In Noah-MP, total soil thickness is set to 2 m by default. This unavoidably causes systematic errors, in particular with regard to the lower boundary condition. We applied a free drainage condition at the bottom of the soil column, precluding a possible influence of capillary rise on SWC. In the study areas, we can exclude capillary rise for the following reasons. SA is a karst region, so the groundwater is deeper than 25 m below the surface. In KR, the elevations of the soil moisture stations are about 10–15 m above the river basin level. In presimulations, we found that setting the soil column thickness to the observed one depleted the soil water profile in the simulations (data not shown). Therefore, we decided not to change the default setup. Second, SA soils are rich in clay and stones. VG is known to be less accurate in clay soils (van Genuchten 1980), and CH has known deficiencies in representing the hydraulic properties of rocky soils (Clapp and Hornberger 1978).
In both regions, Noah-MP poorly represented the spatial variability of SWC. In most cases, the greatest share of simulated σθ–〈θ〉 data points was located below the bottom edge of the envelope. The lower bound of the σθ–〈θ〉 phase space reflects the case of a spatially homogeneous matric potential. A position of σθ below the envelope indicates that Noah-MP smooths the spatial variability of SWC. Introducing site-specific parameters improved model performance, but the spatial variability of SWC was still systematically underestimated. The following factors contribute to the spatial variability of SWC. First, the model does not consider topography. It is well known that elevation, slope, and hillslope position (nonlocal controls) contribute considerably to SWC spatial variability (Grayson et al. 1997; Albertson and Montaldo 2003). Second, because of the restriction of a uniform root density profile, Noah-MP tends to uniformly deplete the soil column (Fan et al. 2006; Ingwersen et al. 2011). This results in weaker vertical SWC gradients compared to observations (Ingwersen et al. 2011). Gayler et al. (2014) have shown that including dynamic root growth may markedly improve the performance of Noah-MP. Third, Ks values observed in the field are lognormally distributed, while those calculated with pedotransfer functions are normally distributed (Wang et al. 2013). Using lognormally distributed Ks values increased the spatial variability of SWC. Fourth, the standard Richards equation as implemented in Noah-MP does not consider water movement along macropores or cracks. Wöhling et al. (2009) concluded in their HYDRUS-3D (Šimůnek et al. 2006) simulation study that model structural inadequacy is directly linked to the uniformity of modeled water flow, which is certainly not realistic (Wiekenkamp et al. 2016). Koch et al. (2016) also emphasized the systematic underestimation of spatial variability of SWC at the Wüstebach catchment site in Germany by three hydrologic distributed models that solve the Richards equation in three dimensions. They explained this failure mainly with the poor performance of the VG function during drought periods. This seems in contrast to our findings. Referring to Fig. 8, the temporal dynamics of σθ simulated with VG approaches the observed dynamics during the long dry period in July 2010 (when variability in SWC is mainly attributed to the variability in soil texture and residual water content), while Noah-MP tended to gradually decrease SWC variability during spring and fall. Fifth, spatially constant atmospheric forcing also leads to the homogenization of spatial SWC by the model (Albertson and Montaldo 2003). Monte Carlo simulations with variable weather conditions made more simulated σθ–〈θ〉 points fall into the envelopes. In the review paper by Crow et al. (2012), the meteorological forcing was also mentioned as one of the factors dominating large-scale soil moisture variability (observed at regional and larger scales). They stressed that the mean σθ for intermediate soil moisture conditions increases with increasing spatial extent. Sensor-to-sensor variability as well as sensor noise may also increase the overall variability of observed SWC. However, based on previous studies, sensor noise is usually much lower than sensor-to-sensor variability, and the latter is substantially decreased by sensor-specific calibration (Rosenbaum et al. 2010; Qu et al. 2013).
Huang and Margulis (2013) studied the impact of soil moisture variability on atmospheric boundary layer (ABL) characteristics, based on a coupled LSM–LES model. Introducing the SWC variability into the land–atmosphere system resulted in the formation of mesoscale circulations. Clouds and cloud-covered areas increased as soil moisture variability increased. The SWC variability length scale (extension of dry and wet regions) amplified this effect. The simulation design contained two configurations of dry–wet pixels, with initial SWC values set to 0.33–0.43 cm3 cm−3 and 0.28–0.48 cm3 cm−3. These dry–wet patterns were organized via chessboard-like surfaces with different length scales (1.25, 2.5, and 5 km) within a 10 km × 10 km domain. Using the input data, one can estimate that the minimum σθ of the domain is 5.0 Vol.%, corresponding to a variability length scale of 1.25 km. The maximum σθ was 11.6 Vol.% at the 5-km length scale (Huang and Margulis 2013). In our study the subgrid σθ of the middle domain (9 km × 9 km, including five stations of the inner domain and eight stations of the middle domain) of KR and SA networks were on average 4.4 and 7.0 Vol.%, respectively, in 2010.
The measuring volume of the TDT sensor used in our investigations is a 3-m-long cylinder with a diameter of several centimeters. The measurement is hence representative for a few square meters. Noah-MP, coupled with WRF, can be used down to the LES scale (Bauer et al. 2016). In most weather and climate simulations, however, Noah-MP has been used with a coarser horizontal grid spacing (12 km or more). On such a grid, σθ–〈θ〉 will take smaller values than at the square-meter scale. The performance of such simulations should be tested against measurements with larger support areas, employing remote sensing or hydrogeophysical techniques (e.g., Binley et al. 2015). Recently, high-resolution soil moisture products have gained significance and demand. Importantly, there is a branch of literature devoted to the development of downscaling methods that increase the resolution of soil moisture within a coarse satellite pixel (Mohanty et al. 2017). In hydrology, for example, Mascaro and Vivoni (2012) demonstrated that disaggregation of remotely sensed soil moisture data (from 25.6 to 0.8 km) improved the performance of a hydrologic model within the Sierra Los Locos catchment in Sonora, Mexico. The downscaled high-resolution product delivered a more realistic spatial SWC distribution than the coarse product. This, in turn, positively influences the accuracy of the model prediction of surface fluxes, particularly evapotranspiration and runoff (Vivoni et al. 2010). Underestimation of (or leveling out) the spatial variability of soil moisture modeled by the Noah-MP needs to be further investigated, if one intends to use Noah-MP on finer scales.
5. Summary and conclusions
We evaluated the performance of Noah-MP against distributed time series from a topsoil water content measurement network in two regions in southwest Germany. The water content was measured by TDT and is representative for an area of a few square meters. With regard to the regional mean soil water content 〈θ〉, Noah-MP performed fairly well for the loess soils of KR, but it has limitations as to the shallow, clayey, and stony soils of SA. The choice of hydraulic functions impacted the model output substantially. The best match was achieved with the VG approach and with site-specific rainfall, soil texture, GVF, and LAI input data. Simulation of 〈θ〉 is the best in seasons with a positive water balance (R > ETc). However, the spatial variability of SWC is systematically underestimated by Noah-MP. The model tends to gradually eliminate the variability, particularly during the wet periods of a year (spring and autumn), but during long dry-down periods the observed σθ dynamics is represented considerably better (particularly the VG configuration). The majority of the modeled σθ–〈θ〉 pairs were located outside the envelopes of the measured data and below their lower bounds. The insufficient representation of the spatial SWC variability by Noah-MP can be mainly attributed to missing topography and terrain information, inadequate representation of the spatial variability of soil texture and hydraulic parameters, and the restriction to a uniform root distribution.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-17-0169.s1.