Choose your preferred view mode

Please select whether you prefer to view the MDPI pages with a view tailored for mobile displays or to view the MDPI
pages in the normal scrollable desktop version. This selection will be stored into your cookies and used automatically
in next visits. You can also change the view style at any point from the main header when using the pages with your
mobile device.

Abstract

:

Severity of drought in California (U.S.) varies from year-to-year and is highly influenced by precipitation in winter months, causing billion-dollar events in single drought years. Improved understanding of the variability of drought on decadal and longer timescales is essential to support regional water resources planning and management. This paper presents a soft-computing approach to forecast the Palmer Drought Severity Index (PDSI) in California. A time-series of yearly data covering more than two centuries (1801–2014) was used for the design of ensemble projections to understand and quantify the uncertainty associated with interannual-to-interdecadal predictability. With a predictable structure elaborated by exponential smoothing, the projections indicate for the horizon 2015–2054 a weak increase of drought, followed by almost the same pace as in previous decades, presenting remarkable wavelike variations with durations of more than one year. Results were compared with a linear transfer function model approach where Pacific Decadal Oscillation and El Niño Southern Oscillation indices were both used as input time series. The forecasted pattern shows that variations attributed to such internal climate modes may not provide more reliable predictions than the one provided by purely internal variability of drought persistence cycles, as present in the PDSI time series.

1. Introduction

Drought is a fundamental feature of the climate of North America, where several regions of western United States (U.S.) have experienced protracted decadal-scale dry periods in the past centuries [1] Hydrologic droughts in western U.S. were already widespread and persistent during the so-called Medieval Climatic Anomaly, roughly in the period 900–1300 AD [2], with mega-drought in southern California during 832–1074 AD and 1122–1299 AD [3,4]. Multi-year droughts have also recurred in more recent times, e.g., in 1818–1824, 1827–1829, 1841–1848 and 1855–1865 (Figure 1), causing tremendous disruption on social, agricultural, ecological and economic fronts [5]. Five major droughts followed, which ended in 1924, 1935, 1950, 1960 and 1977. As well, the one started in 2012 [1] resulted in statewide proclamations of emergency [6]. Much of the water supply for California is derived from the Sacramento-San Joaquin River Delta (located in Northern California) via pumps located at the southern end of the delta. However, in recent times, California’s water resources have been subject to increased stress from a combination of factors including a growing population, groundwater deficit, limitations on extraction of water for the protection of fish, and increased competition for available water [6]. Over 2012–2016, drought conditions impacted surface water supplies, and increased agricultural demand and land subsidence owing to groundwater extraction. These factors inspired the development of legislation to regulate groundwater resources and financially support sustainable groundwater management as well as cleanup and storage [7]. Water management in the state (which has been studied extensively [8]) shows that the California case is exemplary of the preparedness and response measures required to cope with extreme drought events, adapt to them and build long-term resilience [9]. How drought may change in future is of great concern as global warming continues [10]. Yet, how has an extreme drought occurrence over California shifted as a result of the change in climate since historical times? How can we see droughts coming? If we are dry during one drought year, will we likely be dry for other drought years, and then for a decade or more? How cyclical will these patterns be and how are they predictable over multidecadal time-scales? To answer these questions, we examined (with focus on California) uncertainties in estimating the future ramifications of years of drought, and how drought changes may recur in the near future using the Palmer Drought Severity Index (PDSI).

Many of drought indices developed for the purpose of drought monitoring are based on meteorological and hydrological variables, which show the size, duration, severity and spatial extent of droughts. The Palmer Drought Severity Index (PDSI) is such an example. Originally developed by Palmer [11], it is one of the most well-known and widely used drought indices in the U.S. [12,13] and beyond [14,15,16,17] PDSI values are computed along the soil moisture balance that requires time series of temperature, precipitation, ground moisture content (or available water-holding capacity) and potential evapotranspiration. The calculation algorithm of PDSI—either in its original version by Palmer [11] or in modified ones [18] is thus a reflection of how much soil moisture is currently available compared to that for normal or average conditions. The PDSI incorporates both precipitation and temperature data in a simplified, though reasonably realistic, water balance model that accounts for both supply (rain or snowfall water equivalent) and demand (temperature, transformed into units of water lost through evapotranspiration), which affect the content of a two-layer soil moisture reservoir model (a runoff term is also activated when the reservoir is full). Not explicitly bounded, the PDSI typically falls in the range from −4 (extreme drought) to +4 (extremely wet). The PDSI is a dimensionless quantity for comparison across regions with radically different precipitation regimes. This means that there are limitations in the use of this index at specific scales, for which other drought indices have been developed to characterize local agricultural and socio-economic contexts [19].

Land-atmosphere interactions can introduce persistence into droughts because reduced precipitation lowers soil moisture, reduces surface evapotranspiration and, with less vapor in the atmosphere, further reduces precipitation. In this sequence, soil moisture adjustment occurs with a length of time, which introduces a lag and a memory. Depending on situations, there might be a strong coupling between soil moisture and precipitation, and land surface processes can lead to persistence [20]. The calculation of PDSI is intended to model soil moisture persistence (or memory). The combination of past wet/dry conditions with past PDSI data means that the PDSI for a given time step (generally one month) can be seen as a weighed function of current moisture conditions and a contribution of PDSI over previous times [21]. In the light of this persistence structure, PDSI chronologies can be used to reconstruct drought conditions, but persistence can also be a criterion to be used as a measure of predictability [16].

This paper deals with time series analysis (TSA) related to PDSI dynamics. Several statistical TSA approaches were applied to predict climate variables, including their extremes [22,23]. Mossad and Alazba [24] proved the potential ability of these modelling approaches to forecast drought. However, drought forecasts performed at monthly time-scale for early warning [25,26,27,28,29] do not account for long-term patterns of evolution, which are essential to study and monitor drought from a climate perspective [30]. Here, we target annual to decadal time scales. We investigate to what extent TSA model simulations may provide reliable forecasts of future hydrological changes. Although research on meteorological drought (that is, when dry weather patterns dominate) is particularly difficult because of the complex and heterogeneous character of drought processes, their temporal trends respond to climate fluctuations (e.g., large-scale atmospheric circulations). Specifically, the work explores a homogenized long series of annual PDSI data (1801–2014) as derived for California by Griffin and Anchukaitis [1] and accessible at https://www1.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/california/griffin2015drought.txt (identifier ‘precip-ONDJFMAMJ-rec-2rmse’, providing precipitation anomalies serving as reasonable proxy of PDSI data taken at the lower limit of twice the root mean square error). Then the study assesses the response of an exponential smoothing (ES) model, using an ensemble prediction approach. ES [31,32] and autoregressive integrated moving average (ARIMA) models [33] are the most representative methods in TSA. In this study, ES was used because it is known to be optimal for a broader class of state-space models than ARIMA models [32]. ES responds easily to changes in the pattern of time series [34] and is often referred to as a reference model for time-pattern propagation into the future [35,36]. It is also less complex in its formulation and, as such, it was expected to be easier in identifying the causes of unexpected results. The ensemble approach has been adopted as a way to consider uncertainty in hydrological forecasting, and thus enhance accuracy by combining forecasts made at different lead times, as in Armstrong [37] and in previous authors’ papers [38,39,40,41]. A lengthy PDSI series offers a unique opportunity to explore past interannual-to-interdecadal climate variability, under the assumption that the past interannual climate variability, with its internal dependence structure, can be used to replicate future PDSI ramifications at the local scale. This approach was compared with the more traditional TSA approach using transfer function models (TFM), introduced by Box and Jenkins [42] and re-visited by Shumway and Stoffer [43]. In this case, input time series of El Niño Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) were considered as impulse inputs to the output PDSI time series. An approximate ensemble approach was also developed under the transfer function models (TFM) framework for comparison purposes with the ES results. For California, there are now accessible accurate long-time series of PDSI. Our focus is motivated because severe and widespread drought are of particular concern for this U.S. state [44].

2. Materials and Methods

2.1. Environmental Setting and Data

The California’s climate varies widely, from hot desert to subarctic, depending on latitude, elevation, and proximity to the coast. California’s coastal regions, the Sierra Nevada foothills, and much of the Central Valley have a Mediterranean climate, with warm to hot, dry summers, and mild, moderately wet winters (Figure 2a,b). The influence of the ocean generally moderates temperature extremes, creating warmer winters and substantially cooler summers in coastal areas. The rainy period in most of the country is from November to April (Figure 2a). Prevailing westerly winds from the Pacific Ocean also bring moisture. The average annual rainfall in California is about 350 mm, with the northern parts of the state generally receiving higher rain amounts than the south. The reference evapotranspiration follows a more complex pattern, mostly in relation to elevation and distance from sea (Figure 2c). Temperature and evapotranspiration are especially important in California, where water storage and distribution systems are critically dependent on winter/spring rainfall, and excess water demand is typically met by groundwater withdrawal [45]. The PDSI time series derived from Griffin and Anchukaitis [1] reconstructed drought conditions for California.

2.2. Exponential Smoothing

The exponential smoothing (a popular scheme to produce smoothed time series) is a relatively simple prototype model for TSA-based forecasting, analysis and re-analysis of environmental variables [46,47]. It uses historical time series data under the assumption that the future will likely resemble the past, in an attempt to identify specific patterns in the data, and then project and extrapolate those patterns into the future (without using the model to identify the causes of patterns). Compared to other techniques (e.g., moving averages), which equally weight past observations, exponential smoothing apportions exponentially decreasing weights as observations get older. This means that recent observations are given relatively more weight in forecasting than older observations. To compute predictions based on the observed time series of PDSI data, we made use of available knowledge concerning the period of the system under investigation [48]. The following periodic simple exponential smoothing [35] was selected as reference model for time-pattern propagation into the future:

F(X)t+mR=α·StIt−p

(1)

where F(X)t+mR represents the m-step-ahead forecast from the annual series of the variable X (PDSI) on N years for an ensemble of R runs; St is the smoothed PDSI at decadal scale centered on time-year t (Equation (2)); α is the smoothing parameter for the data; It−p is the smoothed cycle index at the end of period t, its number being defined by the periods p in the seasonal cycle (Equation (3))

St=α·XtIt−p+(1−α)·St−1

(2)

It−p=δ·XtS(X)t+(1−δ)·It−1

(3)

where δ is smoothing parameter for cyclical indices.

2.3. Transfer Function Models

Results from seasonal exponential smoothing (that uses the temporal dependence structure of the time series itself to reproduce the time series behavior in the future) were compared to an alternative methodology based on transfer function models (TFM). It represents a linear transfer function approach where input time series potentially impacting the drought behavior at large spatial scales are used as explanatory time series variables in a lagged regression model. The methodology called TFM was introduced by Box and Jenkins [25] and re-visited by de Guenni et al. [49] and Shumway and Stoffer [43] (2017) to forecast monthly rainfall in the coast of Ecuador based on El Niño indices and model the impact of El Niño on fish recruitment, respectively.

In a TFM, the output series (in this case PDSI) can be represented as:

Y(t)=α1(B)·X1(t)+α2(B)·X2(t)+…+αk(B)·Xk(t)+η(t)

(4)

where X1(t), X2(t), …, Xk(t) are the input time series to be considered as explanatory variables contributing to the temporal dynamics of the output series Y(t) and η(t) is a stationary random process. The terms α1(B), α2(B), …, αk(B) are fractional polynomials in the back-shift operator B (such that BS(X(t) = X(t − s)) of the form:

2.4. Model Validation Methods

To ensure the optimal runs over the hold back prediction (testing validation), model parameterization was achieved by minimizing together the Root Mean Squared Error (RMSE) and the Mean Absolute Scaled Error (MASE), and maximizing the correlation coefficient (R). The commonly used RMSE quantifies the differences between predicted and observed values, and thus indicates how far the forecasts are from actual data. A few major outliers in the series can skew the RMSE statistic substantially because the effect of each deviation on the RMSE is proportional to the size of the squared error. The overall, non-dimensional measure of the accuracy of forecasts MASE [50] is less sensitive to outliers than the RMSE. The MASE is recommended for determining comparative accuracy of forecasts [51], because it examines the performance of forecasts relative to a benchmark forecast. It is calculated as the average of the absolute value of the difference between the forecast and the actual value divided by the scale determined by using a random walk model (naïve reference model on the history prior to the period of data held back for model training). MASE < 1 indicates that the forecast model is superior to a random walk. The correlation coefficient between estimates and observations [52] (anti-correlation) (perfect correlation)—assesses linear relationships, in that forecasted values may show a continuous increase or decrease as actual values increase or decrease. Its extent is not consistently related to the accuracy of the estimates. WESSA R–JAVA web [53] was used to assess model simulations with spreadsheet-based support.

In order to quantify long-range dependence and appraise the cyclical-trend patterns in the series, we estimated the Hurst [54] H exponent (rate of chaos), which is related to the fractal dimension D = 2 − H of the series. Long memory occurs when 0.5 < H < 1.0, that is, events that are far apart are correlated because correlations tend to decay very slowly. On the contrary, short-range dependence 0.0 < H < 0.5 is characterized by quickly decaying correlations, i.e., past trends tend to revert in the future (an up value is more likely followed by a down value). Calculating the Hurst exponent is not straightforward because it can only be estimated and several methods are available to estimate it, which often produce conflicting estimates [55,56]. Using SELFIS (SELF-similarity analysis [55], we referred to two methods, which are both credited to be good enough to estimate H [57]: the widely used rescaled range analysis (R/S method) [58], and the ratio variance of residuals method, which is known to be unbiased almost through all Hurst range [59]. Long-memory in the occurrence of PDSI values was also analyzed to see if the memory characteristic is correlated with the length of the time series. To determine whether this characteristic changes over time, the Hurst exponent was not only estimated for the full time series (1801–2014), but also for a shorter series starting in 1901 (the most recent period, which is also the period held out of the calibration process).

3. Results and Discussion

3.1. Data Analysis

The first step in any time-series analysis and forecasting is to plot the observations against time, to gain an insight into possible trends and/or cycles associated with the temporal evolution of datasets. Figure 3a shows that the PDSI time series presents important inter-annual and decadal variability, with smooth changes in its structure and turning points which help in orienting the choice of the most appropriate forecasting method [60]. Two homogeneity tests indicate a stepwise shift in the observational series in the years just before 1920. The Buishand range test [61] places the change point in 1969, whereas the Mann-Whitney-Pettitt test [62] locates it in 1920 but the two tests are not significant (p > 0.10), from which the series can be considered as relatively stationary.

The smoothed periodogram of the PDSI time series (Figure 4) was calculated by using the smoothing method [63] implemented in the R software [64] This estimate shows that most of the total variability in the series is associated with both short and large frequencies. The multiple observed maxima in the power spectrum confirm the complex interactions of several physical drought-triggering processes acting at several time scales. The maximum estimated spectral density occurs at frequency 0.185, which corresponds to a cycle of 5.4 years [65]. This cycle might be associated with El Niño phenomenon, but other frequencies have also an important contribution to the overall series variability.

3.2. Validation Results and PDSI Time Series Predictability

The whole of the PDSI time series (214 years of data from 1801 to 2014) was segregated into sub-sets for the purposes of training and validation (Figure 3a). The choice of 1801 as starting time of the series was driven by the necessity of having a sufficient amount of data for training without laying too long back in time, considering that with at least 50 observations are necessary for performing time-series analysis/modelling [66]. On the other hand, with at least 150–200 observations potentially reliable forecasts can be obtained for 30 to 50 steps ahead [67]. Forecasts were performed for the 40-year follow-up period (Figure 3b). Alternative initial conditions were simulated for each run, taking periods with a different start year (in 10-year steps-up from year 1801 to 1900) and periodical cycles (41, 42 and 43 years) for model training (training datasets).

For 1954–2014 (Figure 3b), the simulation results for validation testing are quite promising, judging by the closeness of ensemble prediction mean (red curve) to the observed 11-year Gaussian Filter (black curve) PDSI evolution. The results indicate that the ES model performs well at both high and low frequency variability, which is consistent with inter-annual to inter-decadal climate-variability. In fact, the residuals between predicted and observed time-series are coherent in the validation period: residual histograms and Q-Q plots do not identify substantial departures from normality in both the official run with the longest training time (Figure 5a,a1) and the average (ensemble mean) of all the runs (Figure 5b,b1).

The data are somewhat right-skewed; however, the right tail of the distribution is fairly closely approximated by the normal distribution, with some high extreme values.

In the validation stage, RMSE and MASE were equal to 1.0 and 0.68, respectively, which indicate a satisfactory performance, and that the forecast model is superior to a random walk. The estimated Hurst (H) exponent values are reported in Table 1.

With the R/S method, the H was found to be greater than 0.6 in both the whole series (0.611) and the sub-set 1901–2014 (0.743), which is around the threshold of 0.65 used by [68] to identify series than can be predicted accurately. In the case of the variance of residuals method, we have a situation in which obtained results are hard for interpretation. With an increase of the number of series terms (amount of observations), the Hurst exponent is expected to get closer to 0 [69], i.e., the memory effect decreases. However, with the variance residuals method, the estimated Hurst exponent moves away from 0.5 with the whole of the time series (0.611 against 0.550 with the sub-set 1901–2014). These apparently contradictory results can be reconciled by considering that a complex concept such PDSI is hardly captured by one metric, the Hurst exponent, which (depending on the estimation method used) may not reflect the changes of heading direction [70]. Indeed, the whole of the series (Figure 3a) shows frequent and sudden pulses of drought, with a change-point in 1917, as identified by the Buishand test [61], observed in coincidence with the early 20th century pluvial centered on 1915, which has received much attention in the western U.S. [21]. By combining these results, it can be stated that the California’s PDSI series is related with either a short-range or a long-range memory (in turn reflecting influences on the occurrence of droughts of both large-scale and small-scale climate systems), which assumes that some dependence structure exists that advocates the foreseeability of the series. We thus performed our forecasting analysis on the original time series of PDSI data.

3.3. Simulation Experiment

Once the performance of the ES model was established, the model trained over 1801–2014 periods was run to produce an ensemble of forecast paths of annual PDSI for 2015–2054. Our major interest was directed towards assessing the predictability of interdecadal variations. Several forecast members show for the coming decades (Figure 6) some trajectories following a cyclical pattern, in which PDSI may fall below and above “incipient drought”, with negligible monotonic, long-term trend. However, moving forward, ongoing changes in atmospheric circulation and associated precipitation and temperature variability in the western U.S. raise questions about the stationarity of extreme drought estimates [71].

When examining the projection of PDSI over four future decades (2015–2054), the ensemble mean value (Figure 6, black bold curve) is observed to roughly lie around the “incipient drought” class, approaching “mild drought” around 2030, although some members push to “extreme drought”. Around 2020 and 2036, PDSI forecasts approach “near normal” with some members which are inclined up to “incipient wet spell”. After the year 2040, the PDSI resumes decreasing and remains below the “incipient drought” for years.

3.4. Comparison with the Transfer Function Modelling Approach

The band of warm ocean water that develops inter-annually in the central and east-central equatorial Pacific, El Niño Southern Oscillation (ENSO), is the major source of climate variability affecting different parts of the world [72,73]. However, the Pacific Decadal Oscillation (PDO), i.e., the variation of sea surface temperatures in the Pacific Ocean north of 20° N with a warm phase and a cool phase can modulate the interannual relationship between ENSO and the global climate [74]. The teleconnection of precipitation in California with climate phases such as the PDO and ENSO are reported in literature. A warm (positive) PDO is thought to have a similar spatial precipitation signature as a positive ENSO (wet in the American Southwest and dry in the Pacific Northwest), and a cool (negative) PDO has a similar signature as a negative ENSO [75]. ENSO has an important influence on the rainfall regime in California and most of the U.S. with most dramatic impacts during the winter season [76]. The PDO is also relevant because its cool phase is linked to dry conditions in Southern California and neighboring states [77]. A plot of all available ENSO time series jointly with the PDO and PDSI time series is shown in Figure 7.

Sample cross-correlation functions (Figure 8) show that, among all ENSO indices, the ONI series produced the highest cross-correlation with the PDSI series at a lag of −1 (=0.43), with the ONI series leading by one time step (one year) the PDSI series. However, since this series is rather short (available from 1950 onwards), we selected the next highly correlated series with PDSI, i.e., el Niño3.4 index (=0.35 at lag −1), with the Niño3.4 series leading the PDSI series. Since the PDO time series is available from year 1900, this was considered the initial year for the analysis. The model training period was the interval 1900–1953 and the model validation period was the interval 1954–2014. The latter coincides with the validation period used for the ES approach (Section 3.3). An ARIMA model was fitted to the PDO series for the training period. An autoregressive model of order 1 (AR(1)) was adequate for the series. Figure 9a2 presents the sample cross-correlation function (CCF) between the PDO series (X) and the PDSI series (Y), and the sample CCF between the pre-whitened X series (residuals after fitting an ARIMA model), with the filtered Y series (after applying the AR(1) filter) presented in Figure 9a1).

Similarly, an ARIMA model was fitted to the El Niño3.4 series for the training period. An autoregressive model of order 2 (AR(2)) was adequate for the series. Figure 10b1 presents the sample cross-correlation function (CCF) between El Niño 3.4 series (X) and the PDSI series (Y), and the sample CCF between the pre-whitened X series (residuals after fitting an ARIMA model) and the filtered Y series (after applying the AR(2) filter) in Figure 10b1. Figure 9b1,b2 show a significant spike at lag = −1 for the CCF between Niño3.4, and PDSI series and PDO and PDSI series. After filtering the input and output time series to discard the autocorrelation effects, Figure 9a1,a2 show the persistent significant leading impact of the Niño3.4 and the PDO series on the PDSI series one year in advance (lag = −1). From the sample CCF functions and according to Box and Jenkins [42], a transfer function model of order (r1,s1,d1)=(1,1,1) for input series X1(t) = Niño3.4, and order (r2,s2,d2)=(1,0,1) for input series X2(t)=PDO, was proposed for this data set. In this case α1(B)=(δ01+δ11B)B1/(1−ω11B1) and α2(B)=(δ02)B1/(1−ω12B1).

The final model to be fitted is of the form:

Yt=α1Yt−1+α2Yt−2+α3Nino3.4t−1+α4Nino3.4t−2+α5PDOt−1+ηt

(5)

Following Shumway and Stoffer [43], this model was initially fitted by least squares and the ARIMA model associated to the estimated residuals η was identified. As a second step, the model was refitted assuming autocorrelated errors following an ARIMA model with order identified in the previous step. Figure 10 shows the autocorrelation and partial autocorrelation function of the estimated residuals, suggesting a white noise structure with no additional refitting required.

Figure 11 compares the observed values (black line) with the fitted values for the training period (blue line) and the observed values with the fitted values for the validation period (red line). The 95% prediction intervals are also shown in the analysis.

RMSE = 1.0 and MASE = 0.95 during the validation period indicate that the TFM provided an improvement over the naive forecast. Considering that in this case the training period (1900–1953) is much shorter in comparison with the training period used for the ES method (1800–1953), the TMS provides a competitive approach as a forecasting method for the PDSI series. The histogram and Q-Q plots of the residuals between the predicted and observed values for the PDSI time series during the validation period shows a satisfactory performance with an approximate normal distribution of the residuals (Figure 12).

3.5. Ensemble Forecast with the Transfer Function Model

Once the adequacy of the model was assessed, the model was trained over the period 1900–2014 to produce a simulation plume of annual PDSI values for the period 2015–2054. El Niño3.4 series and the PDO series were jointly simulated first, by using a multivariate ARIMA [78] model that considers dependence between the two series. The simulated values were included as external covariates for the PDSI model trained over the period 1900–2014. Simulations are presented in Figure 13 from a model of the form:

Yt=α1Yt−1+α2Yt−2+α3Nino3.4t−1+α4Nino3.4t−2+α5PDOt−1+ηt

(6)

The inter-decadal cycles observed in the ensemble forecast from the ES approach (Figure 5) are not present in this case since a seasonal component was not considered in the model. Figure 14 compares the two approaches (ES ad TFM) in the validation and forecast periods. With ES, the projections of PDSI over four future decades (2015–2054) lie around the “incipient drought” class, with episodes of “mild drought”, while the projections of the TFM remain around the “incipient drought” region.

3.6. Limitations and Perspectives

Droughts occur over long-time spans, and their timing is difficult to identify and predict. This paper takes the challenge to examine a strategy for structuring knowledge about drought dynamics, for use in annual PDSI extrapolation for the coming decades. Extrapolation suffers when a time series is subject to shocks or discontinuities. Few extrapolation methods account for discontinuities [79]. Instead, when discontinuities occur, extrapolation may lead to large forecast errors. For example, ENSO and PDO can lead to strong upward or downward trends of drought index values and frequencies [80]. According with Sheffield and Wood [81], it is plausible that thermal impacts on drought frequency in the long term are likely to dominate precipitation changes. We could thus expect a monotonic and positive temperature change with increasing drought frequency across a range of drought metrics by the late 21st century. However, the future direction of PDSI series remains uncertain, because uncertain is the direction of its causal forces (temperature and precipitation). This is a challenge in PDSI future extrapolation. Forecasts from Esfahani and Friedel [82] suggested the likelihood for the current moderate drought in California to shift to a mid-range condition in 2020 and a constant level of PDSI towards 2060. These authors advocate that California might have reached its equilibrium, the end of a long-memory process, which would be an exception in southwestern U.S., where PDSI would increase. Our empirical forecasts are in agreement with this finding.

4. Concluding Remarks

We presented a first observational assessment of Californian drought, with the focus on the PDSI at interannual time scales and the possibility offered by statistical approaches to forecast it. Our results can be useful for water resources management and planning (for instance, under the California Irrigation Management Information System programme, https://cimis.water.ca.gov), and provide basic knowledge to support further predictive studies beyond the use of Global Climate Models (GCMs) [83,84], which vary their parameters for climatic simulations under alternative GHG emission scenarios. So far, ES approaches were applied almost exclusively in econometric and financial domains [36]. The major potential cause of bias in ES is that extrapolative forecasts can differ substantially depending on the training period start date [37]. However, this bias can be dwindled by identifying long time-series and by producing a bouquet of forecasts (ensemble) from different starting points. This is what we have done with an extended series of PDSI data. This study may help in stimulating the debate about less conventional ways of understanding the climate mechanisms behind drought onset and persistence.

Author Contributions

N.D. designed and ran the exponential smoothing study, and wrote the first draft of the manuscript. L.B.d.G. ran the transfer function model. M.G. contributed to data analysis. G.B. participated in data analysis and finalized the writing of the manuscript.

Funding

The authors acknowledge that the above study is an investigators-driven research run without grant support.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Wichard, J.D.; Merkwirth, C. Robust long term forecasting of seasonal time series. In Proceedings of the 8th International Work-Conference on Artificial Neural Networks, Barcelona, Spain, 8–10 June 2005. [Google Scholar]

Figure 4.
Smoothed periodogram of the PDSI time series (bandwidth is a measure of the width of the frequency interval used in the smoothing procedure).

Figure 4.
Smoothed periodogram of the PDSI time series (bandwidth is a measure of the width of the frequency interval used in the smoothing procedure).

Figure 5.
(a,a1) Residual histogram and normal Q-Q plot between PDSI forecasted and observed in validation period for the official run; (b,b1) Residual histogram and normal Q-Q plot between PDSI forecasted and observed in validation period for the ensemble mean.

Figure 5.
(a,a1) Residual histogram and normal Q-Q plot between PDSI forecasted and observed in validation period for the official run; (b,b1) Residual histogram and normal Q-Q plot between PDSI forecasted and observed in validation period for the ensemble mean.

Figure 8.
Sample cross-correlation functions between the PDSI series and all ENSO indices shown in Figure 7. Numbers indicate the associated lag at the peak value.

Figure 8.
Sample cross-correlation functions between the PDSI series and all ENSO indices shown in Figure 7. Numbers indicate the associated lag at the peak value.

Figure 9.
(a1) Sample cross-correlation function (CCF) between the pre-whitened Pacific Decadal Oscillation (PDO) series denoted as X and the filtered PDSI series denoted as Y for the training period; (a2) the sample CCF between the PDO series and the PDSI series; (b1) CCF pre-whitened Niño3.4 and the filtered PDSI; (b2) the sample CCF between the Niño3.4 series and the PDSI series.

Figure 9.
(a1) Sample cross-correlation function (CCF) between the pre-whitened Pacific Decadal Oscillation (PDO) series denoted as X and the filtered PDSI series denoted as Y for the training period; (a2) the sample CCF between the PDO series and the PDSI series; (b1) CCF pre-whitened Niño3.4 and the filtered PDSI; (b2) the sample CCF between the Niño3.4 series and the PDSI series.

Figure 11.
Observed PDSI time series (black line) with the training dataset to build the model (blue line) for the period 1900–1953, including the filtered observed values (navy blue). Also comparison between the observed values (black line) and predicted values with the TF model for the validation period (1954–2014) (red line), including the corresponding 95% confidence intervals (red dash).

Figure 11.
Observed PDSI time series (black line) with the training dataset to build the model (blue line) for the period 1900–1953, including the filtered observed values (navy blue). Also comparison between the observed values (black line) and predicted values with the TF model for the validation period (1954–2014) (red line), including the corresponding 95% confidence intervals (red dash).

Figure 12.
Residual histogram and normal Q-Q plot between the PDSI predicted values and observed PDSI time series for the validation period (1954–2014).

Figure 12.
Residual histogram and normal Q-Q plot between the PDSI predicted values and observed PDSI time series for the validation period (1954–2014).

Figure 13.
Observed PDSI time series (black line) with the simulation plume (grey lines) for the period 2015–2054, including the filtered observed series (navy-blue line), the median of the simulated values (thick grey line) and the 2.5% (bottom dashed line) and 97.5% quantile (top dashed line) of the simulated values.

Figure 13.
Observed PDSI time series (black line) with the simulation plume (grey lines) for the period 2015–2054, including the filtered observed series (navy-blue line), the median of the simulated values (thick grey line) and the 2.5% (bottom dashed line) and 97.5% quantile (top dashed line) of the simulated values.

Figure 14.
Comparison between estimates from the exponential smoothing model (ESM) and the Transfer Function model (TFM) for both the validation period (1954–2014) and the forecast period (2015–2054).

Figure 14.
Comparison between estimates from the exponential smoothing model (ESM) and the Transfer Function model (TFM) for both the validation period (1954–2014) and the forecast period (2015–2054).

Table 1.
Estimated values of the Hurst (H) exponent (with two methods) for the PDSI annual series as a whole and for a reduced number of years.

Table 1.
Estimated values of the Hurst (H) exponent (with two methods) for the PDSI annual series as a whole and for a reduced number of years.