2.1 Data assimilation framework

LMR employs ensemble data assimilation (DA) to blend information from proxies and climate model data. DA is performed using a variant of the ensemble Kalman filter, which for our application appears to perform well compared to alternative PDA methods such as particle filters (Liu et al., 2017). The update equation is given by

Here, x_b is the prior state vector, which contains the climate variables to be reconstructed, averaged over an appropriate timescale (here, annual), and x_a is the posterior state vector (i.e., the reanalysis, or reconstruction). The state vector may include scalars, such as climate indices, and/or grid-point data for spatial fields. Vector y contains the assimilated proxy data (i.e., observations), and y_e is a vector containing estimates of the proxies derived from the prior by

where ℋ is the forward model mapping the prior x_b to proxy space (i.e., the PSM; see Sect. 2.4). The innovation, [y–y_e], is the new information from the proxies not already contained in the prior. This new information is weighted against the prior through the Kalman gain matrix:

where B is the prior covariance matrix, R is the error covariance matrix for the proxy data, and H is the linearization of ℋ about the mean value of the prior. Here, Eq. (1) is solved using the ensemble square-root filter (EnSRF) approach of Whitaker and Hamill (2002), in which the ensemble mean and perturbations about the ensemble mean are solved separately. Moreover, R is taken as a diagonal matrix (uncorrelated observation errors) where the diagonal elements represent the error variance for each assimilated proxy record; details on how these are estimated are provided in Sect. 2.4. This allows for serial processing of observations, in which observations are assimilated one at a time. This greatly simplifies the implementation of covariance localization, which is used to control sampling error in the prior covariance. Solutions for the ensemble mean, ${\stackrel{\mathrm{‾}}{\mathit{x}}}_a$, and perturbations, $\mathit{x}{{}^{\prime }}_a$, for the single kth proxy y_k, are obtained from

where y_e,k is the prior estimate of the proxy from Eq. (2), R_k is the diagonal element of R corresponding to proxy y_k, and cov() and var() represent the covariance and variance functions, respectively. The ensemble of updated states is then recovered by combining the posterior ensemble mean and perturbations:

Covariance localization, given by a Schur product denoted by ∘ in Eq. (4b) (i.e., element-wise multiplication), is a distance-weighted filter w_loc on the prior covariance matrix (see, e.g., Hamill et al., 2001). Sections 4.2 and S5 in the Supplement provide more details on localization.

We also use an “appended state vector” approach that avoids the need to recompute Eq. (2) after each proxy is assimilated. The y_e proxy estimates from each record are appended to the state vector x_b:

where the x₁…x_N elements contain the ensemble grid-point data from model variables included in the state (e.g., temperature, precipitation), with N the sum of the number of variables multiplied by the number of grid points, and the $y_{\mathrm{e}}^{\mathrm{1}}\mathrm{\dots }{y}_{\mathrm{e}}^P$ are the ensemble proxy estimates for each of the P proxy records considered. Each of the x₁…x_N and $y_{\mathrm{e}}^{\mathrm{1}}\mathrm{\dots }{y}_{\mathrm{e}}^P$ elements are of dimensions 1×N_ens, where N_ens is the specified size of the ensemble. Hence, x_b is a matrix of dimension $(N+P)\times N_{\mathrm{ens}}$. With such an appended state, the y_e elements in Eq. (6) are updated through Eq. (4b) as any other state variables, eliminating the need to re-evaluate y_e with Eq. (2) once the state has been updated. This simplification is particularly attractive in the context of LMR updates discussed herein as it enables a straightforward implementation of seasonal PSMs (i.e., forward models more accurately representing the seasonal responses of individual proxy records) as discussed in Sect. 2.4. In our implementation for the reconstruction of annually averaged states, the data assimilation procedure follows this general algorithm.

1. The proxy estimates (y_e) are precalculated using Eq. (2) with either annually or seasonally averaged model data as input (i.e., the x_b in Eq. 2).

2. A sample ensemble of annually averaged model states is randomly drawn from a pre-existing simulation to form the main part of the prior state vector (i.e., the x₁…x_N elements in Eq. 6).

3. The precalculated $y_{\mathrm{e}}^{\mathrm{1}}\mathrm{\dots }{y}_{\mathrm{e}}^P$ proxy estimates are added on to form the appended state as shown in Eq. (6). This appended state becomes the x_b in Eq. (1), which is decomposed into an ensemble mean (${\stackrel{\mathrm{‾}}{\mathit{x}}}_b$) and perturbations about the mean ($\mathit{x}{{}^{\prime }}_b$) as shown in Eq. (4b).

4. Proxies forming the y vector are then serially processed, with the updated state, including the proxy estimates, obtained from Eq. (4b). The reanalysis is completed for 1 year once all proxies have been assimilated.

We note here that with a configuration involving seasonal PSMs without the use of an appended state, the vector x_b has to include states with sufficient temporal resolution to allow the calculation of the updated seasonal $y_{\mathrm{e}}^{\mathrm{1}}\mathrm{\dots }{y}_{\mathrm{e}}^P$ proxy estimates. In this scenario, an additional step to the ones listed above is required, involving Eq. (2) using the appropriate seasonally averaged updated states as input. With proxies characterized by a wide range of seasonal responses, this requirement would impose an x_b composed of monthly data which would greatly increase the computational cost of the reanalysis. Reanalysis results would also likely be adversely affected by the larger noise level characterizing data at shorter (i.e., monthly) timescales through its impact on ensemble estimates of prior covariances (see, e.g., Tardif et al., 2015).

Figure 1Locations (a, c) and temporal (b, d) distributions of proxy records available for assimilation (proxies for which linear PSMs calibrated with GISTEMP version 4 are available). Panels (a, b) are used in the prototype version, (c, d) LMR proxy database updated to PAGES 2k Consortium (2017) proxies.

Download

As in Hakim et al. (2016), an “offline” DA approach is used, where the prior ensemble is formed by random draws of time-averaged states from a pre-existing millennium-long model simulation, with the same randomly drawn ensemble members used for every year in the reconstruction of a given reanalysis realization (see Sect. 2.3 below). We note that in the limit of no proxy information, this approach leads to a posterior that reverts to the prior ensemble (see Eq. 1), which randomly samples the model climate and therefore has no skill over the model climatology. This is in contrast to online DA (e.g., Matsikaris et al., 2015; Perkins and Hakim, 2017), where a numerical model is used to dynamically forecast the evolution of climate states from the latest proxy-informed analysis to the following year, when new proxy observations are assimilated. The “offline” approach, introduced by Oke et al. (2002) and Evensen (2003), and used in an ocean DA system by Oke et al. (2005), offers several practical advantages, particularly from a computational cost perspective (Oke et al., 2007). Its use is further justified when model forecasts have limited skill over timescales corresponding to the time interval between updates, as is the case here with global climate models and proxies assimilated on an annual basis. This scenario is further supported by the PDA results of Matsikaris et al. (2015), who show similar performance is achieved with online and offline approaches. From a cost–benefit perspective, the high cost of running ensembles of comprehensive global climate model simulations does not appear justified. However, ongoing research suggests cost-effective online PDA may be achieved by using simplified climate models (Perkins and Hakim, 2017).

2.2 Climate proxies

Our proxy database is updated to the latest PAGES 2k collection (PAGES 2k Consortium, 2017, hereafter PAGES2k-2017). This data set represents the community standard in global proxy observations covering the Common Era (CE) and serves as the core source of proxy information used in our updated reanalysis. PAGES2k-2017 proxies were screened to retain temperature-sensitive records, extensively quality controlled, and described by more metadata compared to previous collections. The additional records assembled by Anderson et al. (2019)¹, consisting in large part of the tree-ring-width records from Breitenmoser et al. (2014) (hereafter B14), are considered as a potential enhancement to proxy information used in our paleoreanalyses (see Sect. 4.3).

As in the LMR prototype (Hakim et al., 2016, hereafter H16), only records with sub-annual to annual resolutions are considered; sub-annual records are averaged to annual. Figure 1 compares the PAGES 2k Consortium (2013) (hereafter PAGES2k-2013) data set used in H16 and the PAGES2k-2017 update. Only records for which a PSM can be established are shown in Fig. 1, defined by proxy records with at least 25 years of (non-contiguous) overlap with calibration data (see Sect. 2.4). Compared to the proxies assimilated in H16, PAGES2k-2017 data provide enhanced spatial coverage in the tropics with additional coral δ¹⁸O and Sr∕Ca records. Additional tree-ring wood-density records from Europe and western North America are also included. The temporal distribution of the total number of records remains similar, except for significant increases in the number of tree-ring-width and coral proxies during 1800–2000 CE and tree-ring wood-density records during 1500–2000 CE.

2.3 Climate model prior information

For all reconstruction experiments reported in this paper, the prior state vector is formed with data from the Coupled Model Intercomparison Project phase 5 (CMIP5) (Taylor et al., 2012) Last Millennium simulation from the Community Climate System Model version 4 (CCSM4) coupled atmosphere–ocean–sea-ice model. The simulation covers years 850 to 1850 CE and includes incoming solar variability and variable greenhouse gases, as well as stratospheric aerosols from volcanic eruptions known to have occurred during the simulation period (see Landrum et al., 2013). The same “offline” DA methodology as in H16 is used, where the prior ensemble is a random sample of model states, with the same sample used for all years of the reconstruction. The sampled states are deviations (i.e., anomalies) from the temporal mean taken over the entire length of the simulation. Therefore, the prior ensemble mean does not contain time-specific information about climate events (e.g., a volcanic eruption) or trends characterizing specific periods (e.g., 20th century warming). Consequently, all trends and temporal structure in reconstructed fields result from information provided by the proxies. Finally, the spatial resolution of prior state variables is reduced from $\mathrm{0.95}{}^{\circ }\times \mathrm{1.25}{}^{\circ }$ of the Last Millennium simulation to a $\mathrm{4.3}{}^{\circ }\times \mathrm{5.7}{}^{\circ }$ Gaussian grid as in H16.

All reconstruction experiments are composed of 51 Monte Carlo assimilation realizations, each using a different randomly chosen 100-member ensemble and 75 % of available proxy records for assimilation. This Monte Carlo sampling over subsets of prior states and proxy records is designed to incorporate uncertainties in covariance estimates derived from model states and uncertainties associated with proxy error estimates. Moreover, we have found that averaging over ensembles from Monte Carlo realizations leads to more accurate results. This is likely the result of averaging over random errors introduced into the reanalysis from a few randomly chosen proxy records with underestimated observation errors. Little sensitivity to the use of 75 % of the proxies for each realization has been found (not shown), while 100 members have been chosen to maintain consistency with H16. In the following, climate reanalyses are taken as the mean over the 100-member DA ensembles and 51 Monte Carlo realizations (i.e., a 5100-member “grand ensemble”).

2.4 Proxy modeling

A critical component of PDA is the mapping of prior climate state variables (e.g., temperature, precipitation from a climate model) to the assimilated proxies (e.g., tree-ring width). This is expressed mathematically by Eq. (2), Sect. 2.1, where the operator ℋ (i.e., the forward model) ideally represents the complete set of processes associated with proxy values, i.e., a comprehensive physically based PSM. This remains a major challenge as the information archive is often complex, involving physical, biological and chemical processes (Evans et al., 2013). Despite recent progress in the development and use of process-based PSMs (e.g., Dee et al., 2015, 2016; Goosse, 2016; Steiger et al., 2017; Acevedo et al., 2017), the focus here is on statistical PSMs, which offer distinct advantages: (1) ease of implementation and flexibility with respect to forward modeling of multiple proxies, regardless of archive types, measurements, units, etc.; (2) observation error statistics for each assimilated record are well defined from the regression (see below); and (3) regressions are formulated on the basis of deviations from the mean over a reference period (e.g., 1951–1980) of the driving climate variable(s), therefore avoiding issues with absolute calibration where climate model bias is problematic, particularly for PSMs having threshold transitions (see, e.g.,  Dee et al., 2016). Statistical PSMs also have distinct disadvantages: (1) PSMs cannot be calibrated without sufficient overlap with calibration data (a threshold of at least 25 overlapping data is imposed); (2) the accuracy of the models depends on the limitations of the calibration data sets (e.g., less reliable analysis over the Southern Ocean and over high-latitude continental areas due to a lack of observations); (3) possible lack of stationarity of the derived relationships established with instrumental-era data; and (4) lack of representation of nonlinear and/or multivariate influences when PSMs are formulated as linear univariate models. Despite these limitations, statistical PSMs provide advantageous capabilities within the context of the LMR and moreover define a baseline to measure future progress with the development of process-based PSMs.

Here, univariate and bivariate statistical PSMs are considered:

and

where y_k are annualized observations from the kth proxy time series, $X_{\mathrm{1}}{}^{\prime },X_{\mathrm{2}}{}^{\prime }$ are anomalies, with respect to the mean over a reference period, of key climate variables (e.g., near-surface air temperature and precipitation) from calibration instrumental-era data sets, β₀ is the intercept, and β₁,β₂ are the slopes with respect to the X₁^′ and X₂^′ independent variables, respectively, and ϵ is a Gaussian random variable with zero mean and variance σ². The overbar in Eqs. (7) and (8) denotes time averages over annual periods, as in H16, or over appropriate seasonal intervals for the seasonal PSMs. Calibration data concurrent with available proxy observations are taken at the grid point nearest the proxy location and the appropriate least-squares solution determines regression parameters $({\mathit{\beta }}_{\mathrm{0}},{\mathit{\beta }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{2}},\mathit{\sigma })$. In this version of LMR, PSM configuration is the same for each proxy category (e.g., univariate for all coral δ¹⁸O, bivariate for all tree-ring-width records).

With the framework described above, the regression-based approach measures the diagonal elements in matrix R through the variance of regression residuals, i.e., R_k=σ². This is a key parameter in PDA as it determines the extent to which the information provided by the proxy is weighted against prior information in the resulting reanalysis. This method provides a sound basis through which assimilated proxy records influence the reanalysis depending on the strength of their relationship to the dependent climate variables. For example, a record with a poor fit to calibration data will be characterized by larger residuals, hence larger observation error variance, and less weight in the reanalysis relative to a record that has a stronger correlation with climate variables. We note that modestly different results are obtained with different observational calibration data sets (see H16).

The calibration data sets used in this study are the NASA Goddard Institute for Space Studies (GISS) Surface Temperature Analysis (GISTEMP) (Hansen et al., 2010) version 4 for temperature and the gridded precipitation data set from the Global Precipitation Climatology Centre (GPCC) (Schneider et al., 2014) version 6 as the source of monthly information on moisture input over land surfaces. The use of precipitation instead of the more traditional Palmer Drought Severity Index (PDSI) to account for moisture is described in more detail in Sect. S4 of the Supplement.

4.1 Proxy system models

The different PSM configurations described in Sect. 2.4 are used in a series of reconstruction experiments using PAGES2k-2017 proxies exclusively. We note that these records have well-defined seasonal metadata.

The impact of seasonal PSMs is first considered with three experiments performed using univariate temperature regression models for (1) annual-mean calibration; (2) seasonality defined by expert metadata; and (3) objectively determined seasonality. Performance is again measured by correlation and CE scores with verification against the Berkeley Earth analysis. Relative to reconstructions with annual-mean PSMs (Fig. 6a and b), the reconstructions with seasonal PSMs (Fig. 6c–f) show improvements in both measures over nearly the entire globe (Fig. 6g–j). Results show a larger improvement for CE (Fig. 6h and j) compared to correlation (Fig. 6g and i), reflecting improvement in both the amplitude of temperature variability and bias. Noteworthy improvements are found in regions with large numbers of tree-ring proxies, such as the western United States, the region around and including Alaska, northern Canada and the western Arctic Ocean, over Scandinavia and the Norwegian Sea, central Asia and over the southern Pacific west of the Antarctic Peninsula (see Fig. 6h). Comparing the differences of correlations and CE in Fig. 6i and j to those shown in Fig. 6g and h reveals that PSMs with objectively derived seasonality contribute positively to skill for the aforementioned regions, especially where tree-ring-width records are most abundant (e.g., North America and Asia).

Figure 6Verification of LMR temperature anomalies against the Berkeley Earth instrumental-era analysis, for experiments using PAGES2k-2017 proxies and univariate PSMs, with contrasting seasonalities. Shown are time series r and CE for (a, b) experiment 1: annual, (c, d) experiment 2: seasonality from the proxy metadata and (e, f) experiment 3: objectively derived seasonality. Differences in skill metrics are also shown (g, h) between experiments 2 and 1, and (i, j) between experiments 3 and 1.

Download

We turn now to the impact of moisture on seasonal TRW PSMs on the reconstructions. Since objectively defined seasonality performs best (i.e., Fig. 6e and f), reconstructions generated with univariate PSMs are used as the reference for measuring skill improvements for modeling TRW records as univariate in either temperature or moisture (abbreviated as “TorM”) (Fig. 7c and d) and for bivariate “temperature and moisture” PSMs (Fig. 7e and f). Improvement over univariate PSMs is apparent for the bivariate approach compared with the univariate “TorM” approach (cf. Fig. 7g, h with i, j, respectively). In the bivariate approach, regions such as western North America and central Asia, where most of the TRW records are found, improve the most in CE, but also over Australia, likely in response to the improved modeling of TRW records in New Zealand and Tasmania. Improvements are also noticeable, through teleconnections with proxy locations in the central Atlantic and southern Indian oceans, and over the eastern North Pacific Ocean. A decrease in skill is present over the midlatitude Pacific Ocean, but this is smaller in magnitude compared with skill enhancements elsewhere.

Figure 7As in Fig. 6 but comparing experiments performed using PAGES2k-2017 proxies with different PSM configurations for tree-ring-width proxies. (a, b) Experiment 1: univariate on temperature for all proxies, (c, d) experiment 2: univariate with respect to temperature or moisture for TRWs and (e, f) experiment 3: bivariate on temperature and moisture for tree-ring widths. Differences in skill metrics are shown (g, h) between experiments 2 and 1, and (i, j) between experiments 3 and 1. All reconstructions are based on objectively derived seasonal PSMs.

Download

Verification of GMT for reconstructions using seasonal PSMs (Table 3) yields a similar interpretation to the spatial verification results. Compared to the consensus of instrumental-era products, we find that the 20th century trend in GMT is overestimated with the PAGES2k-2017 proxy data set if univariate PSMs are used. This is particularly the case with annual PSMs. Better agreement is obtained when seasonal bivariate PSMs are used to model TRW proxies. The representation of GMT interannual variability as measured by verification of the detrended GMT is also improved with seasonal PSMs, particularly for the CE metric. Similar to spatial verification results, PSMs with objectively derived seasonality and bivariate TRW modeling have GMT reconstructions with consistently higher skill scores.

Table 3Summary of instrumental-era verification results for reconstruction experiments performed with various PSM configurations. Verification scores shown are the trend over the 20th century (in K/100 years), r and CE for the annual GMT and detrended GMT verified against the consensus of instrumental-era analyses. The GMT trend in the consensus of instrumental-era analyses is 0.56 K/100 years.

Download Print Version | Download XLSX

We recognize that the previous evaluation relies on comparisons with observation-based products covering the same time period as the data used to calibrate the statistical PSMs. To test the sensitivity of the results to the calibration period, we conduct additional independent instrumental-era calibration–validation experiments where PSMs are calibrated over a subset of the instrumental-era period and reconstructions are evaluated with data not used in calibration. Results from these experiments, described in Sect. S3 in the Supplement, confirm the main results and conclusions drawn here on the superiority of seasonal PSMs relative to those calibrated with annual averages and the use of bivariate models for TRW proxies.

We now examine results from an evaluation performed in proxy space using proxies withheld from assimilation as in Sect. 3. Results for both the PSM calibration and pre-calibration periods are shown in Table 4. Differences among the various experiments suggest the superiority of the seasonal (with objective seasonality) PSMs as skill scores consistently rank among the highest among all experiments for both calibration and pre-calibration periods. The reconstruction using univariate annual PSMs shows the weakest verification statistics, confirming the verification based on instrumental-era analyses. Finally, use of bivariate seasonal PSMs for TRW records is also suggested from proxy validation results, as larger correlations and ΔCE are obtained with this configuration.

Table 4Verification of LMR reconstructions against independent (withheld from assimilation) proxies for experiments using various PSM configurations. Skill scores shown are the median of distributions for r, the fraction of proxy records characterized by a positive ΔCE (%+CE) and the median of the ΔCE distribution. Statistics are compiled over 51 Monte Carlo realizations for two distinct periods: 1880–2000 (PSM calibration period) and 0–1879 (pre-calibration period).

Download Print Version | Download XLSX

4.2 Covariance localization

One approach to managing sampling error in ensemble data assimilation is through spatial covariance localization. Localization is applied to minimize the adverse impact of spurious covariances at large distances from a proxy location, which results from sample error in finite ensembles (Hamill et al., 2001). If localization is not applied, spurious covariances allow proxies to affect remote locations, which adversely affects the quality of the analysis. On the other hand, too-short localization length scales reduce the useful information that can be derived from the proxies. Therefore, a balance is sought between minimizing sampling noise versus retaining useful proxy information.

We use the Gaspari–Cohn (Gaspari and Cohn, 1999) fifth-order polynomial with a specified cut-off radius for the localization function (w_loc in Eq. 4b). See Sect. S5 for information on the characteristics of w_loc. A series of reconstructions is performed with a wide range of localization length scales. As with previous experiments, 51 Monte Carlo realizations are carried out, each with 100 ensemble members assimilating 75 % of proxy records. Results from the instrumental-era verification scores previously described are summarized in Table 5. We observe that the GMT trend is underestimated and verification scores are significantly reduced when “too-small” localization radii are used, indicating the information on temperature provided by some proxy records is not properly incorporated in the reanalysis. In contrast, the trend is overestimated and verification scores are generally reduced without covariance localization. This is particularly the case for the CE score for the detrended GMT, sensitive to the amplitude in interannual variability. This skill measure is maximized for localization radii within the 15 000 to 25 000 km range. A localization radius at the upper end of this range (25 000 km) is preferable, as results from the other verification scores suggest that a skillful reconstruction is obtained with this covariance localization configuration. See Fig. S4 for an example where the 25 000 km localization function is applied to a proxy record located in California, United States. We note that the optimal localization radius depends on a number of factors, such as ensemble size, the observation network and observation error characteristics.

Table 5The 20th century trend of GMT, r and CE, for the annual and detrended GMTs, as well as the global mean of the spatial (i.e., grid point) r and CE of reconstructed temperature verified against the consensus of instrumental-era analyses for reconstruction experiments performed with covariance localization using various localization cut-off radii (L_R). Verification statistics for an experiment without covariance localization are also shown for comparison. Results from the prototype are shown for reference. The GMT trend in the consensus of instrumental-era analyses is 0.56 K/100 years.

Download Print Version | Download XLSX

4.3 Proxy data sets

Here, we explore the impact of adding the large number of proxies from Anderson et al. (2019) (hereafter A19), which include the tree-ring-width chronologies from Breitenmoser et al. (2014) (hereafter B14), not strictly screened for climate sensitivity in contrast to the PAGES2k collection. Duplicate records between data sets are identified (based on correlation between co-located records and cross-referencing metadata) and eliminated. Priority is given to records found in the PAGES2k collection (see A19 for more details). Figure 8 shows the spatial and temporal distributions of the B14 records, which reveals enhanced coverage over eastern North America, southern Europe, boreal Eurasia and southern South America. Other additions, totaling 94 records, provide additional records in the tropics (23 coral records) and an enhanced number of ice core records concentrated over Greenland and the eastern Canadian Arctic (37 records) and Antarctica (26 records in West Antarctica and Dronning Maud Land). A few lower-latitude ice core records (six records) are also added in the Peruvian Andes and Tibetan Plateau, along with two higher-latitude lake core records. From a temporal perspective, the addition of the B14 tree-ring-width records contributes a notable number of additional proxies back to 1000 CE, more than double the number of records available for assimilation from 1500 CE onward, up to a 4-fold increase during the 19th and 20th centuries.

Figure 8Locations (a) and temporal distributions (b) of the additional proxies from Anderson et al. (2019) considered for assimilation, including the tree-ring chronologies from Breitenmoser et al. (2014). As in Fig. 1, only records available for assimilation (proxies for which regression-based PSMs can be calibrated) are shown.

Download

In order to measure the impact with the best configuration, the reconstruction experiments reported in this section are carried out using seasonal PSMs with objectively derived seasonality for all records, with a bivariate formulation on temperature and precipitation for all TRW proxies and univariate on temperature for all other proxies. The baseline reconstruction uses the PAGES2k-2017 proxies (as in Sect. 3), which we compare to results first obtained with the addition of the B14 TRW records and finally with the further addition of the coral, ice and lake core records from A19 (i.e., the full proxy database). Other trial reconstructions performed with the vastly expanded proxy network, not reported here, show that a well-calibrated GMT ensemble is obtained with a covariance localization cut-off radius of 25 000 km. Next, we compare reconstruction results from this configuration to the baseline reanalysis.

Differences in correlation and CE associated with the addition of the B14 collection over the PAGES2k-2017 proxies show skill improvements in temperature reconstructions over the continental United States and Mexico, Europe and the southern edge of the Tibetan Plateau (see Fig. 9g and h). Through the influence of significant spatial covariances with the added records, assimilation of the additional TRW records also leads to improved temperature skill over remote areas of the midlatitude Pacific and northern Atlantic oceans. The addition of records described in A19 has minimal additional impact overall, with the exception of modest increases in correlation and CE over Greenland (see Fig. 9i and j).

Figure 9As in Fig. 6 but comparing experiments performed with different proxy networks: (a) r and (b) CE for experiment 1: PAGES 2k Consortium (2017) proxies only, (c, d) experiment 2: with the addition of tree-ring chronologies from Breitenmoser et al. (2014) and (e, f) experiment 3: with all proxies in the updated LMR database. The differences in correlation and CE between experiments 2 and 1 are shown in panels (g, h), respectively, and between experiments 3 and 2 in panels (i, j). Notice the latter is different from Fig. 6, where differences between experiments 3 and 1 are shown.

Download

Hydroclimate verification is defined by a comparison of the reconstructed PDSI with the Dai (2011) product. We note here that the reconstruction is not directly related to the PDSI product used for verification, as TRW forward models were calibrated on precipitation and not on PDSI as in Steiger et al. (2018). A comparison of the reconstructed PDSI between the prototype³, the updated reanalysis of Sect. 3 and a reconstruction carried out with the B14 TRW records and the additional coral, ice and lake core records (i.e., the full database) is shown in Fig. 10. The PDSI is slightly improved in the updated reanalysis compared to the prototype (Fig. 10g and h). Enhanced skill is noticeable over western North America and over eastern Europe and Asia to a lesser degree. Decreased skill is found over the central plains of North America and along a narrow band along the Siberian Taiga. The impact of adding the Anderson et al. (2019) records is mostly found over the eastern part of the United States and over western Europe (Fig. 10i and j). Finally, we note that this impact is due entirely to the B14 TRW records, as the additional coral, ice and lake core records from A19 do not significantly affect the PDSI reconstruction skill (from results of additional reconstruction experiments carried out to isolate this impact; not shown).

Table 6As in Table 4 but statistics compiled for tree-ring wood-density (MXD) proxies only, and for experiments using the PAGES2k-2017 proxies only, PAGES2k-2017 with the addition of all proxies from A19 (all proxies) and PAGES2k-2017 plus only a subset of A19 records obtained after removing all but 188 TRW records from B14 (B14 subset). See text for selection details. Skill scores are the median of r distributions, the fraction of proxy records characterized by a positive ΔCE (%+CE) and the median of ΔCE distributions. Statistics are compiled over the 51 Monte Carlo realizations for the following periods: 1880–2000 (PSM calibration period) and 1600–1879 (pre-calibration period with a significant number of MXD records and covering a significant portion of the Little Ice Age).

Download Print Version | Download XLSX

Figure 10Similar to Fig. 9 but comparing PDSI reconstructions against the Dai (2011) analysis for experiments performed with different proxy networks: (a) correlation and (b) CE for experiment 1: prototype reanalysis from H16, experiment 2: PAGES 2k Consortium (2017) proxies, (e, f) experiment 3: with further the addition of tree-ring chronologies from Breitenmoser et al. (2014) and the coral, ice and lake core records from Anderson et al. (2019) (i.e., the full proxy database). The differences in correlation and CE between experiments 2 and 1 are shown in panels (g, h), respectively, and between experiments 3 and 2 in panels (i, j).

Download

Examining the differences between reconstructions over the entire Common Era (Fig. 11), we see a significantly modified Northern Hemisphere temperature (NHMT) resulting from the assimilation of the additional proxies. A generally warmer NHMT is obtained throughout the Common Era but most significantly during the LIA, worsening the agreement with reconstructions from other studies shown in Fig. 2. A noticeable loss of variability is observed, confirmed by comparing spectra from both experiments (Fig. 11c). This loss of variability in the reconstruction using all proxies occurs at nearly all scales, underlining an adverse impact from assimilating B14 tree-ring-width proxies.

Figure 11(a) Northern Hemisphere temperature (NHMT) grand ensemble mean (solid lines) and 5th-–95th percentile range (shading) from experiments performed with PAGES2k-2017 proxies (in blue) and with the addition of proxies from Anderson et al. (2019) (in red). (b) Spectra of NHMT grand ensemble mean from both experiments (solid lines), along with the χ² 95 % highest density regions (shading).

Download

We now turn to verification in proxy space, which is the only source available prior to the instrumental period. Proxy estimates from reanalyses (estimated using the appropriate PSM) are compared directly to proxy observations. Here, reanalysis skill is assessed using independent (the 25 % withheld from assimilation) proxies. We further restrict our analysis to verification against tree-ring wood-density proxies, as they are among the most reliable recorders of temperature in our database, as evidenced by the generally better fits to calibration temperature data obtained when calibrating the univariate PSMs. Also, these proxies provide good temporal coverage of the latter portion of the LIA into the industrial period, as shown in Fig. 1. The results, presented in Table 6, show distinctly larger skill scores for the experiment using PAGES2k-2017 proxies only compared to when all proxies are assimilated. Improved skill is observed for both periods of interest. Results from a third reconstruction experiment are also presented, where only a small fraction of B14 records are assimilated (B14 subset experiment in Table 6). A total of 188 records (out of the 2156 available) have been selected on the basis of their strong relationship to calibration temperature and precipitation data as determined from the correlation coefficient characterizing bivariate PSMs. Records with a calibration correlation above 0.6 are found to be located for the most part over the United States. Proxy verification results indicate an increase in skill in the representation of tree-ring wood-density proxies, as indicated by skill metric values only slightly lower than in the PAGES2k-2017 experiment. Spatial verification of temperature and PDSI (not shown) also suggests that some of the skill enhancements shown in Fig. 10i and j are retained even when this small fraction of the B14 records is considered. This suggests that the issues with the assimilation of the B14 records identified above can possibly be mitigated while maintaining some of the skill they provide toward enhanced temperature and hydroclimate reconstructions in local regions. Optimal selection of these records requires further careful attention and could serve as the basis for future efforts.