Climatic Inﬂuences on Agricultural Drought Risks Using Semiparametric Kernel Density Estimation

: A bivariate kernel density estimation (KDE) method was utilized to develop a stochastic framework to assess how agricultural droughts are related to unfavorable meteorological conditions. KDE allows direct estimation of the bivariate cumulative density function which can be used to extract the marginal distributions with minimal subjectivity. The approach provided excellent ﬁts to bivariate relationships between the standardized soil moisture index (SSMI) computed at three- and six-month accumulations and standardized measures of precipitation (P), potential evapotranspiration (PET), and atmospheric water deﬁcit (AWD = P − PET) at 187 stations in the High Plains region of the US overlying the Ogallala Aquifer. The likelihood of an agricultural drought given a precipitation deﬁcit could be as high as 40–65% within the study area during summer months and between 20–55% during winter months. The relationship between agricultural drought risks and precipitation deﬁcits is strongest in the agriculturally intensive central portions of the study area. The conditional risks of agricultural droughts given unfavorable PET conditions are higher in the eastern humid portions than the western arid portions. Unfavorable PET had a higher impact on the six-month standardized soil moisture index (SSMI6) but was also seen to inﬂuence three-month SSMI (SSMI3). Dry states as deﬁned by AWD produced higher risks than either P or PET, suggesting that both of these variables inﬂuence agricultural droughts. Agricultural drought risks under favorable conditions of AWD were much lower than when AWD was unfavorable. The agricultural drought risks were higher during the winter when AWD was favorable and point to the role of soil characteristics on agricultural droughts. The information provides a drought atlas for an agriculturally important region in the US and, as such, is of practical use to decision makers. The methodology developed here is also generic and can be extended to other regions with considerable ease as the global datasets required are readily available.


Introduction
Agricultural droughts have resulted in billions of dollars of losses in recent years and have caused other socio-ecological impacts [1,2]. In addition to devastating fragile rural economies, agricultural droughts have the potential to disrupt global food security [3]. Droughts are also known to decrease structural sugars and lignin content and thus affect the ethanol yields in biorefining operations [4]. Thus, drought-related water deficits can also impact renewable energy production. Therefore, understanding agricultural droughts is of paramount importance for sustaining rural economies and fostering food security and energy independence.
Agricultural droughts arise when the soil moisture is affected by enhanced atmospheric dryness brought forth by precipitation deficits and associated temperature increases over an extended period. This latter situation is referred to as meteorological drought. Soil moisture deficits affect plant growth, focused on the use of stochastic frameworks to elucidate how changes in meteorological forcings affect agricultural drought risks and not on developing a joint drought index perse. The developed framework is presented next and illustrated by applying it to the US High Plains region overlying the Ogallala Aquifer to understand factors affecting agricultural droughts in this agriculturally important region.

Soil Moisture Modeling
While the proposed stochastic framework is model agnostic, the leaky bucket model [15,22] is adopted here for illustration. This model is known to provide suitable soil moisture predictions for quantifying agricultural droughts and is widely used for monitoring soil moisture anomalies [23]. Furthermore, the model, while built using the fundamental principle of mass balance, is also parsimonious and can be run with readily available data. As such, it provides a suitable conceptualization to illustrate the developed methodology.
Details of the leaky bucket model can be found in Huang et al., (1996) [15]. Briefly, the model evaluates the soil moisture dynamics in the upper 1-2 m (root zone) of the soil over an area A. The basic governing equation is given as: where W is the soil moisture analogue (in mm) aggregated over the area A and root zone depth; P, E, R, and G represent precipitation, evapotranspiration, runoff, and deep percolation fluxes, respectively, and take the units of mm/month. E, R, and G are assumed to be a function of soil moisture W. The relationships between the hydrological fluxes and soil moisture are stated using four parameters-µ,α,m,W max where W max is the maximum soil moisture that can be held in the root zone and µ,α,m are essentially fitting parameters. The actual evapotranspiration (E) is also a function of potential evapotranspiration (PET), which is computed as a function of temperature using the Thornthwaite PET model. These unknown model coefficients must be obtained via calibration. Huang et al., (1996) provide estimates for these parameters based on calibration studies carried out using watersheds in Oklahoma [15]. The calibrated model was shown to possess reasonable skill for use in drought monitoring studies based on comparisons from other locations [22], and is particularly suited for use in the illustrative case study, which also shares similar geographic characteristics to which the model was calibrated.

Standardized Soil Moisture Index (SSMI)-An Agricultural Drought Indicator
The soil moisture W computed using the model can be averaged over a time period T and used to calculate the agricultural drought index (referred to as the standardized soil moisture index or SSMI here) using the same steps that are used to calculate the standardized precipitation index, SPI [9,12,21]. For example, the time period T can be taken as three months and six months to represent intra-season and full season droughts [10]. Once the SSMI is computed, the severity (S) of the drought indicator can also be obtained via summation over a fixed duration (D) of interest (e.g., annual, summer, winter). The intensity (I) can then be calculated as the ratio of the severity (S) over (fixed) duration (D). While this approach does not give the actual duration of a drought event, it is nonetheless useful to simultaneously evaluate droughts of different severities or intensities over a duration of interest [32]. As growing seasons can generally be defined as a set of specific months, this approach of evaluating drought severity and intensity over a specified duration is advantageous in agricultural applications and as such adopted here. Meteorological forcings responsible for agricultural droughts can be obtained by summing up the atmospheric water fluxes over the adopted duration (D).

Relationship Between the Agricultural Drought Indicator (SSMI) and Meteorological Fluxes
Let P D and PET D be the precipitation and potential evapotranspiration fluxes aggregated over a duration D. The atmospheric water deficit (AWD) can then be written as [33]: In this study, P D , PET D , and AWD D are used as indicators of meteorological conditions. To facilitate a consistent comparison across space and time, these two indicators are standardized in a manner analogous to the SSMI. These standardized indicators are like the standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI) drought indicators, but only account for the total precipitation, potential evapotranspiration, and atmospheric water deficit that are summed over a duration of interest (i.e., summer, winter, annual) and unlike the SSMI are not accumulated using a rolling average of three and six months because, while the SSMI is known to show greater persistence over time due to moisture stored in the soil, the atmospheric conditions are assumed to be more volatile as the storage in the atmospheric compartment is not usually significant over a parcel of land.
Let S T,D be the severity of the SSMI corresponding to an averaging time period of T and summed over a duration of interest D. Let F(P D ), F(AWD D ), and F(S T<D ) denote the marginal distributions of standardized precipitation and standardized atmospheric water deficit (AWD) and severity of SSMI T , respectively. Similarly, F(S T,D , P D ) and F(S T,D , AWD D ) represent the joint distributions between the agricultural drought severity and meteorological fluxes-precipitation and atmospheric water deficit-computed over time T.
For a given set of cutoffs, the relationship between marginal and joint distributions of drought severity and meteorological flux (M D = {P D , AWD D }) can be specified using the relationships presented in Table 1. Similarly, the conditional distribution between agricultural drought severity given some information about the meteorological flux can be obtained using the Kolmogorov rule: For example, the cut-offs (−0.8: −1.2 (moderate drought); −1.3: −1.5 (severe drought); −1.6: −1.9 (extreme drought) and ≤ −2.0 (exceptional drought)) were adopted here based on recommendations of the US Drought Monitor (USDM) for meteorological droughts [23].

Estimation of Marginal and Joint Distribution Functions
The estimation of marginal and joint distributions is an important aspect of stochastic risk assessment. There are several approaches developed in the literature to model these distributions. Parametric distributions (e.g., the Gumbel logistical model) have been adopted in some flood frequency studies [34]. Copula-based approaches have been widely used to develop joint distributions from the knowledge of marginals [35]. Empirical joint distribution functions have also been used to model droughts [9]. Kernel density estimation (KDE) is used here to model the required joint distributions. KDE is a nonparametric method that has found usage in hydrology to model marginal [29] as well as joint distributions [30]. KDE procedures offer greater flexibility and accuracy in modeling distributions, especially the tail behavior. In addition, KDE does not make assumptions about the underlying distributions. It is always difficult in hydrology to know a priori the underlying distribution of the data, and the final distribution is selected from a small set of candidate distributions which may not be optimal [29,36]. The KDE approach can be used to directly fit cumulative distribution functions [37], which can overcome difficulties that may be associated with numerical integration of certain probability density functions. Rather than fitting two separate marginals and a copula for joint distribution, the cumulative distribution function (CDF) of the joint distribution can be directly integrated out, thus minimizing the subjectivity and error propagation during the estimation process. KDE also extends the range over which the distribution is fit and therefore is capable of extrapolating values out of the range of observations, much more so than empirical plotting position formulas. Therefore, KDE offers several advantages over other existing methods used for obtaining joint and marginal distributions.
The basic idea of KDE is to approximate the probability density function as a summation of a set of kernel functions. Mathematically for a d-variate dataset, the KDE of the probability density function (PDF) can be written as: where n is the number of data points, K is a kernel density function, H is the bandwidth matrix (which need not be a constant), x = X = (x 1 , x 2 , . . . , x d ) T are the scale parameters, and X i is a d-dimensional data vector (X 1,I , X 2,I , . . . , X d,i ) T given i = 1, 2, . . . , n. The kernel K is a symmetric probability density function, usually taken to be Gaussian, and the bandwidth matrix H is symmetric and positive-definite. Parameter estimation of the kernel density functions is carried out by minimization of the mean integrated squared error (MISE). The MISE is the sum of the squared bias and variance. The MISE depends upon the selected density (typically Gaussian) and the chosen bandwidth (which need not be fixed) [38].
The cumulative distribution function (CDF) can also be written in a manner analogous to the KDE as [39,40]: where the function W H satisfies the following property: According to Equation (5), W H is the cumulative distribution function corresponding to the probability density function K H .
The estimation procedures detailed in Duong (2007) [41] are adopted here, as it can directly estimate both the marginal and joint cumulative distribution functions necessary for risk analysis. This approach also considers the off-diagonal matrix elements in the bivariate variance-covariance matrix, which improves the ability to capture the correlation structure between the random variables. The Wand and Jones (1994) [42] plug-in estimator was used to guide the bandwidth selection process.
The marginal CDF of individual random variables can be obtained from the joint CDF. For the bivariate case (X = X 1 and X 2 ), the relationship between marginal and joint CDF is as follows: F(X 1 ) = F x1,x2 (X 1 , ∞) and F(X 2 ) = F x1,x2 (∞, X 2 ) where F is the marginal univariate CDF and F x1,x2 is the bivariate joint distribution of X 1 and X 2 .
Kernel density estimation can be carried out using software packages in R [43][44][45]. Once the distribution functions are estimated using KDE, visual comparisons with empirical bivariate CDF [46] and other goodness of fit measures can be calculated [47] to evaluate the reasonableness of the distribution functions.
A flowchart depicting the workflow of the proposed risk assessment methodology is presented in Figure 1.
Water 2020, 12, x FOR PEER REVIEW 6 of 28 A flowchart depicting the workflow of the proposed risk assessment methodology is presented in Figure 1.

Illustrative Case Study
The Ogallala Aquifer is the largest aquifer in the US and spans over an area of nearly 450,658 square kilometers (174,000 square miles) extending over eight states. The area overlying the aquifer (see Figure 2a) is predominantly rural and produces about 30% of the nation's agricultural output [48]. The region is comprised of five different climate types ranging from arid to humid condftabeitions ( Figure 2b). There is a pronounced precipitation gradient moving east to west, and the average annual precipitation declines from nearly 750 mm/year to less than 340 mm/year (Figure 2c). A temperature gradient can be seen moving along the north-south transect of the study area ( Figure 2d).
Water 2020, 12, x FOR PEER REVIEW 7 of 28

Illustrative Case Study
The Ogallala Aquifer is the largest aquifer in the US and spans over an area of nearly 450,658 square kilometers (174,000 square miles) extending over eight states. The area overlying the aquifer (see Figure 2a) is predominantly rural and produces about 30% of the nation's agricultural output [48]. The region is comprised of five different climate types ranging from arid to humid condftabeitions ( Figure 2b). There is a pronounced precipitation gradient moving east to west, and the average annual precipitation declines from nearly 750 mm/year to less than 340 mm/year ( Figure  2c). A temperature gradient can be seen moving along the north-south transect of the study area ( Figure 2d).  Portions of the Ogallala Aquifer have undergone severe depletion [49], which has increased the reliance on precipitation by farmers. The region is known for its erratic climate characterized by hot summers and mild to cold winters. Droughts are a recurring phenomenon in this region and their effects are often exacerbated by lack of suitable drought adaption [50]. Groundwater production from the already stressed Ogallala Aquifer is also known to increase during droughts [51]. Understanding how droughts propagate through agricultural production systems is not only important to develop suitable short-term contingency measures, but also for promoting effective water conservation programs that prolong the useful life of the aquifer. Sustaining the Ogallala Aquifer is not only important to the viability of the rural economies of the region, but also is critical for ensuring food security of urban areas that depend on it [3]. Therefore, the Ogallala Aquifer provides an excellent test bed to explore how precipitation and temperature impact agricultural droughts.
Precipitation, temperature, and PET data necessary for calculating AWD (Equation (2)) were obtained from the Climate Research Unit (CRU, version 4.03), University of East Anglea, Norwich, UK from 1901 to 2018. This dataset has been widely used in many climate related studies, and has been developed using several thousand weather stations across the world [52] and is available on a 0.5 • × 0.5 • grid which yields 187 locations across five different climate zones within the study area. PET data in the CRU dataset are based on the modified Penman-Monteith equation and account for both thermal and wind effects. The standardized precipitation and atmospheric water deficit were computed using this dataset. Please refer to Table S1 in Supplementary Material for additional details on data and their sources used in this study.
The leaky bucket model of Huang et al., (1996) [15] and van der Dool et al., (2003) [22] was extended to cover the entire period of interest and used to compute monthly SSMI values for three-month (intra-seasonal) and six-month (inter-seasonal) accumulations. To be consistent with the original formulation, the PET values (Equation (1)) were computed from temperature and latitude data using the Thornthwaite equation. The computed monthly SSMI values were then summed over summer (April-September), winter (October-March), and annual timescales to obtain agricultural drought severities at each station.
Joint KDE distributions were fit separately at each station, which were then used to obtain marginal distributions using Equation (6) and the procedures outlined by Duong, (2007) [41]. The one-dimensional and two-dimensional Gringorten plotting position formulas (Gringorten, 1964 [53] and Yue et al., 1999 [46]) were used along with the Kling-Gupta Efficiency (KGE) measure [47] to evaluate the KDE fits. All computations were performed within an R statistical and programming environment [45].

Evaluation of KDE for Fitting Joint and Marginal Distributions
The goodness of fit of bivariate joint distributions among standardized precipitation (P), standardized PET (PET), and standardized AWD (AWD) were evaluated using exploratory data analysis (EDA) methods and KGE metrics. Figure 3 presents a comparison of empirical and theoretical joint probabilities and KGE metrics for the SSMI3-P pair at six representative stations in five different climate zones in the region (see Figure 2b for station locations). Plots for other bivariate fits between SSMI3 and PET and AWD and SSMI6 with P, PET, and AWD are presented in Supplementary Information in the interest of brevity (see Figures S1-S17). Tables 2 and 3 present the goodness of fit metrics for SSMI and climate variables used in the study.
Water 2020, 12, x FOR PEER REVIEW 9 of 28  These plots and KGE metrics presented in Tables 1 and 2 capture and bracket the range of variability noted within the study area. As marginal distributions are directly obtained from bivariate CDFs, a good fit of the bivariate distribution is also an indication of a good fit of the marginals. The P-P plots (not presented in the interest of brevity) and KGE statistics for marginals summarized in Tables 2 and 3 indicate that the fits to marginals of both agriculture drought indicators and standardized meteorological variables were also excellent, with all KGE values being greater than 0.90 and most of them greater than 0.925. The correlation, bias, and variability terms were also reasonable, typically within ±5% of ideal values, and showed no systematic errors. The KDE models also exhibited a high degree of fidelity in capturing the lower tail behavior, which is critical to drought studies.
Based on the visual comparisons and goodness of fit results, it was concluded that the KDE provided excellent bivariate and univariate fits to the drought severity and corresponding standardized meteorological fluxes at all sites and over a range of climate conditions. Therefore, the KDE-based marginal and joint distributions were deemed suitable to conduct drought risk analysis and evaluate how meteorology affects agricultural droughts. Table 2. Goodness of fit statistics for the agricultural drought indicator (SSMI) at six representative stations across climate zones (see Figure 2b for the locations of stations and corresponding climate zones).

Variations of Agriculture Drought Risks with Precipitation
The joint distribution between SSMI3 and P for summer (critical growing season in the area) is shown in Figure 4. Similar distributions for SSMI3 and P for winter and annual timescales and SSMI6 and P for summer, winter, and annual periods can be found in Supplementary Information (Figures S18-S22). The distributions show similar trends across the region, but subtle differences can be seen between different stations. Stations S1, S3, and S4, which are in semi-arid and arid portions, show similar trends compared to S2, S5, and S6, which represent the more humid portions of the study area, indicating that variations in precipitation affect agricultural droughts. Figure 5 depicts a spatial variation of risks of having moderate or higher intensity agricultural droughts (SSMI ≤ −0.8) given a moderate or higher meteorological drought (Figure 4a,c,e), or when there is no meteorological drought as measured using standardized precipitation (Figure 4b,d,f) for summer, winter, and annual time series. Figure 6 shows similar variations for SSMI6, which represents full season drought. Variations of SSMI3 and SSMI6 for other cut-offs (severe, extreme, and exceptional drought) are presented in Supplementary Information  (Figures S23-S28) where additional summary statistics for these cutoffs can also be found in Tables S2  and S3. Table 3. Goodness of fit statistics for marginal distributions of meteorological fluxes at six representative stations across climate zones (see Figure 2b for the locations of stations and corresponding climate zones).      The results presented in Figures 5 and 6 indicate that in the central portions of the study area, precipitation has a relatively stronger influence on agricultural droughts, with an over 50% likelihood of observing an agricultural drought given precipitation deficits. The risks are slightly lower in the northwestern and southern portions where agricultural droughts are likely controlled by other meteorological factors. In addition, agricultural droughts in the winter given precipitation deficits exhibit a much greater variability across the study area, with risks ranging from nearly 11.6% to over 40%. This variability points to the moisture holding capacity of the soils, especially during the cooler The risk computations in Figure 5 suggest that agricultural drought risks in the central portions of the Ogallala Aquifer correspond strongly with summer and winter rainfall deficits. There is a stronger correlation between meteorological and agricultural droughts during winter than during summer in the southern portions. The northeastern portions (cold humid) of the study area also show a moderate to strong correlation between agricultural and meteorological droughts. However, the northwestern portions (cold arid) show a lesser coincidence between meteorological and agricultural droughts. Agricultural droughts tend to occur and persist in this region even when there are no meteorological droughts. The southern portions of the aquifer also exhibit a relatively stronger correlation between agricultural droughts and precipitation-related meteorological non-drought states. Figure 6, which presents similar risk profiles using SSMI6, also shows similar trends but the risk values are somewhat subdued as the accumulation period of six months dilutes the seasonal effects by aggregating non-seasonal data.

Station
The results presented in Figures 5 and 6 indicate that in the central portions of the study area, precipitation has a relatively stronger influence on agricultural droughts, with an over 50% likelihood of observing an agricultural drought given precipitation deficits. The risks are slightly lower in the northwestern and southern portions where agricultural droughts are likely controlled by other meteorological factors. In addition, agricultural droughts in the winter given precipitation deficits exhibit a much greater variability across the study area, with risks ranging from nearly 11.6% to over 40%. This variability points to the moisture holding capacity of the soils, especially during the cooler winter months. Furthermore, while the likelihood of agricultural drought occurring under wet conditions is lower than that under dry (precipitation deficit) conditions, the chances of having an agricultural drought without precipitation deficits are higher in summer than in winter, pointing to other meteorological factors influencing agricultural droughts as well.

Variation of Agriculture Drought Risks with Potential Evapotranspiration
The joint distribution of SSMI3 and PET for the summer months is shown in Figure 7 (the distributions of SSMI3-PET for other time periods and SSMI6-PET for all times can be found in Supplementary Information (Figures S29-S33). Unlike precipitation whose deficits (or negative values) indicate drought, positive values of standardized PET correspond to higher temperatures and are indicative of atmospheric deficits. Again, while trends of joint distributions are similar, variations in joint probabilities can be seen across stations. In the semi-arid to arid portions, PET influences become less strong moving northwards (stations S1, S3, and S4), which is to be expected given that the climate becomes cooler up in the north. A similar trend of less strong influence between SSMI3 and summer PET can also be seen in the humid portions of the study area (stations S2, S5, and S6) but with much lower variability. The SSMI3-PET distributions for winter and annual timescales presented in Supplementary Information (Figures S29 and S30) indicate similar trends in arid portions, but PET has a slightly higher influence on agricultural droughts in the drier northern portions during winter months. The relationships at an annual scale are largely controlled by summer PET values, which tend to be higher than winter PET values. The risks of agricultural droughts due to given unfavorable PET conditions are much less variable across the study area during the winter months, with risks ranging from 22.9% to slightly over 44%, especially in comparison to precipitation deficit conditions (see Figure 5). Figures 8 and 9 indicate that agricultural drought risks (as measured using SSMI3 and SSMI6) during summer months are also strongly correlated to PET (see Figures S34-S39 for other cut-offs). Higher values of PET generally caused by higher temperatures affect soil dryness. A comparison of probabilities from PET (Figures 8 and 9) and standardized precipitation (Figures 5 and 6) indicates that the impacts of PET on agricultural drought risks during summer are generally less compared to risks arising from precipitation deficits, but are nonetheless significant.      Figure 10 illustrates the joint distribution of SSMI3 with standardized atmospheric water deficit (AWD). Based on Equation (2), negative values indicate below normal atmospheric water (PET > P) and therefore drought conditions. The shape of the KDE distribution more closely resembles the In contrast, the soil moisture deficit risks computed using SSMI6 (presented in Figure 9) are more greatly influenced by PET than precipitation ( Figure 6). Thus, precipitation plays a greater role in controlling short-term intra-seasonal risks, but PET plays a bigger role over longer periods. This result likely arises because seasonal variations in temperature and winds (factors controlling PET) exhibit greater persistence than precipitation, which tends to exhibit greater seasonal volatility. Therefore, PET plays a much bigger role in controlling soil dryness in the southern portions of the study area, which are characterized by lower rainfalls and higher temperatures. PET has a significant influence in the central and more humid northeastern and northcentral portions of the study area during winter months. The study area is characterized by cooler winters in these regions, so even small changes in temperature can have a large effect on soil dryness.

Variations of Agriculture Drought Risks with Atmospheric Water Deficit
The right panel of Figures 8 and 9 evaluates the risk of observing agricultural droughts when the standardized PET values are close to or less than the historical means. These graphs, therefore, depict the conditions other than PET in driving agricultural droughts. The risks of droughts are much smaller when the PET values are below their statistical averages. However, these values exhibit a clear increasing trend moving from east to west, much like the decreasing precipitation gradient depicted in Figure 2c. Therefore, the results again point to the importance of precipitation in controlling the agricultural droughts within the study area. Figure 10 illustrates the joint distribution of SSMI3 with standardized atmospheric water deficit (AWD). Based on Equation (2), negative values indicate below normal atmospheric water (PET > P) and therefore drought conditions. The shape of the KDE distribution more closely resembles the bivariate joint distribution functions of SSMI3-P shown in Figure 4, more so than the SSMI3-PET shown in Figure 7, once again highlighting the relative importance of P over PET in defining agricultural droughts during summer months. Bivariate distributions for SSMI3-AWD for winter and annual timeframes and SSMI6-AWD for all three time periods are presented in Supplementary Information  (Figures S40-S44).

Variations of Agriculture Drought Risks with Atmospheric Water Deficit
Variations between bivariate distributions at different stations can be seen in Figure 10. The SSMI3 risk diminishes somewhat going northward in the arid portion of the study area (stations S1, S3, and S4) as well as in the more humid portions (stations S2, S4, and S5) albeit with lower sensitivity. This northward variation is consistent with the temperature gradient shown in Figure 2d, indicating that while the general shape of the SSMI3-AWD is influenced by precipitation, PET affects the regional variability of these distributions and therefore agricultural drought risks to net atmospheric water deficit. Figure 11 depicts the short-term agricultural drought risks given that the atmospheric moisture deficit is dry or unfavorable (AWD ≤ −0.8) or not dry (AWD > −0.8). Similar maps for other cut-offs can be found in Supplementary Information (Figures S45-S47). A comparison of AWD ( Figure 11) with P ( Figure 5) and PET ( Figure 8) indicates that AWD risk profiles are largely consistent with the patterns noted in analogous risk calculations based on standardized precipitation alone ( Figure 5). However, the magnitudes of the risks are somewhat higher (~49-~60%) when AWD is used compared to P or PET alone, especially for summer and winter seasons. This result again highlights that the regional-scale agricultural drought trends are largely controlled by precipitation, but PET variations affect the magnitude of the drought risks. Similar trends can also be seen for SSMI6 ( Figure 12); the regional trends in long-term droughts are largely conditioned by precipitation, but PET also plays a role (albeit small). However, PET represents a larger role in controlling the magnitude of the risks. Refer to Supplementary Information (Figures S48-S50) for SSMI6 maps corresponding to other cut-offs. Boxplots summarizing variability of all drought risks corresponding to different cut-offs can be found in Figures S51-S52 for SSMI3 and SSMI6.    Figures (a,c,e) and not dry states in Figures  (b,d,f)).  Figures (a,c,e) and not dry states in Figures (b,d,f)).  Figures (a,c,e) and not dry states in Figures  (b,d,f)). Figures 11 and 12 also show that agricultural risks when AWD is not under dry conditions are higher in the northwestern and southern portions, which are characterized by higher aridity during summer months; this extends over the entire northern portion of the study area during the winter months, as these areas are characterized by cold but dry winters. Furthermore, the risk of agricultural droughts given no atmospheric dryness (as measured using AWD) is generally the lowest in comparison to similar risks of agricultural droughts given above moderate precipitation or below  Figures (a,c,e) and not dry states in Figures (b,d,f)). Figures 11 and 12 also show that agricultural risks when AWD is not under dry conditions are higher in the northwestern and southern portions, which are characterized by higher aridity during summer months; this extends over the entire northern portion of the study area during the winter months, as these areas are characterized by cold but dry winters. Furthermore, the risk of agricultural droughts given no atmospheric dryness (as measured using AWD) is generally the lowest in comparison to similar risks of agricultural droughts given above moderate precipitation or below moderate PET. This result follows from the fact that the magnitude of agricultural drought risks that correlate more strongly to AWD than P or PET alone and include both of these factors therefore accounts for a higher level of agricultural risk. However, factors in addition to AWD affect the likelihood of occurrence of agricultural droughts but to a much smaller degree. Soil characteristics are implicitly accounted for using the scaling coefficients in the leaky bucket model. Intrinsic soil properties are not fully captured by the model, and the agricultural drought risks when AWD is not dry can be due to variations in soil properties across the study area.

Summary and Conclusions
The soil moisture is the master variable controlling agricultural droughts. As long-term records of soil moisture are not readily available in most areas, evaluation of multi-decadal and century-scale soil moisture drought evaluations rely on water balance models. These models use climatic fluxes (i.e., precipitation and evapotranspiration) to estimate soil moisture (referred to here as soil moisture analogs). Soil moisture analogs and agricultural drought indicators derived from them provide a convenient way to evaluate how meteorological droughts propagate through soil systems and affect agricultural droughts.
Stochastic relationships between the standardized soil moisture index (SSMI) computed at three-and six-month accumulation scales, standardized precipitation, standardized PET, and standardized AWD values for summer (April-September), winter (October-March), and annual (October-September) timescales were developed using kernel density estimation (KDE) protocols. KDE offers several advantages, as bivariate distributions can be direct fit and marginals can be integrated out. This approach, therefore, avoids the need to select different distribution functions for marginals and joint distribution and thus minimizes the subjectivity associated with the model selection process. The approach was used to fit bivariate and marginal distributions at 187 sites in the High Plains region that is underlain by the Ogallala Aquifer (the largest aquifer in the US). Excellent fits were noted based on visual inspections and use of Kling-Gupta Efficiency metric and its components, adding confidence to this model building approach.
The bivariate joint distributions and marginal distributions were used to then compute conditional probabilities and evaluate agricultural drought risks given P, PET, or AWD being in unfavorable and favorable states. The results indicate that precipitation is a major influence on agricultural droughts and the likelihood of an agricultural drought given a precipitation deficit could be as high 40-65% within the study area during summer months and between 20-55% during winter months. The relationship between agricultural drought risks and precipitation deficits is strongest in the central portions of the study area. The conditional risks of agricultural droughts given unfavorable PET conditions are higher in the eastern humid portions than the western arid portions of the study area. In general, unfavorable PET had a higher impact on SSMI6 but also seemed to influence SSMI3 to some degree. Dry states as defined by atmospheric water deficits (AWD) provided higher risks than P and PET, suggesting that both P and PET influence agricultural droughts. Agricultural drought risks under favorable conditions of AWD were much lower than when AWD was unfavorable. The agricultural drought risks were higher during the winter when AWD was favorable and point to the role of soil characteristics on agricultural droughts.
The information presented in this paper and the associated Supplementary Information provides a comprehensive risk-based drought atlas for an agriculturally important region in the US, and as such is of practical use to decision makers. The methodology developed here is also generic and can be extended to other regions with considerable ease, as global datasets required are readily available.  Figure S1: P-P Plot of SSMI3-P for Winter Months to Visually Evaluate Goodness of Fit (GOF) of KDE, Figure S2: P-P Plot of the SSMI3-P for Annual Timescale to Visually Evaluate Goodness of Fit (GOF) of KDE, Figure Figure S18: Bivariate Relationship between SSMI3-P using KDE for Winter Months, Figure S19: Bivariate Relationship between SSMI3-P using KDE for Annual Timescale, Figure S20: Bivariate Relationship between SSMI6-P using KDE for Summer Months, Figure S21: Bivariate Relationship between SSMI6-P using KDE for Winter Months, Figure S22: Bivariate Relationship between SSMI6-P using KDE for Annual timescale, Figure