A pseudoproxy assessment of data assimilation for reconstructing the atmosphere–ocean dynamics of hydroclimate extremes

Abstract

Because of the relatively brief observational record, the climate dynamics that drive multiyear to centennial hydroclimate variability are not adequately characterized and understood. Paleoclimate reconstructions based on data assimilation (DA) optimally fuse paleoclimate proxies with the dynamical constraints of climate models, thus providing a coherent dynamical picture of the past. DA is therefore an important new tool for elucidating the mechanisms of hydroclimate variability over the last several millennia. But DA has so far remained untested for global hydroclimate reconstructions. Here we explore whether or not DA can be used to skillfully reconstruct global hydroclimate variability along with the driving climate dynamics. Through a set of idealized pseudoproxy experiments, we find that an established DA reconstruction approach can in principle be used to reconstruct hydroclimate at both annual and seasonal timescales. We find that the skill of such reconstructions is generally highest near the proxy sites. This set of reconstruction experiments is specifically designed to estimate a realistic upper bound for the skill of this DA approach. Importantly, this experimental framework allows us to see where and for what variables the reconstruction approach may never achieve high skill. In particular for tree rings, we find that hydroclimate reconstructions depend critically on moisture-sensitive trees, while temperature reconstructions depend critically on temperature-sensitive trees. Real-world DA-based reconstructions will therefore likely require a spatial mixture of temperature- and moisture-sensitive trees to reconstruct both temperature and hydroclimate variables. Additionally, we illustrate how DA can be used to elucidate the dynamical mechanisms of drought with two examples: tropical drivers of multiyear droughts in the North American Southwest and in equatorial East Africa. This work thus provides a foundation for future DA-based hydroclimate reconstructions using real-proxy networks while also highlighting the utility of this important tool for hydroclimate research.

Introduction

Hydroclimate extremes, including persistent droughts and pluvials, can have extensive impacts on human welfare and agricultural production. While the frequency of such extremes is estimated to increase with global warming e.g.,, it is unclear if the underlying climate dynamics of such events are accurately produced in climate model simulations e.g.,. The approximately 100-year observational record also limits our ability to understand the mechanisms of these events because this time interval does not capture the full range of hydroclimate extreme events and it poorly characterizes multidecadal to centennial variability. A canonical example is the multidecadal “megadroughts” that were a prominent though infrequent feature of the climate system over the past 2000 years in western North America and Europe . Of particular concern is that many of these drought events were more persistent than any of the droughts observed over the historical era . Reconstructing hydroclimate variables along with their driving dynamics would therefore provide a critical perspective and allow for improved characterizations of such high-impact events.

While hydroclimate reconstructions exist over many different timescales and resolutions (e.g., ; ), hydroclimate reconstructions over the Common Era (the past two thousand years) are particularly important because this time period allows for seasonal- and annual-timescale reconstructions over much of the globe . Hydroclimate reconstructions for the Common Era have been performed for localized regions, such as for snowpack and precipitation in specific mountain ranges e.g., or for the streamflow of particular rivers e.g.,. On continental scales, gridded “drought atlases” targeting the Palmer drought severity index (PDSI) have been derived over portions of North America , Europe , Southeastern Asia , and Australia and New Zealand . Despite these regional successes, no truly global hydroclimate reconstructions have been attempted.

The goal of this work is to demonstrate and test the application of data assimilation (DA) for global hydroclimate reconstructions. DA is a method that optimally fuses proxy information with the dynamical constraints of climate models . We lay this groundwork by performing a set of synthetic “pseudoproxy” reconstructions that estimate a realistic upper bound on reconstruction skill: given the current proxy networks and state-of-the-art climate models, we ask whether or not it is possible to skillfully reconstruct hydroclimate variables and their related dynamics. Pseudoproxy experiments provide a test of the reconstruction methodology while controlling for uncontrolled and sometimes unknown factors that pervade real-proxy reconstructions seefor a review. Pseudoproxy experiments have been successfully used to test the skill and characteristics of several reconstruction methodologies, including DA, and reconstructed variables such as temperature and precipitation e.g.,. The present pseudoproxy reconstruction work will help to focus future DA reconstruction efforts on the variables, timescales, and geographic areas that are most likely to be successfully reconstructed. This paper also showcases the new kinds of analyses that are possible with the dynamical atmosphere–ocean variables from a DA-based climate reconstruction. We first analyze the pseudoproxy reconstructions by focusing on global-scale temperature and hydroclimate indices. We then move into specific analyses of the mechanisms of drought in the North American Southwest and equatorial East Africa.

Experimental framework

DA methodology

We employ a DA technique that optimally combines observations (in this context, proxy data) with climate model states. The model provides an initial, or prior, state estimate that is updated in a Bayesian sense based on the observations and an estimate of the errors in both the observations and the prior. The prior may contain any climate model variables of interest. The updated prior, called the posterior, is the estimate of the climate state given the observations and the error estimates. The basic state update equations of DA e.g.,Chap. 5 are given by $${\mathit{x}}_{\text{a}}={\mathit{x}}_{\mathrm{b}}+\mathbf{K}[\mathit{y}-\mathcal{H}({\mathit{x}}_{\text{b}}\left)\right],$$ where ${\mathit{x}}_{\text{b}}$ is the prior (or “background”) estimate of the state vector and ${\mathit{x}}_{\text{a}}$ is the posterior (or “analysis”) state vector. Observations (or proxies) are contained in vector $\mathit{y}$. The observations are estimated by the prior through $\mathcal{H}\left({\mathit{x}}_{\text{b}}\right)$, which is in general a nonlinear vector-valued observation operator that maps ${\mathit{x}}_{\text{b}}$ from the state space to the observation space. The Kalman gain matrix $\mathbf{K}$ can be written as $$\mathbf{K}={\mathbf{BH}}^{\text{T}}[{\mathbf{HBH}}^{\text{T}}+\mathbf{R}{]}^{-\mathrm{1}},$$ where $\mathbf{B}$ is the prior covariance matrix, $\mathbf{R}$ is the error covariance matrix for the proxy data, and $\mathbf{H}$ is the linearization of $\mathcal{H}$ about the mean value of the prior. Specifically, we use an ensemble square root filter to compute Eq. (). We also note that Eq. () assumes that $\mathit{y}$, ${\mathit{x}}_{\mathrm{b}}$, and $\mathcal{H}\left({\mathit{x}}_{\text{b}}\right)$ are Gaussian distributed and that their errors are unbiased. The reconstruction process works by essentially computing an optimal linear fit between the initial guess of the climate state (the prior ${\mathit{x}}_{\mathrm{b}}$) and the proxies ($\mathit{y}$). This process involves iteratively computing Eq. () for each year of the reconstruction to arrive at ensemble estimates (the posterior ${\mathit{x}}_{\mathrm{a}}$) of the climate for each year. In any given year of the reconstruction, ${\mathit{x}}_{\mathrm{a}}$ is the probability distribution of states that are consistent with the proxy observations and errors and the prior distribution; therefore it is the full probability distribution of ${\mathit{x}}_{\mathrm{a}}$ that represents the probabilistic reconstruction of the climate state. In the time series figures herein, we show the mean of ${\mathit{x}}_{\mathrm{a}}$, with uncertainty estimates derived from the posterior distribution. For more mathematical details of the reconstruction methodology and the precise calculation procedure, see the Appendix of .

As in previous studies we use an “offline” DA approach in which the prior distribution is drawn from existing climate model simulations. For this approach, the ensemble members are seasonally or annually averaged climate states instead of an ensemble of independently running model simulations, as in “online” DA. In principle, the ensemble members can be drawn from a single long simulation, multiple simulations, or even from simulations of a collection of climate models; the only important requirement is that the prior be climatologically representative of what one is trying to reconstruct (e.g., to reconstruct a year with a large volcanic eruption, the prior should contain ensemble members that come from simulation years with large volcanic eruptions). Because of how the prior is constructed here, it does not contain year-specific forcings or boundary conditions, information which appears to be superfluous according to many previous reconstructions experiments . The reconstruction process is also performed for each year independently such that no information is propagated forward in time. As discussed in and , the offline approach is well suited to a paleoclimate context in which it has been shown to be highly skillful at annual timescales without the immense computational costs of a traditional online approach. Moreover, tests of offline vs. online approaches for paleoclimate have so far shown no improvement in reconstruction skill with an online method .

Climate model simulations

Here we explore hydroclimate reconstructions using DA with a series of pseudoproxy experiments. For these experiments, we employ two full-forcing simulations of the Community Earth System Model from the Last Millennium Ensemble Project (CESM LME; ). These CESM simulations used $\sim$ 2${}^{\circ }$ atmosphere and land components and $\sim$ 1${}^{\circ }$ ocean and sea ice components. The simulations were run from the years 850 to 1850 CE using estimates of the transient evolution of solar intensity, volcanic emissions, greenhouse gases, aerosols, land-use conditions, and orbital parameters. The simulations were given identical forcings but differed by round-off error in the initial atmospheric state; this difference was sufficient to generate simulations with different time histories (e.g., Nino 3.4 indices from the two simulations we use here have an annual average correlation of 0.097). These simulations provide the climate states for our prior ensembles and the climate inputs for the pseudoproxies: the pseudoproxies, $\mathit{y}$, are generated for the CESM LME simulation 10, while the prior, ${\mathit{x}}_{\mathrm{b}}$, and the prior estimate of the pseudoproxies, $\mathcal{H}\left({\mathit{x}}_{\mathrm{b}}\right)$, are from simulation 9. In other words, simulation 10 plays the role of the actual climate in a real-world reconstruction setting; the pseudoproxies are all derived from that simulation, while simulation 9 generates the ensemble of states that define the prior. We use two different simulations in this way so that there is no overlap between the specific prior states and the reconstructed climate states.

Pseudoproxy construction

We generate realistic pseudoproxies via process (or forward) models for each proxy type seefor a review, which transform the simulated climate signal (e.g., temperature, precipitation) into synthetic proxy observations. For the pseudoproxy experiments, we construct pseudo-tree-ring width and an approximation of pseudo-coral ${\mathit{\delta }}^{\mathrm{18}}$O. We only consider these two proxy types because they have established proxy system models and because they represent the majority of proxies in global databases (e.g., 92 % of proxies in the PAGES2k network; ). Because our experimental construction defines an upper limit on reconstruction skill, we do not consider age uncertainties, which would act to degrade reconstruction skill and are not applicable for the tree-ring proxies that constitute the majority of the proxy network. We use the pseudoproxy network shown in Fig. , which is based on the updated PAGES2k network of real, global proxy data . In this pseudoproxy network, we created only one pseudoproxy per climate model grid cell, though in real-proxy reconstructions there would be the possibility of having multiple proxies within a given grid cell. Reconstructions were performed using both pseudo-tree and pseudo-coral proxies together and also using pseudo-trees alone.

We generated 160 pseudo-tree-ring widths using the model VS-lite , which uses monthly temperature and precipitation from the CESM LME simulations to calculate tree-ring width time series at each pseudo-tree-ring location. This model accounts for the seasonal dependence of tree growth via simple parameterizations of temperature and moisture threshold responses at monthly scales and latitudinal light availability. We estimated the four growth parameters of VS-lite for each pseudo-tree location using the real PAGES2k tree-ring width data , the monthly observational temperature and precipitation data of CRU TS3.23 , and the VS-lite parameter estimation code from . Note also that because this proxy system model takes care of the averaging timescale of each individual tree, there is no need to define a rigid timescale imposed on all tree sites (e.g., all trees in the Northern Hemisphere reflect only June–July–August temperature).

We generated 55 pseudo-coral ${\mathit{\delta }}^{\mathrm{18}}$O values through the linear proxy system model of , which for a given coral proxy can be written as $${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\text{pseudo-coral}}=a_{\mathrm{1}}\text{SST}+a_{\mathrm{2}}\text{SSS}+\mathit{ϵ},$$ where SST is local, annually averaged sea surface temperature and SSS is local, annually averaged sea surface salinity. We estimated the parameters $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$ based on the real coral proxies in the PAGES2k database and annual (defined as April to the next March) SST and SSS from the ocean reanalysis dataset SODA . Note that the CESM simulations are not isotope enabled and the coral proxy system model only requires SST and SSS to create coral ${\mathit{\delta }}^{\mathrm{18}}$O pseudoproxies.

Within the formulations for both the pseudo-tree-ring widths and pseudo-coral ${\mathit{\delta }}^{\mathrm{18}}$O, each have realistic parameters estimated with real-proxy and historical data but are driven offline by climate time series inputs from the CESM LME simulations, viz., 2 m air temperature, total precipitation, SST, and SSS. While the combination of the climate simulation output and the proxy system models produce single proxy values for a given year, these values represent the many short timescale processes included in the climate model simulations, such as daily weather events.

Pseudoproxy network used in this study based on the PAGES2k global database of tree-ring and coral records , and the drought regions explored in this study, with the North American Southwest in blue and equatorial East Africa in purple.

[Figure omitted. See PDF]

We additionally note a special issue regarding the use of the VS-lite proxy system model within our Kalman-filter-based DA method: we find the output of VS-lite for some locations to be strictly non-Gaussian. This is important because Eq. () assumes that $\mathit{y}$ and $\mathcal{H}\left({\mathit{x}}_{\mathrm{b}}\right)$ are Gaussian. In preliminary work, we observed “filter divergence” (in which the state estimate diverges unrealistically away from the true state) for small prior ensembles when using VS-lite pseudoproxies. We have not observed this effect with the proxy system models for corals or ice cores . A common though ad hoc way of dealing with filter divergence is to employ covariance localization that smooths out spurious long-range correlations, acts to boost the “effective” ensemble size, and tends to reduce non-Gaussian issues . Using covariance localization is unattractive, however, because it is ad hoc and ensemble sizes can be very large in an offline DA approach. We therefore found that to avoid filter divergence we needed over 500 ensemble members in our prior. We did not find much sensitivity to the results with ensemble sizes between about 600 and 1000 members (e.g., maximum changes of 0.05 in the correlation skill of the global mean temperature reconstructions). We therefore performed the reconstructions with 1000 members, the maximum available years for the CESM LME simulations. This does not limit future real-proxy reconstructions because of the many millennial-length climate model simulations available , but it did require that we use the CESM LME in the present pseudoproxy experiments because it is the only high-resolution, CMIP5/PMIP3-class climate model with multiple millennial-length simulations that are publicly available.

A further idealization in these reconstructions is that we consider only measurement error. Specifically, we add white noise with a standard deviation of 0.01 (approximate width measurement uncertainty) to the pseudo-tree-ring time series that already have a variance equal to 1 and 0.1 (approximate isotope measurement uncertainty) to the pseudo-coral time series that are in units of per mil; the diagonal elements of $\mathbf{R}$ are ${\mathrm{0.01}}^{\mathrm{2}}$ for pseudo-trees and ${\mathrm{0.1}}^{\mathrm{2}}$ for pseudo-corals. This choice is motivated by our interest in the best-case scenario so as to estimate an upper bound for DA-based reconstructions. Our choice is also motivated by the fact that quantifying errors in nonlinear proxy system models like VS-lite is nontrivial , and thus we can have more confidence in an upper-bound reconstruction skill estimate than a “most-realistic” skill estimate. However, in the context of traditional regression-based reconstruction approaches, multivariate signals in proxies are “noise.” For reference, when the no-noise-added tree-ring and coral pseudoproxies are compared against the climate model PDSI and SST, respectively, the tree rings have a mean signal-to-noise ratio (SNR) of 0.46 with a standard deviation of 0.50 and the corals have a mean SNR of 0.41 with a standard deviation of 0.23 (calculated following , $\text{SNR}=\left|\mathit{\rho }\right|/\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}$, where $\mathit{\rho }$ is the correlation between the proxy and the nearest grid point climate variable time series). These SNRs are therefore typical of estimates for real proxies and are comparable to the SNRs adopted in pseudoproxy contexts that test traditional regression approaches .

Reconstruction procedure and experiments

Here we employ a “perfect model” framework in which the same underlying climate model is used for both the prior and the “true” climate model state. The prior, ${\mathit{x}}_{\mathrm{b}}$, consists of all the climate states at their relevant resolutions (including monthly, seasonal, and annual resolutions as defined below) from the CESM LME simulation number 9. Using this prior, we reconstruct the corresponding annual and seasonal states of simulation 10 based only on our prior and the pseudoproxies that are generated from simulation 10. Despite the idealization of the perfect model framework, previous DA-based reconstruction work has shown similar results using different climate models in both pseudoproxy and real-proxy reconstructions that did not assume a perfect model framework . These studies saw similar basin and continent-wide patterns and magnitudes of reconstruction skill.

For the reconstructions, the annual states are constructed out of monthly CESM LME output averaged from April to the next calendar year's March. This annual average is chosen to account for the seasonal cycle of a global network of proxies and climate phenomena like the El Niño–Southern Oscillation (ENSO), the continuity of which would be ignored with a calendar year average. The seasonal states we use are June–July–August (JJA) and December–January–February (DJF) averages of the monthly climate variables.

We include the following global variables in our state vectors: 2 m air temperature, the Palmer drought severity index (PDSI), and the standardized precipitation evapotranspiration index (SPEI) using a 12-month decaying exponential weighting kernel chosen to closely resemble the timescale of PDSI. Both PDSI and SPEI were computed using the Penman–Monteith equation for potential evapotranspiration and monthly climate model output of precipitation, 2 m temperature, vapor pressure, net surface radiation, surface pressure, and surface wind (estimated from 10 m down to 2 m using the wind profile power law). Here we focus on PDSI and SPEI because they are so widely used for current drought monitoring and historical drought reconstructions. In addition to these global variables, we also include the following index variables: the location of the intertropical convergence zone (ITCZ) defined by the precipitation centroid (center of mass) of the zonal mean tropical precipitation (20${}^{\circ }$ N to 20${}^{\circ }$ S; ) and the monthly Nino 3.4 index. Note that both of these variables are computed upfront and included as index variables in the prior state vector rather than being post-processed from reconstructed spatial climate fields (both approaches give equivalent results); this was done simply to save computational memory and speed up the reconstructions.

As in previous work , we quantify uncertainty in our reconstructions through a combination of Monte Carlo sampling of the proxy network and accounting for the spread in the posterior distribution, ${\mathit{x}}_{\mathrm{a}}$. Each 1000-year reconstruction is repeated 25 times in a Monte Carlo fashion by sampling 75 % of the pseudoproxy network for each reconstruction iteration (we only perform 25 iterations because the small pseudoproxy errors lead to very similar reconstructions). This gives us 25 equally likely reconstructions that each come with a 1000-member posterior ensemble estimate for each year of the reconstruction. Because this generates very large data files, we save the ensemble mean for each year of the Monte Carlo iterations and the full posterior ensembles for only one of the iterations. This yields two distributions for each year based on the proxy network sampling and the spread in climate model states. We compute the standard deviation of each of these distributions and add them in quadrature, giving us the final uncertainty estimates for all of the variables in space and time. For this particular experimental setup, we note that it is actually only necessary to save one full posterior ensemble because the posterior spread depends only on the spread in the prior, the spread in the prior estimate of the proxies, the proxy error estimates, and the number of proxies (see, for example, Eqs. A2–A4 in ), each of which are identical for all the Monte Carlo iterations. The upshot is that the ensemble spread (but not the ensemble mean) is identical for each year and each reconstruction iteration.

All of the reconstruction experiments here have been designed to show upper-bound skill estimates. This experimental framework allows us to highlight what is theoretically possible while also making falsifiable claims about what is impossible to reconstruct. To summarize, experimental idealizations discussed previously include the following: a “perfect model” setup in which the covariance structures of $\mathbf{B}$ are derived from the same climate model (but not the same simulation) as the truth and $\mathcal{H}$ is known perfectly; and low experimental error, where $\mathbf{R}$ represents only measurement error, is based on white noise and is assumed to be diagonal (proxy errors are uncorrelated). In a real reconstruction, the climate model covariances will not be identical to the covariance structures of Earth's climate, the proxy system model $\mathcal{H}$ will be imperfectly known, and proxies will include “errors” that do not reflect climate variability and are larger than measurement errors . Using this DA approach, previous studies have looked at the effects of each of these idealizations within a pseudoproxy framework . These studies show that the uncertainties related to the imperfect climate models, proxy system models, and proxies will likely act to degrade reconstruction skill but not to the point at which DA-based reconstructions provide no meaningful information. However, because such experiments are synthetically constructed, there is significant uncertainty about their fidelity to the real climate reconstruction problem; for example, used a model of intermediate complexity to reconstruct a simulation from a state-of-the-art coupled model to test the effect of inexact model covariance structures of the prior, but it is unclear how representative this comparison is to a real-world reconstruction. Working within a more idealized context circumvents such issues but also affects the framing and interpretation of the results. These upper-bound skill estimates indicate what is theoretically possible and are most useful in showing which regions and variables lack skill and are thus likely unreconstructable.

Global reconstruction skill assessment

We first assess the spatial temperature reconstructions in Fig. , which shows the skill of the reconstruction at each grid point using the metrics of correlation ($r$) and the mean continuous ranked probability skill score (CRPSS). Seasonal (JJA and DJF) and annual reconstructions are organized by column. Correlation is computed using only the reconstruction mean time series at each grid point, while the CRPSS metric accounts for both the mean grid point time series and the grid point uncertainty estimates. CRPSS is based on the continuous ranked probability score (CRPS), which is a “strictly proper” scoring rule that accounts for the skill of the entire posterior reconstruction distribution ; because the posterior ensemble estimates are normally distributed, we use Eq. (5) from , which is given by $$\text{crps}=\mathit{\sigma }\mathit{\left\{}{y}_n\right[\mathrm{2}\mathrm{\Phi }\left(y_n\right)-\mathrm{1}]+\mathrm{2}\mathit{\varphi }(y_n)-{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathit{\pi }}}}\mathit{\}},$$ where $y_n=(y-\mathit{\mu })/\mathit{\sigma }$, with $y$ being the true value, $\mathit{\mu }$ is the mean of the posterior ensemble estimate, $\mathit{\sigma }$ is the standard deviation of the posterior ensemble, and $\mathit{\varphi }\left(y_n\right)$ and $\mathrm{\Phi }\left(y_n\right)$ are respectively the normal probability density function and the normal cumulative distribution function of $y_n$. All of our uses of Eq. () are for time series, either individual or grid point. We therefore compute the mean of all the time-step values and denote it as CRPS. The skill score version, CRPSS, is the reconstructed CRPS computed with respect to the CRPS of a reference distribution, $\text{CRPSS}\equiv \mathrm{1}-{\text{CRPS}}_{\text{rec}}/{\text{CRPS}}_{\text{ref}}$, here the initial uninformed prior; positive CRPSS indicates that the reconstructed distribution is more skillful for this metric than the uninformed prior. CRPSS is a much more stringent skill metric than correlation, so we focus most of our attention on CRPSS.

Returning again to Fig.  and focusing on the CRPSS results, we note that skill is highest in the densest proxy regions. Comparing the first and second columns in Fig. , skill is also much higher near the pseudo-trees during their primary growing season (note that the full growing season for each pseudo-tree site is not preset but determined by VS-lite based on the local temperature and moisture conditions that allow or suppress growth). Despite these seasonal dependencies, the annual mean temperature reconstructions in the third column of Fig.  are skillful (positive CRPSS) despite the fact that the pseudo-trees are not responsive to the full annual cycle of climate. This is possible because the DA algorithm can exploit climate field covariance information across timescales; see and note that, for example, ENSO and the approximately 500-day thermal timescale of the global ocean mixed layer are physical features of the climate system that can retain information across seasonal and annual timescales. The skill of the temperature reconstructions in Fig.  is high across the tropics and is not solely dependent on the presence of corals there, though they do improve tropical reconstructions (see Sect. ).

Reconstruction skill for 2 m air temperature. Rows show the skill metrics of correlation ($r$) and the mean continuous ranked probability skill score (CRPSS). These skill metrics are computed for the entire 1000-year reconstruction against the true grid point time series, CESM LME simulation 10. Columns show the reconstruction skill for June–July–August (JJA) and December–January–February (DJF) seasonal means and the annual mean from April to the next calendar year's March. Pseudo-tree sites are indicated by blue squares, and pseudo-corals are indicated by blue circles.

[Figure omitted. See PDF]

One more subtle feature of Fig.  is that not all proxy locations show similar local temperature reconstruction skill, particularly for the pseudo-tree rings. This suggests that certain pseudoproxy locations provide more local information about temperature than other locations. Importantly, the tree-ring proxy system model VS-lite explicitly models how much of a given site's growth is affected by either temperature or moisture limitations (the two growth-controlling factors modeled by VS-lite). Following , we used VS-lite's growth responses (optional outputs of VS-lite) to calculate which locations were either moisture or temperature limited, expressed here by a “moisture limitation fraction,” as shown in Fig. . Specifically, we computed the fraction of summer months over the entire simulation in which the growth response to soil moisture was less than the growth response to temperature; summer is defined as JJA in the Northern Hemisphere and DJF in the Southern Hemisphere. Values of this fraction approaching zero indicate a temperature-limited site, while values approaching one indicate a moisture-limited site.

Comparison of moisture limitation fractions with the reconstruction skill estimates in Fig.  shows that the locations with the highest CRPSS skill values correspond to locations that are most temperature limited. For instance, the four tree-ring sites in Arctic Russia reveal a range of temperature- or moisture-limited responses. The CRPSS regional values in JJA follow these sensitivities, namely the most temperature-limited proxy generates the highest regional skill, while the easternmost proxy that is moisture limited is associated with the lowest regional skill. This generally holds for other isolated locations or clusters of similarly limited trees, such as for the cluster of temperature-limited sites in Fennoscandia or the cluster of moisture-limited sites in the Pamir mountain range. Note that the identification of certain moisture- or temperature-limited sites is primarily due to the CESM model simulation output, so the identifications here may not correspond perfectly to the related real-world proxy sites.

Moisture limitation fraction for each pseudo-tree ring, with green values indicating more moisture-limited sites and red values indicating more temperature-limited sites. A value of 0 indicates a pseudo-tree site that is completely temperature limited, while a value of 1 indicates a completely moisture-limited site.

[Figure omitted. See PDF]

Figures and show the skill of the reconstructed PDSI and SPEI. Both indices behave very similarly, with slightly higher skill values for SPEI (e.g., over boreal Canada and Russia for JJA correlation and CRPSS). Similar to the temperature reconstruction, the higher skill values tend to be local to the proxy sites, especially for the CRPSS metric. This result particularly holds for the areas of single moisture-limited proxies or clusters thereof. These dependencies are again illustrated by comparing Figs. , , and around the four Arctic Russia sites where skill is increased around the moisture-limited locations, particularly for the far eastern proxy in that region. Similarly, the cluster of moisture-limited sites in the Pamir mountain range is associated with increased skill there and westward to the Caspian Sea. In contrast, the arc of Himalayan pseudo-tree rings is largely temperature limited, and despite their density they have low CRPSS skill values for the PDSI and SPEI reconstructions. We also note that the skill over Australia and the equatorial regions is not dependent on the presence of the pseudo-corals (not shown), which is consistent with some of the findings of real-proxy drought reconstructions that indicate the long covariance length scales in these regions. Similar to the temperature reconstructions, Figs.  and show that skillful annual reconstructions of hydroclimate are in principle possible in this DA framework despite the fact that the pseudo-trees are primarily responding to growing season conditions. This result may be due in part to the persistence built into the drought indices.

Reconstruction skill, correlation ($r$), and the mean continuous ranked probability skill score (CRPSS) for the Palmer drought severity index (PDSI); cf. Fig. .

[Figure omitted. See PDF]

Because the reconstructions shown here provide an upper-bound skill estimate, we expect Figs. , , and to be most useful in showing where and why the reconstructions lack skill. If skill is limited under idealized conditions, it is unlikely that increased skill would be achieved in real reconstructions. In this case CRPSS is more useful because correlation is less stringent in that it only tracks the phase of the mean reconstruction and results in generally quite high correlations for most locations and seasons. The comparisons we have so far made between Fig.  and the skill metrics in Figs. , , and illustrate a general rule that without moisture-sensitive trees skillful PDSI or SPEI reconstructions are unlikely; similarly, without temperature-sensitive trees skillful temperature reconstructions are unlikely. This might at first appear to be an obvious observation. It is nevertheless important to recall that PDSI and SPEI are both dependent on temperature through their incorporation of potential evapotranspiration. Therefore it might be possible in principle for tree rings to reconstruct both temperature and moisture indices equally well regardless of whether the trees are moisture or temperature limited. To the contrary, our experiments suggest that this is not the case. Based on these experimental results, we recommend that future DA-based hydroclimate reconstructions use a spatial mixture of temperature- and moisture-sensitive trees if the goal is to reconstruct both temperature and hydroclimate indices (as happens to be the case over western North America in our experiments). Where this is impossible, our pseudoproxy experiments suggest that reconstruction skill will be significantly diminished for the target variable not represented by the proxy network.

Regional hydroclimate reconstructions

Given the above demonstrations of global reconstruction skill, we now showcase two regional examples that explore the driving dynamics of multiyear droughts in the North American Southwest and equatorial East Africa within the pseudoproxy reconstructions.

Reconstruction skill, correlation ($r$), and the mean continuous ranked probability skill score (CRPSS) for the standardized precipitation evapotranspiration index (SPEI) with a 12-month decaying exponential weighting kernel; cf. Figs.  and .

[Figure omitted. See PDF]

North American Southwest area-averaged JJA PDSI in the reconstruction mean (purple) and the climate model truth (black). Both time series have been smoothed with an 11-year moving average . The skill metrics of correlation ($r$) and the mean continuous ranked probability skill score (CRPSS) for the smoothed time series are shown in the lower left corner. The five most severe multiyear droughts over this period are highlighted in corresponding colors. Reconstruction error estimates, similar to those shown in Fig. , have been left off to clearly show the severe drought highlighting.

[Figure omitted. See PDF]

North American Southwest multiyear drought

We first focus on the dynamical mechanisms of multiyear droughts in the North American Southwest. The choice of this region is motivated by the local reconstruction skill seen in Figs.  and and because several established theories exist about the dynamical drivers of drought in this region e.g.,. Figure  shows the time series reconstruction of JJA PDSI in the North American Southwest region defined by the land area within the latitudinal range of 25–42.5${}^{\circ }$ N and the longitudinal range of 125–105${}^{\circ }$ W as previously used in and indicated by the blue box in Fig. . In Fig. , the five most severe droughts are highlighted for both the reconstruction and the truth; overlap is indicated by the slightly darker purple highlights. The cumulative drought severity is ranked by the sum of persistently negative PDSI values relative to an 11-year moving average . Very similar results were found using an alternative definition of drought that commences after two consecutive years of negative PDSI values and continues until two consecutive years of positive PDSI values (). Consistent with the spatial drought reconstructions shown in Figs.  and , this North American Southwest regional average reconstruction has high skill, $r$ $=$ 0.90 and CRPSS $=$ 0.57 (we note that for a single time series verification, CRPSS reduces to the skill score version of the mean absolute error, which we have computed here); the reconstruction using pseudo-trees alone has skill values of $r$ $=$ 0.89 and CRPSS $=$ 0.54 (not shown). Despite the high skill of the overall reconstruction, it is challenging to consistently identify the rank of the most severe droughts based on the imperfect reconstruction; this is true not only for the top 5 droughts, but also for the top 10 and 15. This is due to the fact that drought identification and ranking procedures tend to be sensitive to the precise zero crossings of the hydroclimate time series that set the beginning or ending of a drought (thus affecting both the identification and subsequent ranking); for example, the mid- to late-1300s-era drought in Fig.  is cut short in the reconstruction because of a zero crossing of the PDSI values near the middle of the true drought.

Given the established teleconnection between Pacific SST anomalies and hydroclimate variability in North America , we reconstruct the Nino 3.4 index and assess its relationship to North American Southwest multiyear droughts in CESM. As part of the reconstruction experiments, we reconstructed monthly Nino 3.4 with both pseudo-trees and pseudo-corals and as a separate reconstruction just with pseudo-trees. Figure  shows the power spectra of the true Nino 3.4 time series along with the two reconstruction types. As a representative example of the reconstruction, Fig.  shows the monthly Nino 3.4 time series during the most severe drought during the reconstruction (this short interval was shown for clarity). Skill scores for the entire monthly time series for the experiment assimilating both pseudo-trees and pseudo-corals were $r$ $=$ 0.97 and CRPSS $=$ 0.78, while for the experiment with pseudo-trees alone the skill values were $r$ $=$ 0.87 and CRPSS $=$ 0.53. Taking both Fig.  and Fig.  together, the pseudoproxy framework indicates that it is possible to skillfully reconstruct the monthly Nino 3.4 index and to potentially do so without coral proxies. The monthly skill associated with the Nino 3.4 reconstruction is likely because monthly values sufficiently covary with the annual proxy observations such that the reconstruction can provide meaningful monthly states for this variable . The results of the reconstructions using pseudo-trees alone is also particularly important given that coral records are generally short and not extant over more than a few centuries.

As an analysis of the relationship between North American Southwest drought variability and ENSO, Fig.  shows the monthly percent occurrence of El Niño or La Niña states during the 15 most severe multiyear droughts (spanning droughts of 9 to 38 years in length). The El Niño and La Niña states are defined as the months for which the SST anomaly exceeds 0.5 and $-$0.5 ${}^{\circ }$C, respectively. Figure a shows how the percent occurrence for each drought reconstruction compares to the true percent occurrence, while Fig. b shows the distribution of all 15 droughts for the truth and the two reconstructions. This figure illustrates that in the CESM LME simulation, multiyear droughts are associated with an above-average occurrence of La Niña states. This relationship is consistent with observations and historical model simulations, which show that cooler La Niña-like SSTs in the tropical Pacific tend to shift the storm track over North America polewards, leading to a drying of the American Southwest . In a paleoclimate context, it has also been inferred that increased incidences of La Niña events have likely driven megadroughts in western North America over the last 1000 years . Additionally, last-millennium simulations of the CESM model have also been shown to have the strongest relationship between ENSO and drought in the American Southwest compared to other CMIP5-class models . These results thus collectively indicate that drought reconstructions and specific driving dynamics can be successfully investigated within a DA-based reconstruction.

Power spectral density of the true model Nino 3.4 index time series along with the reconstructions based on pseudo-trees and pseudo-corals (purple) and pseudo-trees alone (green).

[Figure omitted. See PDF]

Monthly reconstruction of Nino 3.4 index during the most severe drought period in the reconstruction (a) using both pseudo-trees and pseudo-corals and (b) using only pseudo-trees. The $\mathrm{2}\mathit{\sigma }$ uncertainty estimates are derived from both the Monte Carlo reconstruction iterations and the posterior distribution; see Sect.  for details.

[Figure omitted. See PDF]

(a) Monthly percent occurrence of El Niño or La Niña states during the 15 most severe multiyear droughts. The occurrences are computed as a percentage of the total number of months in a given drought. Dashed lines indicate the average percent occurrence of El Niño (red) or La Niña (blue) in the CESM LME simulation 10. (b) Box plot distribution summaries of the monthly percent occurrence of El Niño or La Niña states shown in panel (a). Light gray lines indicate the means and the darker colors indicate 1 standard deviation of the data, while the lighter colors indicate the 95 % confidence intervals of the mean. The letters LN and EN stand for La Niña and El Niño, while the letters T and C indicate experiments with pseudo-trees or pseudo-corals. “LN” and “EN” by themselves are the true La Niña and El Niño occurrences (gray box plots) corresponding to the distribution of filled circles in panel (a).

[Figure omitted. See PDF]

The latitude locations of the zonal mean multiyear drought and precipitation centroids (PrCent) for both the truth and the reconstruction over equatorial East Africa. Panel (a) shows the centroid locations as a function of drought severity rank, while panel (b) has the same information but more clearly shows the negative correlation between the locations of drought and precipitation centroids.

[Figure omitted. See PDF]

Equatorial East Africa multiyear drought

The next reconstruction illustration we consider is multiyear drought in equatorial East Africa. Specifically, we explore where drought is occurring and its plausible large-scale drivers according to the CESM LME simulations. Similar to the results shown in Sect. , we use drought dynamics in East Africa as an example of what can be reconstructed using the DA method. We additionally provide a dynamical analysis that would be the starting point for a region-specific analysis that could be performed based on a real-proxy reconstruction. The choice of this region is motivated by the observation that the spatial drought index reconstructions, Figs.  and , show skillful drought reconstructions in Eastern Africa and because tropical dynamics and variables are generally better reconstructed than midlatitude variables and dynamics in global, DA-based reconstructions (see, for example, Fig.  and ).

For this analysis we assess how shifts in East African tropical circulation relate to multiyear droughts. An extensive number of paleoclimate studies have broadly interpreted tropical or even extratropical hydrologic changes to shifts in the ITCZ seefor a review. Here we explore how a common indicator of the location of the ITCZ, the zonal precipitation centroid , is related the location of East African droughts. We calculate the precipitation centroid as the center of mass of the zonal mean precipitation within the latitudinal range of 20${}^{\circ }$ N–20${}^{\circ }$ S and the longitudinal range of 28–52${}^{\circ }$ E; as seen in Fig.  this boxed area is centered over Eastern Africa and includes precipitation over both land and ocean. We construct a similar drought location index by computing the centroid of the zonal mean negative SPEI values over the continental African land area within this East African box. We use SPEI for this region because some work indicates that it may be better suited to drought assessments in East Africa than PDSI , though the results shown in this section are not dependent on the choice of drought variable.

Figure a shows the mean latitudinal locations of the 15 most severe multiyear droughts in equatorial East Africa for both the reconstruction and the model truth. Additionally, this figure shows the mean latitudinal location of the precipitation centroid averaged over each of the drought periods. The locations of the two centroids are well matched to the truth in the reconstruction. Furthermore, Fig. b shows that these two centroids are negatively correlated (in both the truth and the reconstructions), indicating that as the average precipitation shifts south in CESM, drought occurs in the north and vice versa. Because the zonal mean precipitation is strongly peaked in the annual mean with more variability in the wings of the zonal mean distribution, small changes in the precipitation centroid can lead to larger changes in the location of drought (Fig. a). Similar to the results presented for North America, these collective results for East Africa further illustrate the utility of DA-based reconstructions and their potential to test specific hypotheses, such as the hypothesis that shifts in the ITCZ can drive multiyear droughts.

Conclusions

Using DA for hydroclimate reconstructions is a new area of research. The aim of this paper was therefore to explore some foundational aspects of global hydroclimate reconstructions using DA. Through a set of pseudoproxy reconstruction experiments, we estimated a realistic upper bound on hydroclimate reconstruction skill for a global pseudoproxy network based on pseudo-tree rings and corals. Unlike most pseudoproxy experiments that use statistical noise perturbations for proxy construction, the pseudo-tree rings and corals were generated here using established proxy system models. These models are physically based and account for multiple climatic factors in proxy development.

Because these reconstruction experiments have been designed as upper-bound skill estimates, the results are most useful in showing where and for which variables the reconstructions lack skill. This framing allows us to highlight what is theoretically possible while also making falsifiable claims about what is impossible to reconstruct given current proxy networks. In our experiments we find that skill is highest in the tropics and in regions local to proxy sampling, in accordance with many previous pseudoproxy and real-proxy reconstructions that focused on temperature or other non-hydroclimate variables e.g.,. For tree rings (which constitute the majority of annually resolved proxies) we also find that local reconstruction skill depends on the moisture or temperature sensitivity of particular tree-ring proxies: moisture-sensitive trees are necessary for skillful PDSI or SPEI reconstructions, and temperature-sensitive trees are similarly necessary for skillful temperature reconstructions. The covariability of moisture and temperature is not sufficiently strong for our DA approach to reconstruct both variables equally well from a given tree. Based on these results, we recommend that future DA-based hydroclimate reconstructions use a spatial mixture of temperature- and moisture-sensitive trees if the goal is to reconstruct both temperature and hydroclimate indices. Importantly, we also performed seasonal and monthly reconstructions and found that it is in principle possible to reconstruct non-growing season or annual mean climate variables in many regions using DA. We also found that it is possible to reconstruct a monthly Nino 3.4 index.

Two regional climate dynamics analyses also explored the key drivers of hydroclimate extremes in the North American Southwest and equatorial East Africa. These regional examples highlight how DA-based reconstructions could in general be used to find the dynamical drivers of hydroclimate extremes. This information-added benefit from DA-based reconstructions is very important given that dynamical information is difficult to simultaneously obtain from a traditional reconstruction. Based on these results, the application of DA to the PAGES2k network may yield skillful hydroclimate reconstructions that will deepen our current understanding of decadal to centennial hydroclimate variability. We find that this approach can potentially yield dynamical insights about regions that have not been previously well characterized, such as Africa.

The code used in the production and analysis of the reconstructions is available upon request from [email protected].

Monthly CESM LME data are available from http://www.earthsystemgrid.org. The pseudoproxies and reconstructions from this paper are publicly available at 10.5281/zenodo.1035402 ().

NJS and JES designed the experiments and prepared the paper. NJS developed the code, ran the experiments, and analyzed the experimental results.

The authors declare that they have no conflict of interest.

Acknowledgements

We acknowledge the CESM1 (CAM5) Last Millennium Ensemble Community Project and the supercomputing resources provided by NSF/CISL/Yellowstone. This work was supported by the NOAA Climate and Global Change Postdoctoral Fellowship Program administered by UCAR's Visiting Scientist Programs. This work was also supported in part by the National Science Foundation under grants AGS-1243204, AGS-1401400, AGS-1602581, and AGS-1602920. LDEO contribution number 8153. Edited by: Stefan Bronnimann Reviewed by: two anonymous referees