Simulated soil moisture is increasingly used to characterize agricultural droughts but its parametric uncertainty, which essentially affects all hydrological fluxes and state variables, is rarely considered for identifying major drought events. In this study, a high-resolution, 200-member ensemble of land surface hydrology simulations obtained with the mesoscale Hydrologic Model is used to investigate the effects of the parametric uncertainty on drought statistics such as duration, extension, and severity. Simulated daily soil moisture fields over Germany at the spatial resolution of 4 × 4 km2 from 1950 to 2010 are used to derive a hydrologically consistent soil moisture index (SMI) representing the monthly soil water quantile at every grid cell. This index allows a quantification of major drought events in Germany. Results of this study indicated that the large parametric uncertainty inherent to the model did not allow discriminating major drought events without a significant classification error. The parametric uncertainty of simulated soil moisture exhibited a strong spatiotemporal variability, which significantly affects all derived drought statistics. Drought statistics of events occurring in summer with at most 6 months duration were found to be more uncertain than those occurring in winter. Based on the ensemble drought statistics, the event from 1971 to 1974 appeared to have a 67% probability of being the longest and most severe drought event since 1950. Results of this study emphasize the importance of accounting for the parametric uncertainty for identifying benchmark drought events as well as the fact that using a single model simulation would very likely lead to inconclusive results.
Drought is a recurrent and extensive climatic phenomenon characterized by below-average water availability whose duration might last for several years. It is considered to be one of the most costly natural disasters because it often induces huge socioeconomic losses (Wilhite 2000) as well as environmental degradation. During the summer of 2003, for instance, several parts of Europe endured the highest temperatures of the last 500 years (Luterbacher et al. 2004; Fink et al. 2004) and one of the most extensive and severe drought in records. In Germany alone, the estimated loss in the agricultural sector was 1.5 billion Euros (COPA-COGECA 2003). In extreme cases, prolonged drought spells might lead to unprecedented environmental disasters often associated with the decline of human societies (Hodell et al. 1995; Haug et al. 2003) or the trigger for mass migrations and famine (Field 2000). Droughts occur indifferently in high and low rainfall areas and in virtually all climatic zones (Dracup 1991; Mishra and Singh 2010), although the most severe human consequences happen in arid regions.
Currently, hydrometeorological mechanisms originating droughts are relatively well understood. In general, droughts are driven by extreme macroclimatic variability originated by atmospheric interactions and feedback between the atmosphere, the oceans, and the land surface (e.g., Nicholson 2000; McCabe and Palecki 2006). This variability is, in turn, related to the solar activity as well as atmospheric composition, and strongly affected by anthropogenic activities (Sheffield et al. 2009).
Our ability to make reliable drought predictions, however, is not satisfactory (Wilhite 2000), although there is vast scientific literature on this topic. One of the main reasons is related to the insufficient knowledge regarding the processes controlling drought development and persistence, as well as its spatiotemporal variability (Sheffield et al. 2009). Another reason stems from the fact that there is no clear definition of this phenomenon (Wilhite and Glantz 1985) since it depends upon the variable that is used for its quantification.
Droughts have been mainly classified into three types: 1) meteorological drought, usually defined as an extreme anomaly of precipitation; 2) hydrological drought, which is related to a deficit in the supply of surface and subsurface water; and 3) agricultural drought, being a combination of meteorological and hydrological droughts leading to deficits in root zone soil moisture available to vegetation (Wilhite and Glantz 1985). Since precipitation and discharge data are widely available, a plethora of drought indices have been proposed in the scientific literature to quantify meteorological and hydrological droughts—for instance, the Palmer drought severity index (Palmer 1965), the standard precipitation index (McKee et al. 1993), and the regional deficiency index (Stahl and Demuth 1999), among others.
It is widely accepted, however, that these empirical indices are not adequate to represent extreme water stress conditions that would lead to a significant reduction of biomass and crop yield (Keyantash and Dracup 2002; Mishra and Singh 2010). In Germany, for example, Döring et al. (2011) have shown that empirical drought indices based only on available data such as precipitation and temperature do not constitute adequate measures to describe agricultural drought stress because they do not explicitly account for the available water stored in the root zone, which is ultimately the plant’s life-supporting substance.
Direct soil moisture observations, on the other hand, are not available at the regional level because measuring this variable at large scales is not logistically and economically feasible (Vereecken et al. 2008). This implies that hydrologic or land surface models would have to be employed for the estimation of the soil water content. Soil moisture, in contrast to precipitation or discharge, constitutes a good index for quantifying agricultural drought because it controls the proportion of the rainfall that percolates, runs off, or evaporates from the earth surface (i.e., root zone). Concisely, it integrates precipitation and evapotranspiration as well as the delays introduced by interception, snow accumulation, and melting over periods of days to weeks. In other words, soil moisture in the root zone is a governing factor sustaining vegetative growth and thus it is a direct indicator of agricultural drought (Keyantash and Dracup 2002). Land surface models such as the three-layer Variable Infiltration Capacity model (VIC-3L; Liang et al. 1996) and Safran–Isba–Modcou (SIM; Soubeyroux et al. 2008), for example, have been used recently to assess agricultural drought characteristics in the United States and France, respectively (Sheffield et al. 2004; Andreadis et al. 2005; Vidal et al. 2010). There are, however, several key issues that should be considered if simulated soil moisture is chosen for quantifying agricultural droughts.
Modeling soil moisture dynamics at large scales (e.g., grid cells greater than 500 m) is difficult and uncertain as was demonstrated by the Project for the Intercomparison of Land-Surface Parameterization Schemes (PILPS) project (Chen et al. 1997). In this project, 23 land surface models (LSMs) exhibited significant differences between modeled and measured soil moisture (among other variables) although all models were based on fundamental principles of mass and energy conservation and forced with identical atmospheric conditions. This experiment also indicated the existing interplay between this state variable and other fluxes such as latent heat as well as the substantial parameter uncertainty that is related with these physical processes. At larger scales, the subgrid variability of the variables involved and the nonlinearity of the processes make the modeling of soil moisture even more complicated because parameterization schemes might become scale dependent (Nykanen and Foufoula-Georgiou 2001). It should be noted that effective model parameters (e.g., saturated soil water content or porosity) at large scales can only be estimated but not measured. This, in turn, constitutes a new source of uncertainty that should be taken into account when modeling soil moisture dynamics. Consequently, a drought monitoring and early warning system based on a soil moisture index, which does not fully take into account the predictive uncertainty of the simulation model, might be inadequate for real applications and/or for impact assessment.
Most of the soil moisture drought studies (Andreadis and Lettenmaier 2006; Vidal et al. 2010; Shukla et al. 2011) found in the literature have not addressed the epistemic uncertainty related to parameterization, model structure, and input data. More recently, Wang et al. (2011) argued that state variables, such as soil moisture, are strongly dependent on the parameterization of the LSMs and the quality of the meteorological forcing data. Similar results have been found by Mo et al. (2012), who concluded that the primary source of uncertainty between two drought monitoring systems operated in the United States is originated from precipitation data and, in a minor degree, from air temperature, shortwave and longwave radiation, and wind speed. As a result, substantial discrepancies with in situ measurements have been found (Entin et al. 2000), which are mainly attributed to the variability of topography, soil, vegetation, and root structure, but could also stem from uncertainty sources mentioned above. Specifically, finding a robust parameterization scheme for a LSM or a hydrological model, which is able to produce reliable estimates of water fluxes at high spatial resolution over large domains, is one of the grand challenges of contemporary hydrology (Beven and Cloke 2012).
It has been noted, however, that multimodel ensembles are able to describe the anomalies and seasonal variability of soil moisture. Wang et al. (2009, 2011) successfully applied this technique to reproduce agricultural drought characteristics in the continental United States and China. In both studies, six LSMs were used to generate soil moisture fields for a period of almost 100 years in the United States and 56 years in China. However, in those studies, only a single simulation for each LSM was used.
In this study, we argue that a unique parameter set for a given LSM is inadequate to estimate water fluxes and related state variables at high spatiotemporal resolutions, considering that both inputs and model parameters over large modeling domains are subject to considerable uncertainties because of the reasons mentioned above (see also Rosero et al. 2011). Thus, we hypothesize that any drought characteristic (e.g., severity and duration) based on simulated soil moisture is prone to large variability due to parametric uncertainty, which, if it is not taken into account, will lead to incorrect estimates of drought characteristics.
The main objectives of this study and the rationale behind them are summarized below: 1) To obtain a consistent ensemble of daily soil moisture fields for Germany since 1950 at a spatial resolution of 4 × 4 km2. Such reconstruction is fundamental to characterize historical drought events and their related characteristics. To the best of our knowledge, this is the first study to perform nationwide agricultural drought reconstruction for Germany. Long-term soil moisture simulations are also fundamental for initializing hydrologic or regional climate models and the basis to fulfill the remaining objectives. 2) To develop a reliable soil moisture drought index (SMI) for Germany at a high spatial resolution. Such SMI is key for implementing a monitoring system and adaptation strategies at regional scale. Available global soil moisture analyses have a spatial resolution 0.5° or larger, which is too coarse for a regional drought analysis. 3) To identify benchmark agricultural drought events occurring in summer and winter in Germany during the last 60 years and the uncertainty of their main statistical characteristics. These exceptional events are necessary to identify potential climate change effects on the hydrological cycle. The uncertainty associated with drought characteristics—such as coverage area, duration, and severity—will be quantified by means of a Monte Carlo method. Ensemble model simulations would allow us to assess the reliability of the predictions which, in turn, will lead to minimize the number of false positive drought events (i.e., cases in which the SMI indicates that a given event is below a certain threshold for a given characteristic when in fact it is not). Additionally, the effect of the ensemble size on the false positive rate will be investigated. 4) To identify regions in Germany prone to strong drought persistence as well as areas exhibiting significant trends in monthly soil moisture fields. These insights would provide hints for mitigation and adaptation measures at regional scale.
2. Soil moisture data
Soil water availability in the root zone is a direct indicator of agricultural drought because it constitutes a governing factor of the state of vegetative growth through the availability of water for transpiration (Keyantash and Dracup 2002). Measuring soil moisture content over the entire domain of Germany at a spatial resolution of 4 × 4 km2, for example, is logistically and economically infeasible (Vereecken et al. 2008). LSMs or hydrologic models are therefore often employed to estimate this key variable over large spatial domains and longer periods (Andreadis and Lettenmaier 2006; Sheffield and Wood 2007; Wang et al. 2009; Mishra et al. 2010; Wang et al. 2011).
In this study, the mesoscale Hydrologic Model (mHM; Samaniego et al. 2010) was used to generate a large ensemble of daily soil moisture fields for the period 1950–2010. A three-layer soil scheme was used to model the soil moisture dynamics over the entire root zone depth (i.e., approximately up to 2 m below ground). The depth of the first two layers was fixed to 5 and 25 cm, whereas the depth of the last one was variable according to soil characteristics provided by the soil texture map. The spatial resolution of each grid was 4 × 4 km2 (level 1). A short description of mHM and the generation of ensemble soil moisture fields are given below.
a. The mHM
The mesoscale Hydrologic Model is a process-based water balance model (Samaniego et al. 2010) that has been developed over the last 5 years at the Helmholtz Centre for Environmental Research–UFZ. This spatially explicit model does not differ significantly from existing large-scale hydrologic models [e.g., the Hydrologiska Byråns Vattenbalansavdelning (HBV) and the VIC-3L model] on how dominant hydrologic processes at the meso- and macroscales are conceptualized, but on how the effective parameters of the model are quantified at a selected modeling scale and on how the subgrid variability of physiographic characteristics provided at level 0 is taken into account for the estimation of these effective parameters. These two fundamental differences constitute the core of the multiscale parameter regionalization technique (Samaniego et al. 2010) that is embedded into mHM. Extensive numerical experiments have shown that this technique is capable of coping with the large spatiotemporal variability of the input data and, as a result, mHM is able to produce quite good performance at multiple spatial resolutions and locations other than those used during calibration (i.e., proxy basin and flux-matching tests).
Currently, mHM has been evaluated in more than 100 basins in Germany ranging from 4 to 47 000 km2 (Samaniego et al. 2010; Kumar et al. 2010, 2012, unpublished manuscript). This model is driven by disaggregated fields of daily forcings such as precipitation, temperature, and potential evapotranspiration. It accounts for the following hydrological processes: canopy interception, snow accumulation and melting, evapotranspiration, infiltration, soil moisture dynamics in three layers, surface runoff, subsurface storage, discharge generation, percolation, baseflow, and flood routing within the river reaches. Readers may refer to Samaniego et al. (2010) for a detailed model description as well as its parameterization.
The morphological and physiographic data required for setting up mHM include a digital elevation model (50 × 50 m2) acquired from the Federal Agency for Cartography and Geodesy, a vector soil map containing information on soil textural properties such as sand and clay contents of different soil horizons, and a vector map of hydrogeological formations containing properties such as saturated hydraulic conductivity. Both vector maps at a scale of 1:1 000 000 were obtained from the Federal Institute for Geosciences and Natural Resources of Germany. Three Corine land cover seamless vector data (http://www.eea.europa.eu) for the years 1990, 2000, and 2006 were employed to account for the changes in states of land cover over the simulation time period (1950–2010). Land cover states, prior to the year 1990, were inferred from the Corine 1990 map. Monthly variability of the leaf area index was estimated for each land cover class with Moderate Resolution Imaging Spectroradiometer (MODIS) scenes from 2001 to 2009. These data are freely available from https://lpdaac.usgs.gov/get_data. For a detailed description on data processing and setting up mHM in several river basins, interested readers may refer to Samaniego et al. (2010) and Kumar et al. (2010). Previous datasets were resampled on a common spatial resolution of 100 × 100 m2 denoted as level 0. This level of information provides the subgrid variability of all morphological and physiographic variables required to run the model at any coarser resolution denoted as level 1 (e.g., 4 km). The time series of discharge data across several gauging stations were acquired from the European Flow Regimes from International Experimental and Network Data (Euro-FRIEND) program (http://ne-friend.bafg.de) and the Global Runoff Data Centre (http://www.bafg.de).
Gridded fields of daily average precipitation as well as maximum, minimum, and average air temperatures at 4 × 4 km2 spatial resolution (level 2) were estimated from their respective point measurement data from about 5600 rain gauges and 1120 meteorological stations, operated by the German Meteorological Service (DWD). Two interpolation techniques were used to derive the daily fields of precipitation, which are detailed in section 2b. Gridded estimates for temperature fields were obtained with external drift kriging, wherein the terrain elevation was used as a drift variable. The daily fields of potential evapotranspiration were estimated with the Hargreaves and Samani method (Hargreaves and Samani 1985) and were subsequently corrected to account for the spatial variability of the terrain aspect.
b. Ensemble description and experimental design
Two major sources of parametric uncertainty were identified through sensitivity analysis. The most important one is related with the variability of the global calibration parameters of mHM (i.e., space and time independent), and the second one is related with the parameters required for the rainfall interpolation method. Consequently, the uncertainty tree was divided into two main branches, each one driven by two independent interpolation methods but both based on the same rainfall measurements. These two branches were denoted as DWD1 and DWD2. Other meteorological variables such as daily, minimum, and maximum temperature required in both branches were kept the same. This assumption was taken considering 1) that precipitation interpolation is one the most important source of error in the input data (Mo et al. 2012) and 2) that the areal coverage of snow-dominated areas in Germany is geographically limited.
The DWD1 branch was created with external drift kriging using terrain elevation as a drift and a combined variogram that comprised a nugget and an exponential part. The resolution of this product was 4 × 4 km2, with daily time steps from 1950 to 2010. The best-fit parameters (i.e., nugget, range, and sill) were found through a cross-validation procedure.
The DWD2 branch was obtained by resampling the original daily Regionalisierung von Niederschlagshöhen (REGNIE; www.dwd.de) product available at 1 × 1 km2 into a regular grid similar to that of the DWD1 dataset. The k-nearest neighbor technique and a standard geo-referencing algorithm were employed for this purpose. The DWD1 data were used to complete this set with daily fields from 1950 to 1959 since the REGNIE dataset is only available from 1960 to 2010. The REGNIE data are based on multiple linear regression having elevation, geographic location, and aspect as predictors.
Within each branch, the propagation of the parameter uncertainty into the soil moisture simulations was evaluated by an ensemble of 100 best parameter sets of mHM. The following procedure was implemented for their selection. First, in every major river basin depicted in Fig. 1, the dynamically dimensioned search algorithm (Tolson and Shoemaker 2008) was employed to find good sets of global parameters which exhibit an acceptable model efficiency (e.g., Nash–Sutclife efficiency of at least 0.75) during the evaluation period (for details refer to Kumar et al. 2010, 2012, unpublished manuscript). In the next step, all parameter sets found for a given basin were transferred to the remaining ones. Finally, only those sets exhibiting a model efficiency greater than or equal to 0.65 at recipient locations were retained as members of the best global parameter sets. This implies that these super sets of global parameters are able to reproduce water fluxes in all major river basins in Germany with an efficiency of at least 0.65. It may be noted that a single set of VIC-3L model parameters for a large domain in the midwestern United States was used in a study by Mishra et al. (2010) for assessing historical drought events. In contrast to that, in this study the ensemble of 200 model realizations was used for the subsequent analysis of historical drought events in Germany including both uncertainty branches.
In general, mHM requires at least 5 years of spinup time to equilibrate. To minimize the influence of initial conditions, all state variables (e.g., water content at a given soil layer) in each ensemble member were initialized with their climatological averages corresponding to the precise time of year at the initialization (Rodell et al. 2005). The climatological average was estimated as the long-term mean of a given state variable within a 7-day window around 1 January. The DWD1 precipitation estimate was employed to estimate the long-term mean. This procedure allowed us to reduce the spinup time to 1 year without inducing large bias due to inappropriate initial conditions. Thus, model simulations during the starting year 1950 were discarded from the following analysis.
3. The mHM soil moisture index
The absolute values of the soil moisture states estimated with mHM do not allow a direct comparison of derived drought indices across the study domain because anomalies in absolute terms reflect climatological and morphological characteristics (Andreadis et al. 2005) rather than strong deviations from the respective normal conditions, which is the main characteristic that defines a drought event. Instead of absolute values, agricultural droughts can be quantified as “deficit of soil moisture relative to its seasonal climatology at a location” (Sheffield et al. 2004). The main idea behind this definition is to develop an index that varies between 0 and 1, which indicates drier to wetter conditions, respectively. The apparent selection for such an index is the conditional cumulative distribution function of the soil water content in the root zone at a given location i and time of the year m. This kind of normalization is inspired by the standardized precipitation index (McKee et al. 1993). The procedure to estimate an SMI based on mHM soil moisture simulations is described next.
a. Aggregation and normalization
Daily mHM soil moisture from three soil layers was averaged for every grid cell to obtain monthly states. These monthly values were, in turn, normalized with respect to the corresponding total root zone saturated water content (i.e., porosity times the total depth of the soil layers) to estimate the monthly soil moisture fraction (x) of the total soil column; namely,
where xl is the monthly soil moisture at root zone layer l (mm), and is porosity or the saturated water content of root zone layer l (mm). In the present study l = 3. In this case, the indexes i and m are omitted to ease the notation.
b. Estimation of the SMI
The monthly soil moisture fraction [Eq. (1)] may exhibit heavily skewed, nongaussian distributions (Koster et al. 2009) whose shape varies depending on climate and soil characteristics. The distribution of this random variable can also be multimodal (Vidal et al. 2010), which is an indication of preferential states of seasonal soil moisture (Rodriguez-Iturbe et al. 1991; Laio et al. 2002). Consequently, describing this random variable with unimodal theoretical distribution (e.g., the beta distribution; Sheffield et al. 2004) is not appropriate. Instead of making assumptions regarding the theoretical distribution of this variable, which would induce an additional source of uncertainty, a nonparametric technique was adopted to estimate the probability density function of the monthly soil moisture fraction at every cell within the domain, denoted hereafter as . The estimation procedure is as follows.
Given a set of data from one of the ensemble members x1, x2, … , xn that corresponds to the monthly soil moisture fractions of a given cell within the domain during month m (e.g., January), the kernel density estimate at a given value x can be obtained by
where K(x) is the smoothing kernel, n the sampling size, and h the bandwidth. The sampling size in this case is equal to 60. There are various possibilities to select K(x) (Wilks 2011); however, the Gaussian kernel is appealing in this case because of its unlimited support. The optimal selection of the bandwidth can be obtained by minimizing the unbiased cross-validation criteria (Scott and Sain 2005) given by
where is the leave-one-out density estimate at x when observation xk is not taken into account. This optimization was performed with a generalized reduced gradient algorithm. Once the optimal bandwidth is found, the best fit of the empirical distribution function can be estimated .
Finally, the mHM soil moisture index for a given cell and month, which denotes the quantile at the soil moisture fraction value x, can be obtained by numerically integrating the expression
c. Identification of drought events
Droughts are regional phenomena covering large contiguous areas over long periods. Understanding the spatial–temporal patterns and their relationships with other variables is therefore a fundamental step for drought prediction. Previous drought studies carried out in Germany, however, have been focused on statistical analysis of readily available point observations such as river discharge or precipitation data (Demuth and Heinrich 1997; Stahl and Demuth 1999; Franke et al. 2004; Schindler et al. 2007; Schindler and Mayer 2007), and in general, they are limited to a regional scale rather than to the national scale. To the best of our knowledge, studies investigating the spatial–temporal drought variability over the whole German territory are not available in the scientific literature.
The retrospective reconstruction of soil moisture analysis in Germany provides a unique dataset to estimate fundamental characteristics (e.g., severity and areal extent) of the major agricultural droughts occurred in Germany since 1950 at a high spatial resolution. Drought events were identified in this continuous spatiotemporal dataset with the method proposed by Andreadis et al. (2005).
First of all, regions under drought stress were identified with the threshold method (Dracup et al. 1980). This implies that cells fulfilling SMIt < τ were selected as potential regions under drought at the monthly time step t. The selection of the truncation level τ is fundamental for this method. A common value adopted in the literature is τ = 0.2 (Andreadis et al. 2005; Vidal et al. 2010). This threshold indicates that a given cell is enduring a soil water deficit occurring less than 20% of the time.
In the second step, drought clusters at every monthly time step have to be consolidated in space. This means that all clusters whose area is less than a minimum threshold area will be excluded from further analysis. This step is necessary to eliminate small isolated areas that are suffering a drought but are too small to be considered as a regional event. In this study the minimum cluster area was set to 640 km2 (i.e., 40 cells).
The final step of the drought event identification consists of consolidating independent spatial clusters over successive time steps into a regional, multitemporal cluster. This kind of clustering is necessary because the spatial variability of a drought event is vast, composed of many branches that can either merge together or split over time. The only condition to join clusters over time is that the overlapping area should be larger than 6400 km2 (i.e., 400 cells). Overlapping areas less than this threshold area was considered as independent drought events.
Both threshold areas (i.e., the minimum cluster area and the overlapping area) were determined though sensitivity analysis but primarily based on rules of thumb often followed in the literature (e.g., Andreadis et al. 2005; Vidal et al. 2010). The main criterion for the selection of these parameters was the stability of drought characteristics described in the following section. It should be noted that the selection of smaller areas, enduring drought conditions, leads to the proliferation of smaller clusters that are not contiguous over time and hence cannot be considered as part of a regional phenomenon.
d. Quantification of drought characteristics
Drought characteristics such as mean duration, mean areal extent, total magnitude, intensity, and severity–area–duration curves were quantified for every drought event and every ensemble member. The mean duration (D) of a spatiotemporal drought event is defined as the average of the drought duration of every cell within a drought event. This statistic is given in months. The mean areal extent (A) is defined as the average of a region under drought from the onset until the end of the drought event, expressed as percentage of the total German surface area. The total magnitude (M) is defined as the spatiotemporal integral of the SMI below the threshold value τ (i.e., the deficit) over those areas that are affected by the drought event, or explicitly
where t0 and t1 denote the onset and the ending months of a given drought event; At is the area under drought at a given time step t, expressed as the percentage of total German surface area; i denotes a given location within the domain At; and (·)+ is the positive part function. Thus, M is expressed in months times percentage of total German surface area.
The above-described three statistics are useful to rank drought events based on the overall impact but they do not allow us to estimate the impact of the drought after some months from the onset. This could be better quantified with the drought intensity (Id) at a given duration d from the onset of the event. This statistic can be estimated as
This statistic also allows us to estimate the impact of various events during summer and winter by discriminating the time step t0 + d to a corresponding season.
Another commonly used method to benchmark drought events is based on the severity–area–duration (SAD) curves proposed by Andreadis et al. (2005). The severity (Sd) for every cell for a given duration d in months can be estimated as
The SAD curves for durations of 3, 6, 9, and 12 months for a given ensemble realization were constructed as follows: Firstly, the grid cells were ranked according their severity. The procedure starts with those cells having the maximum severity. Then, the severities of the adjacent cells were summed up progressively until a threshold area is reached. Afterward, the average severity is estimated for those selected cells. The cumulative area and the average severity constitute the abscissas and ordinates of the SAD curves for a given duration. In this study, regular area intervals equivalent to the area of 20 grid cells were selected (i.e., every 320 km2). This procedure is repeated until the whole area of a given drought event is covered.
The monthly evolution of these statistics was estimated for every member of the ensemble. Based on the ensemble simulations, the uncertainty of the four selected statistics was analyzed.
4. Results and discussion
a. mHM evaluation
The performance of mHM was evaluated against observations of daily streamflow, latent heat, and soil moisture measured at various eddy covariance (EC) stations acquired from www.fluxdata.org as well as with soil moisture observations obtained with a cosmic ray neutron probe (Rivera Villarreyes et al. 2011). Seven large river basins in Germany were selected to cross validate mHM performance with respect to observed daily streamflow. In this proxy basin test, global calibration parameters of mHM obtained at every river basin were transferred to the remaining test basins [e.g., from Neckar to Danube, Main, Ems, Saale, Mulde, and Weser basins (Fig. 1)]. The procedure to find the best 100 global parameter sets is described in section 2b.
High efficiency in this kind of test is a good indication of model performance in ungauged locations. The ensemble mean Nash–Sutcliffe efficiency (NSE) obtained with mHM using the best 100 global parameter sets at proxy basins during the validation period from 1965 to 1999 varied from 0.50 to 0.88, which is quite acceptable considering that these basins have significantly different hydrologic regimes. Model evaluation on those basins with at-site calibrated parameter sets during the same period exhibited on average an NSE value ranging from 0.74 to 0.93. During the calibration period (2000–04), the NSE varied from 0.84 to 0.96. These tests indicated that mHM can be used for hydrological predictions within Germany.
The coefficient of determination between the simulated latent heat fluxes against observations across several EC stations varied between 0.50 and 0.74 during the period 2000–02. The model domain in this case was reduced to a cell size of 100 × 100 m2. Considering the various factors that influence the EC measurements and the fact that mHM is driven by disaggregated hourly values of precipitation and temperature as well as known scaling issues with EC measurements, these results can be regarded as satisfactory. The soil moisture anomalies estimated with mHM were able to explain up to 75% of the variance of their observed anomalies at various EC stations during the same period. Soil moisture estimates were obtained with standard time-domain reflectometry (TDR) probes. The model at the EC sites was forced with observed hourly precipitation and hourly temperature instead of the interpolated data as used for running the model over the whole domain.
The cosmic ray neutron probe, on the other hand, is a promising alternative because it allows an estimate of the soil water content over a control volume with a diameter of approximately 600 m and a depth of 0.3 m, which, in this case, corresponds to the tillage depth setup in mHM. The coefficient of determination (r2) between the mHM soil moisture anomaly and the cosmic ray neutron probe, reported by Rivera Villarreyes et al. (2011), was 0.57 for the period August–September 2011. Correspondingly, the r2 between the simulated and the mean of soil moisture anomalies measured with 16 frequency domain reflectometry probes located within the same control volume was 0.79.
b. Retrospective reconstruction of soil moisture fields
The basis for the analysis of agricultural drought analysis in Germany was the reconstruction of daily soil moisture fields since 1950. Two hundred realizations of these fields were estimated for the whole of Germany at an hourly basis based on the premise that a single simulation is not sufficient for such analysis because of parameter uncertainty. For the subsequent analysis, simulated hourly fields were aggregated to daily and monthly time steps. Monthly soil moisture values were then normalized as indicated in Eq. (1) to ease comparison across locations. The ensemble long-term mean for each month (Fig. 2) is the most evident statistic to evaluate these results and to verify whether the annual variability of soil moisture corresponds to the known climatology of major geographic regions in Germany. The variability of the spatial patterns shown in Fig. 2 indicates almost saturated conditions the whole year round in mountainous areas such as the Black Forest, the Harz Mountains, and the Bavarian Alpine Foreland. Quasi-permanent dryer conditions have been observed on the North German Plain. The variability of the long-term mean of the soil moisture fraction x with respect to its standard deviation indicates a clear seasonality describing wetter and less variable conditions in winter opposed to less wet but highly variable conditions in summer (Fig. 3).
Results indicated that not only the ensemble monthly climatology of the soil moisture fraction x, depicted in Fig. 2, but also other statistics such as the 10th and 90th percentiles of x (P10, P90) exhibit seasonality and strong dependency to geographic location. The annual variability of these two percentiles for selected cells within Germany is depicted in Fig. 4. The geographic location of the selected cells is shown in Fig. 1. The variability and the value of both percentiles indicate marked hydroclimatic regimes in Germany—for instance, humid regions with moderate seasonality on the North Sea (cells 1 and 2), very humid regions with very little seasonality in the alpine regions (cells 18–20), very humid regions with marked seasonality on the Black Forest (cells 13 and 17), moderately dry regions with marked seasonality in the North German Plain (cells 7 and 8), and regions with large seasonality in the prealpine regions (cells 14 and 15). In general, the standard deviation of the 90th percentile of x is less than that of the 10th percentile based on the 200-member ensemble. This corroborates the findings of Schaake et al. (2004) and Meng and Quiring (2008), which point out that the parametric uncertainty in drier regions (cells 7, 8, 11, and 15) is much higher than in humid regions (cells 17–20). The standard deviation of both percentiles exhibits not only seasonal variability, clearly depicted in cell 15 shown in Fig. 4, but also strong geographic dependency. This indicates that there is a complex interplay between climatic conditions and parameterization of the soil moisture processes.
c. Comparison with other indices
The same method proposed in section 3 to estimate the SMI can be used to estimate drought indices based on precipitation and surface runoff generated at each cell before routing (Shukla et al. 2011). The results of these three drought indices are shown in Fig. 5 for one of the ensemble realizations obtained with DWD1. The three upper panels of this figure indicate how different the spatial distribution of the drought index might become depending on the variable used to describe a drought event. Among the three variables, the drought index based on precipitation exhibits the largest spatiotemporal variability because of the lack of memory of the precipitation process, which is one of the main reasons for considering it not appropriate for describing water stress in vegetation (Döring et al. 2011). The drought index based on surface runoff is correlated to the SMI but still quite weak because of fast runoff generation processes. The SMI, as compared with the other two indices, exhibits the largest persistence.
d. Sensitivity of the parameter uncertainty related to precipitation interpolation
Among the two sources of parametric uncertainty investigated in this study, the first one was related to the interpolation methods used to regionalize rainfall point data. For this purpose, two methods were employed to estimate the gridded fields of precipitation data, as denoted by DWD1 and DWD2 (see section 2b for details). Since both methods use the same input data, any possible variation in soil moisture simulations—ceteris paribus—could be attributed to the kriging weights and the variogram parameterization used in DWD1 or the linear weights of the multilinear regression method employed in DWD2. In this respect, two questions were pursued in this study: 1) How important is this source of uncertainty for the estimation of soil moisture? And, 2) how is this uncertainty distributed over space?
To answer these questions, the Pearson correlation coefficient (r) of the monthly soil moisture fractions at every grid cell obtained with both precipitation products (i.e., DWD1 and DWD2) were estimated separately for all 100 global parameter sets. From these r values, the ensemble mean and the coefficient of variation of r were calculated for every cell within the domain. These statistics are depicted in Figs. 6a and 6b, respectively.
In general, most of the grid cells within Germany exhibit a value greater than 0.98, which indicates a high degree of agreement between any pair of simulations driven by DWD1 and DWD2 forcings but having the same global model parameters. There are very few places where this statistic is less than 0.98, but in every case greater than 0.95. This finding along with the very low coefficient of variation indicated a quite low sensitivity of the monthly soil moisture fraction to the precipitation interpolation parameters. The lower values of were obtained mainly in cells located in and around mountainous regions such as the Harz, the Alps, and the Swabian Jura (Fig. 6).
e. Overall parameter uncertainty of the SMI
The two major sources of parametric uncertainty described above induced considerable variability into the SMI as shown in Fig. 7, which depicts the areal average of the SMI over major German river basins, denoted hereafter as 〈SMI〉. It can be noticed from this figure that the overall parameter uncertainty of 〈SMI〉 is neither constant in space nor over time. The 〈SMI〉 obtained with each ensemble member exhibited a large variability within the interquartile range of SMI but a relatively small one at its extreme quartiles (Fig. 7). This behavior is closely related with the high variability of the standard deviation of the soil moisture fraction around the middle ranges of its mean value (e.g., between 0.6 and 0.8 as depicted in Fig. 3).
For further analysis, the temporal variability of 〈SMI〉 within the ensemble simulations is estimated by its range R(t) = 〈SMI(t)〉max − 〈SMI(t)〉min, at every point in time t. Here R(t) denotes the ensemble uncertainty of the soil moisture index over a given domain at time t. The long-term average of R(t) is approximately 0.124 with a standard deviation of 0.014. The correlation coefficient estimated between the range of time series R(t) for every pair of major basins, depicted in Fig. 1, varied from 0.25 to 0.88. This implied that the uncertainty of the SMI is not only the result of independent errors arising from model parameterization, but also the result of systematic interdependencies between soil moisture and climatic variables such as precipitation (P) and potential evapotranspiration (Ep). Based on these results, it was determined that the standard deviation of R(t) tends to decrease as the ratio Ep/p increases. Moreover, given the data provided for each major basin (Fig. 7), the null hypothesis that the time series of the ensemble uncertainty R(t) constitutes white noise can be safely rejected provided that the p value of the Fisher’s kappa statistic was less than 0.001.
The 12-month moving average of 〈SMI〉 depicted in Fig. 8a over the reconstruction period (1951–2010) showed a considerable reduction in uncertainty compared with the monthly values of 〈SMI〉 but still not small enough to be considered negligible. The 12-month moving average of the percentage of area under drought (with respect to the surface area of Germany) exhibited a considerable uncertainty at the peaks of the events (Fig. 8b). This result, however, allows preliminary identification of major drought events covering at least 50% of the German territory, namely those in the periods 1953/54, 1959/60, 1964/65, 1972/73, 1976/77, 1992/93, and 2003/04.
The parametric uncertainty of the SMI also has a strong influence on drought severity classes commonly used for monitoring purposes. Figure 9a depicts the probability of finding a cell, at a given point in time, under one of the five drought severity classes used by the United States Drought Monitor (http://droughtmonitor.unl.edu). These classes denote abnormal (D0), moderate (D1), severe (D2), extreme (D3), and exceptional (D4) dry conditions, which correspond to 0.2 < SMI ≤ 0.3, 0.1 < SMI ≤ 0.2, 0.05 < SMI ≤ 0.1, 0.02 < SMI ≤ 0.05, and SMI ≤ 0.02, respectively. This figure shows also that there are areas in which no unique drought class can be assigned because of parametric uncertainty. A possibility to assign a unique class to a cell is to choose a class with the largest probability, as shown in Fig. 9b for May 1976.
f. Identification of major drought events based on mean duration, mean areal extent, and total magnitude
Major drought events were found in this study using the technique described in section 3c. These benchmark events are required for the future analysis of possible consequences of climate change on agricultural droughts. The drought clustering algorithm was applied to every ensemble realization to find the spatiotemporal evolution of all drought events during the reconstruction period from 1951 to 2010. For every event, drought characteristics such as mean duration (D), total magnitude (M), and mean areal extent (A), among others, were evaluated using the procedure illustrated in section 3d. The ensemble average of these characteristics (i.e., , , and ) are depicted in Fig. 10. The corresponding uncertainty of these characteristics is presented in Table 1.
The eight largest drought events identified during the last 60 years in Germany are the following periods: 1962–65, 1971–74, 1975–78, 1959/60, 1953/54, 1991–93, 2003–05, and 1995–97. It is worth noting that the event from 2003 to 2005 appears in this overall ranking in the seventh position. Vidal et al. (2010) also noticed this fact and concluded that 2003 hardly appears as a benchmark event in France. This is a rather controversial conclusion because in this year the highest temperatures during the last 500 years were recorded (Luterbacher et al. 2004). In Germany alone, great losses in the agricultural sector (COPA-COGECA 2003) were reported. A likely explanation for this paradox is provided in section 4g.
The three drought characteristics D, M, and A depicted in Fig. 10 are highly correlated with each other. The Pearson correlation coefficient between D and M is the highest, and equal to 0.97, whereas those between (D and A) and (M and A) are 0.80 and 0.87, respectively. This indicates that this triplet has low dimensionality. In fact, the first eigenvector of the correlation matrix of this triplet alone explains 92% of the total variance.
Using the k-means cluster analysis, three main groups of drought events were distinguished: 1) events with a large areal extent and duration (i.e., events 1962–65, 1971–74, 1975–78, and 1959/60), 2) events with the largest areal extent and moderate duration (i.e., 1953/54), and 3) events with moderate areal extent and duration (i.e., 1991–93, 2003–05, and 1995–97). Based on the ensemble SMI mean , the event from 1971 to 1974 exhibited the longest duration, and the event from 1953 to 1954 covered the largest area. The events from 1962 to 1965 and 1971 to 1974 reached the two largest magnitudes.
The absolute ranking of these extreme drought events is rather difficult because of the parameter uncertainty as illustrated in Table 2. This table presents an estimate of the probability to order every event into the eight top ranks using a linear, equal-weighted, normalized indicator composed of D and A, as an example. The results presented in this table indicate that the maximum probability of finding an event in one of the top ranks is not greater than 0.67. The ranking of a given event spans at least over three categories. Low-ranking events tend to have a much larger ranking spread than the top ones, though.
The size of the ensemble also played a very important role to estimate the probability of finding an event in a given rank (1 − α), where α denotes the false positive rate. Figure 11, for example, shows the probability of not identifying the event from 1971 to 1974 as the largest since 1951. This figure clearly shows that the variance of the false positive rate is strongly dependent on the ensemble size. These results were obtained by bootstrapping the 200 ensemble simulations without replacement and limiting the number of realizations to 1000 for a given sample size. This figure showed also that the first two moments of α tend to stabilize with ensemble sizes larger than 50. Consequently, it is safe to conclude that small ensemble sizes would lead to misleading results. An ensemble with 200 members, as realized in this study, would lead to safer results. These Monte Carlo realizations clearly highlighted the role of parametric uncertainty in identifying the benchmark drought events that should be handled carefully.
The spatial distribution of severity (Sd) based on at the peak of the eight largest drought events is shown in Fig. 12. It can be observed from this figure that each event has its own peculiarities with respect to the spatial distribution of the affected areas. The drought event during December 1954 has the largest areal coverage, with 93.5% of the German territory under water stress, whereas the event during April 1996 had the lowest coverage with 46.5%. The latter drought event at its peak was particularly concentrated on the northwest part of Germany. The event of 1976, with its peak in August, had spread over the whole of Germany with the exception of the Alpine Foreland. The latter areas endured the highest severity during August 2003.
g. Uncertainty of large events occurring in summer and winter
As mentioned before, the ranking of drought events based on ensemble characteristics (D, M, and A) does not allow the identification of their impact at a given point in time from their onset nor to differentiate them according to their level of incidence in a particular season. The drought intensity proposed in Eq. (6) enables the estimation of the transient evolution of a drought event from its onset, and by so doing, it allows quantifying how fast a drought event covered a given area and by what magnitude. Figure 13a shows the results of plotting drought intensity versus duration from the onset (d) of a given event for the 10 largest events since 1950. Figure 13b depicts the results obtained by ranking the drought intensities of all events at various durations from their onsets (e.g., 3, 6, … months). The classification of an event into summer or winter was estimated with the procedure illustrated in section 3c [Eq. (6)]. The ensemble SMI mean (i.e., ) was used instead of individual realizations for both analyses because the former is an unbiased estimate of the SMI, and thus leads to a robust estimate of the evolution of the drought intensity.
Based on the results described above and shown in Fig. 13, it was found that at a 3-month duration, summer events have much larger drought intensity than the corresponding ones in winter. At a 6- and 9-month duration, the opposite happens. The events with more than a 9-month duration mostly reach their higher intensities during summer as compared to winter ones. However, droughts having a duration of 30 months or more are more intense during winter months. The event 1953/54 not only exhibits the largest intensities at 6- and 9-month durations during winter months (November–April) but also the largest intensity in summer at a 12-month duration. The event 2003–05 is, according to these results, the summer event with the largest intensity at a 6-month duration. Among the 10 largest drought events in Germany during last 60 years, the 1953/54 event had the largest intensity, peaking within a relatively short period of time (less than 12 months). This event, however, lasted for only 1.5 years. Four drought events (i.e., 1962–65, 1971–74, 1975–78, and 1991–93) spanned a period of more than 30 months (i.e., 2.5 years). According to this analysis, the decade of 1970 could be regarded as the most severe drought period in Germany. The drought events 1962–65 and 1971–74 clearly exhibited more than one peak over their whole life span. The analysis also indicated that most of the historical drought events in Germany have their peaks during 6 to 12 months of duration.
The empirical bivariate density function between the average drought area (A) and the total magnitude (M) was constructed to analyze the uncertainty in overall drought characteristics (D, M, and A) based on the ensemble realizations. The large number of model runs also allowed us to assess the uncertainty in time evolution of these characteristics. The four most intense drought events with 6 months and at least a 30-month duration after its onset were selected to illustrate this procedure—namely, the events 1953/54 and 2003–05 for shorter duration, and the events 1975–78 and 1962–65 for longer duration (Fig. 14). It is worth noting that the events 1953/54 and 2003–05 are classified as winter and summer events, respectively, at a 6-month duration (Fig. 13b). Likewise, the events 1975–78 and 1962–65 peaked in winter and summer, respectively. Droughts that are peaking within a relatively short time (up to 6 months) from their onset are quite relevant because they have large repercussions on socioeconomic activities.
Based on the ensemble results, the density function for each event was estimated independently with a bivariate Gaussian kernel smother algorithm. The estimation of the bandwidths in both directions was carried out in a similar way as presented in section 3b. The results of this analysis are depicted in the top panels of Figs. 14a–d, which clearly supports the research hypothesis that the parametric uncertainty of soil moisture has a strong implication for drought characterization. Most events exhibit multimodal behavior, which is the combined result of the uncertainty of the model parameterization and drought identification (e.g., clustering and threshold).
Events having shorter durations and peaking in winter (1953/54) appear to be more certain than those peaking in summer (2003–05) as can be noted by the larger spread of the respective distribution (Figs. 14a,b). Consequently, the probability density values for the summer event are lower than those of the winter event. However, at longer durations no conclusive comparison could be drawn from this analysis because longer events experience various seasons over many years. The time evolution of the area under drought A(t) for each events, as depicted in bottom panels of Fig. 14, also supports the assertion that a single model realization would lead, very likely, to a high rate of false alarms for drought monitoring.
h. Uncertainty of the severity–area–duration curves
SAD curves obtained with the ensemble SMI mean for the eight largest drought events in Germany at duration 3, 6, 9, and 12 months are depicted in Figs. 15a–d. From this analysis, the event from 1975 to 1978 appears to be the most severe and extensive event at durations ranging from 3 to 9 months. Based on this measure, the 2003–05 event, however, hardly appears as a benchmark event at longer durations and area coverage. The event from 1953 to 1954 is quite severe at 3 and 6 months, but not at longer durations. The apparent contradiction of these results can be clarified with the individual evolution graphs presented in Fig. 13.
SAD curves have often been used to rank drought events (Andreadis et al. 2005; Sheffield et al. 2009). Because of parametric uncertainty, however, they exhibit large variability as shown in Fig. 15e. This, again, corroborates our hypothesis that a single model run would lead to unsatisfactory conclusions and event ranking. These results indicate that the SAD variability increases as the area under drought and duration increase. The variability of the SAD curve with a 12-month duration is almost twice as much as that for 3 months. The variability of SAD curves for summer events is higher than that estimated for winter at any duration.
i. Drought persistence and SMI trends
Characterizing areas prone to remain under severe drought conditions when they are already suffering one constitute a relevant piece of information for water resources planning. The level of persistence of the severe drought events can be quantified with a two-state Markov chain with two states: SMI ≤ 0.2 and 0.2 < SMI ≤ 1. The persistence of severe drought can be estimated for each ensemble member as the probability π00 = Pr[SMI(t + 1) ≤ 0.2 | SMI(t) ≤ 0.2], for all t. The ensemble mean of π00 is depicted in Fig. 16a for the whole of Germany. This figure indicates that most of the Northeast German Plain—comprising the area of the Elbe, Saale, and Mulde River basins, as well as large extensions along the Main and Rhine Rivers—exhibits drought persistence greater than 0.8. The Northwest German Plain, comprising the Ems and Weser River basins, tends to have lower drought persistence than the eastern part of Germany, with an average value of π00 less than 0.7. The Alpine Foreland located within the Danube basin and areas in and around the Black Forest, on the contrary, exhibit the largest variability in drought persistence within Germany ranging from less than 0.4 to 0.8. It is worth nothing that those areas exhibiting large drought persistence have been also classified as areas with medium to high agricultural suitability according to a recent study conducted by Umweltbundesamt–Potsdam-Institut für Klimafolgenforschung (UBA-PIK; Research Report 20141253 UBA-FB 000844/e, available online at www.umweltbundesamt.de). These regions comprise large plains within the Saale River basin around the cities of Halle and Magdeburg and flood plains of the Rhine River on the western side of the Black Forest.
Mann–Kendall tests on monthly SMI indicate that there are large extensions of the German territory showing positive trends (i.e., getting wetter) during winter months but negative trends in summer months at a 5% significance level. The largest areas exhibiting significant trends were detected in March and August as depicted in Figs. 16a and 16b, respectively. It is worth noting that positive SMI trends tend to occur in areas with low persistence and negative trends in areas with high persistence. These trends are, in turn, related, with observed trends in temperature and precipitation. Further details on this aspect are beyond the scope of this paper.
5. Summary and conclusions
In this study we have presented a method to derive a soil moisture index based on a process-based hydrological model. This model uses a multiscale parameterization method that goes beyond standard calibration approaches. Great emphasis has been put on testing this model in all major river basins in Germany, especially with respect to the transferability of global parameters across locations and scales. Ongoing tests with the global flux network (FLUXNET) and cosmic ray neutron probe data have also been presented. Using this model a consistent ensemble of high-resolution daily soil moisture fields for Germany since 1950 at a spatial resolution of 4 × 4 km2 were obtained.
Based on this soil moisture reconstruction, a soil moisture index (SMI) representing the corresponding monthly quantile was estimated with the kernel density approach. The derived SMI exhibits high correspondence with total grain yield of Germany and allows us to identify major drought events in Germany that have also been identified using other techniques (e.g., tree rings) and reported in the literature (Büntgen et al. 2010). This approach has advantages over standard empirical approaches or those obtained from satellite-derived products, which are too coarse to account for soil moisture at high spatiotemporal resolutions and quite uncertain because the algorithms used to infer soil moisture do not take into account the water balance of large river basins. Consequently, the proposed technique has a large potential to be used as a monitoring tool in the future. More research is, however, needed to evaluate the SMI against times series of annual crop yield at regional scale. Further research is also required to identify potential driving mechanisms, the feedback effects, and the spatiotemporal correlations of soil moisture with other hydrological state variables such a snow depth and climatic variables.
The effects of other sources of uncertainty stemming from model structure and quality of meteorological data on the soil moisture index should be further investigated. Potential benefits of using ensembles of multimodel, multiparameter soil moisture simulations should be also carried out. Both issues, however, are out of the scope of this study.
Based on the results of this study, the following conclusions were drawn: 1) The main source of parametric uncertainty of the soil moisture index is related with global model parameters. This uncertainty is seasonally and regionally varying. This corroborates findings of other researchers who have advocated for multimodel ensembles to account for model uncertainty. In summary, one single model run is not enough for estimating benchmark events. 2) The uncertainty of overall statistics used for estimating drought events are highly sensitive to this kind of uncertainty. This sensitivity is the result of nonlinear relations and branching effects caused by the clustering method. 3) Events peaking during summer with at most 6 months duration tend to exhibit a much larger uncertainty than those peaking during winter. 4) The SMI is not a stationary variable. Many regions in Germany exhibited significant trends during the study period. Potential triggering mechanisms and drivers behind these trends might be the observed changes of precipitation and temperature as well as other feedback mechanisms. A detailed trend attribution, however, is out of the scope of this study. 5) The identification of benchmark drought events should be based on combined criteria such as SAD or intensity duration curves. Robust estimates can only be made with an ensemble SMI because of the uncertainty mentioned before.
This study was carried out within the context of the REKLIM project. The authors would like to acknowledge the collaboration of the Federal Institute of Hydrology (BfG) supporting the European Water Archive (EWA) and the Global Runoff Data Centre (GRDC), the Federal Institute for Geosciences and Natural Resources (BGR), the Federal Agency for Cartography and Geodesy (BGK), the German Weather Service (DWD), the European Environment Agency (EEA), FLUXNET, and the NASA Land Processes Distributed Active Archive Center (LP DAAC) for providing the data employed in this study. We are also very grateful with Mr. Rivera Villarreyes for kindly providing cosmic ray neutron probe data for validating our model. We would also like to thank Ms. Elisabeth Krüger, Mr. Pablo Beltran, and two anonymous reviewers for their valuable comments.