Using a process-based dendroclimatic proxy system model in a data assimilation framework: a test case in the Southern Hemisphere over the past centuries

Currently available data-assimilation-based reconstructions of past climate variations have only used statistical proxy system models to make the link between climate model outputs and indirect observations from tree rings. However, the linearity and stationarity assumptions of the statistical approach may have limitations. In this study, we incorporate the process-based dendroclimatic model MAIDEN into a data assimilation procedure using the reconstruction of near-surface air temperature, precipitation and winds in the midlatitudes of the Southern Hemisphere over the past 400 years as a test case. We compare our results with a data assimilation approach including a linear regression as a proxy system model for tree-ring width proxies. Overall, when compared to instrumental data, the reconstructions using MAIDEN as a proxy system model offer a skill equivalent to the experiment using the regression model. However, knowing the advantages that a process-based model can bring and the improvements that can still be made with MAIDEN, those results are promising.

1 Introduction

Indirect observations of climate from proxies such as tree-ring width (hereafter TRW) or isotopic content in ice cores and coral provide information on past climate variability beyond the instrumental era (Jones et al., 2009). Those proxy records have been used in numerous studies devoted to climate variability over the last millennium. This has been classically achieved by using the empirical (i.e., statistical) relationship between climate variables of interest (e.g., surface air temperature) and proxy observations (e.g., Mann et al., 1999; Luterbacher et al., 2004; Mann et al., 2008; Büntgen et al., 2011; Neukom et al., 2011; Tingley and Huybers, 2013; Wilson et al., 2016). Tree rings represent one of the main sources of information for those studies as they are the most available and spatially distributed proxy to reconstruct past climates at high temporal resolution (e.g., Fritts, 1976; Jones et al., 2009; Cook et al., 1999; Esper et al., 2002; Wilson et al., 2016; Anchukaitis et al., 2017; Esper et al., 2018).

Although the network of TRW records is well developed, mostly in the extratropical latitudes of the Northern Hemisphere (PAGES 2k Consortium, 2017), these records only provide insight about past climates at a specific location. In order to spread the local information to the large spatial scale, statistical methods have been developed to provide spatially gridded climate reconstructions of the past centuries, generally on the basis of linear regressions between the network of proxy records and the climate field of interest, calibrated over the instrumental period (e.g., Fritts et al., 1971; Mann et al., 1998, 2007; Cook et al., 2010; Wang et al., 2015; Anchukaitis et al., 2017). In this so-called transfer function approach (Fritts et al., 1971), the climate field of interest is expressed as a function of the tree-ring proxy, i.e., in the inverse mode (as opposed to the natural direction in which tree growth responds to climate). However, while this assumes an invertibility of the relationship between climate and tree growth, this may not always be the case. Additionally, the use of linear regression is based on the assumption that the relationship between climate and tree growth is linear and stationary over time. However, many studies have highlighted the shortcomings associated with these assumptions for TRW records (Wilmking et al., 2020; Briffa et al., 1998; Wilson and Elling, 2004; Wilson et al., 2007; D'Arrigo et al., 2008; Guiot et al., 2014; Babst et al., 2018).

On the other hand, global climate model simulations offer a complete spatial picture of the past climate variability. Nevertheless, these simulations, if used alone, cannot follow the observed climate evolution because of the large contribution of internal climate variability that is unpredictable over those timescales and thus cannot be reproduced by models, as well as because of potential biases. Over the last decades, data assimilation (DA) has been increasingly used in paleoclimatology to bring together the information from both indirect climate records from proxies and climate models to provide a more complete picture of the climate changes (Goosse et al., 2010; Widmann et al., 2010; Goosse et al., 2012; Steiger et al., 2014; Franke et al., 2017; Steiger et al., 2018; Tardif et al., 2019).

Initially, DA-based reconstructions were developed using temperature reconstructions derived from the statistical relationship between temperature and proxy observations, not the proxy time series themselves (e.g., Crespin et al., 2009; Goosse et al., 2010). More recently, the use of forward proxy system models (hereafter, PSMs) has emerged in order to directly assimilate the proxy time series. This gives the opportunity to consider the relationship between the proxy and the climate in the natural direction, as opposed to inverse models, and to take into account the combined influence of multiple climate variables on the proxy data. Specifically, PSMs – one specific form of the observation operator in the DA framework – make the link between the outputs of the climate model included in the DA procedure and the assimilated proxy observations, for example TRW (Evans et al., 2013; Dee et al., 2016). To date, reconstructions based on data assimilation using observed tree-ring data have only included statistical PSMs (univariate or multivariate regressions; e.g., Franke et al., 2017; Steiger et al., 2018; Tardif et al., 2019). Contrary to statistical models, process-based dendroclimatic PSMs are able to account for the complexity of the relationship between climate and tree-ring proxy data by explicitly simulating the biological processes governing the climate dependency of tree-ring growth. In this regard, they may overcome some of the limitations of statistical models. However, so far, they have never been used in a DA procedure with actual tree-ring proxy data, but only with pseudo-proxies (i.e., synthetic proxy data; Dee et al., 2016; Acevedo et al., 2017; Steiger and Smerdon, 2017). In these studies, the VS-Lite model (Tolwinski-Ward et al., 2011) is forced by climate model outputs in order to provide synthetic tree-ring proxy data. These synthetic proxy data are then incorporated in a DA framework, providing a controlled environment to evaluate the DA procedure and resulting reconstruction skill. In those experiments, using the nonlinear VS-Lite model as a PSM for tree-ring data has the potential to improve the quality of the reconstruction compared to the traditional statistical linear approach of tree-ring models. However, VS-Lite cannot be considered entirely process-based as it does not include any biological processes to drive tree growth but is rather based on the implementation of the principle of limiting factors with threshold growth response functions (Tolwinski-Ward et al., 2011).

In this study, we introduce for the first time a process-based dendroclimatic model in a data assimilation procedure with actual tree-ring width records. Among all the available mechanistic tree-growth models (e.g., Misson, 2004; Dufrêne et al., 2005; Vaganov et al., 2006; Drew et al., 2010; Tolwinski-Ward et al., 2011), we chose to work with the ecophysiological model MAIDEN (Modelling and Analysis In DENdroecology; Misson, 2004). MAIDEN includes the atmospheric CO₂ concentration among other inputs, which represents a major advantage for taking into account the potential effect of the recent exponential increase in CO₂ on tree growth (Myhre et al., 2013; Boucher et al., 2014). In Rezsöhazy et al. (2021), the MAIDEN model has been successfully applied over the last century to the PAGES 2k TRW database by comparing the simulated annual quantity of carbon allocated to the stem with the observed tree-ring widths (both normalized) using high-resolution climate data as inputs of the model. The model has been calibrated over the second half of the 20th century using the Bayesian procedure with Markov chain Monte Carlo sampling elaborated in Rezsöhazy et al. (2020) and successfully verified over the first half of the century at 20 % of the PAGES 2k tree-ring sites. The studies have highlighted the potential of MAIDEN to be used as a PSM for DA-based reconstruction of past climate variability compared to a simple model like VS-Lite (Tolwinski-Ward et al., 2011).

Here, we perform a comparative analysis on the impact of using a simple regression model or a process-based dendroclimatic model as PSMs for TRW records in the DA framework. The goal is to evaluate if and how a complex process-based model of tree growth like MAIDEN can contribute to the improvement of large-scale DA-based reconstructions of past climate variability compared to the regression model classically used so far. Because MAIDEN has been shown to be skillful in different regions of the Southern Hemisphere (SH; Rezsöhazy et al., 2021), we focus on that region and provide climate reconstructions of near-surface air temperature, precipitation and winds in the middle to high latitudes of the SH over the past 400 years. The Southern Hemisphere also constitutes an interesting study case, as few reconstructions have been performed in this region so far (Neukom et al., 2011, 2014). The DA experiments are based on both TRW and ice core (δ¹⁸O and snow accumulation) proxy data to ensure the consistency of the reconstructed large-scale circulation pattern by assimilating proxy records at middle and high (i.e., Antarctica) latitudes of the Southern Hemisphere and thus avoiding only reconstructing small-scale features. They both represent the best available continental annually resolved proxies in the SH (PAGES 2k Consortium, 2017). While the TRW records require the use of a PSM, the ice core records can be directly compared with the outputs of an isotope-enabled climate model.

First, DA experiments are performed and evaluated against state-of-the-art gridded climate datasets over the last century to evaluate the impact of the assimilation of the tree-ring data and the associated PSMs on the skill of the reconstruction. In Sect. 4.1, we compare the performance of the DA-based reconstructions using a regression- or a process-based dendroclimatic PSM in our framework over the last century. Sensitivity experiments are then performed in Sect. 4.2 to identify the contribution of tree-ring width proxies in our data assimilation framework using a regression- or a process-based dendroclimatic PSM as well as to quantify the impact of the uncertainty of the records. Secondly, the DA-based reconstructions performed with a statistical or process-based dendroclimatic PSM are compared to each other over the past 4 centuries in Sect. 4.3. To match the seasonality of tree growth in the Southern Hemisphere, all the analyses in this study are performed on a July–June year (e.g., the year 1992 corresponds to July 1992–June 1993) instead of the usual January–December year, which is suited for trees growing in the Northern Hemisphere.

2 Data assimilation

In paleoclimatology, the objective of data assimilation is to optimally combine the information from the physics of the climate system as included in climate models and indirect observations of climate from proxies (Kalnay, 2003; Goosse et al., 2010; Widmann et al., 2010; Steiger et al., 2014; Franke et al., 2017; Hakim et al., 2016). More specifically, the DA procedure is based on Bayes' theorem (van Leeuwen, 2009): starting from the prior estimate of the state of the climate system provided by the model, DA produces a reconstruction of the climate system that is as accurate as possible (i.e., the posterior) given the information provided by the climate observations. The updated estimate of the climate system is often called a reanalysis. One of the advantages of DA is its ability to propagate the information brought by the observations both in space and between the climate variables by relying on the spatial covariance of a specific variable and on the covariance between different variables in the model. In other words, DA gives the opportunity to reconstruct the climate at places where there are no proxy sites and for variables for which we do not have any observations to assimilate.

In this study, we use an offline DA approach in opposition to the online DA approach. The offline approach samples a pre-existing ensemble of climate model simulations to build the prior state of the climate system. This has the great advantage of being considerably computationally time-saving. In this configuration, the different years of the sampled ensemble are assumed to be independent of each other. This assumption is reasonable if the predictability of the variables of interest is smaller than the temporal resolution of DA (i.e., 1 year in our study), as is the case for atmospheric variables (Okazaki et al., 2021). With the online approach, the prior state of the climate system is sequentially updated. In this case, the years are not assumed to be independent of each other. Instead, the update made at the previous time step is transferred to the next one. When focusing on the atmosphere, as in our study, the online approach is not expected to outperform the offline approach (Matsikaris et al., 2015; Okazaki et al., 2021).

Here, the offline approach is based on a particle filter method described in Sect. 2.1. In addition, the framework of DA includes three main components: (1) prior climate model simulations (Sect. 2.2), (2) indirect observations of climate provided by proxy data (Sect. 2.3), and (3) proxy system models to make the link between climate model outputs and indirect observations of climate (Sect. 2.4). The error associated with the data used in the DA procedure is also considered and defined in Sect. 2.5. Finally, in Sect. 3, we elaborate the DA experimental design and diagnostics of this study.

2.1 Particle filter

In this study, the offline DA procedure employs a particle filter as in Dubinkina et al. (2011). The particle filter (van Leeuwen, 2009) describes the prior of each climate variable by a probability density function, which is represented discretely on the basis of independent states of the model (Sect. 2.2), i.e., the particles (Fig. S1 in the Supplement). These model states are sampled at a chosen frequency in the model simulations that determines the number of sampled particles. For instance, for a chosen frequency of 1 year, we sample a particle or model state in the model simulation every 1 year. It will be characterized later in Sect. 3. The prior stays the same throughout the DA procedure. As the observations are annually resolved (see Sect. 2.3), the time step of the DA is also annual. As a consequence, the particles are defined as the annual mean of the model state variables.

For each year of the reconstruction period (i.e., the period covered by observations), all the particles are compared to the available observations through the use of a PSM (Sect. 2.4). This allows estimating the likelihood of each particle knowing the observations at that time. This computation also takes into account all known sources of errors (see Sect. 2.5). The likelihood then determines the weight given to each particle. The weight given to a particle will be higher if the particle is closer to the observations and lower if further from the observations. In this way, starting from a prior distribution wherein all the particles have the same weight, the particle filter produces a posterior distribution wherein the weights of the particles are redistributed according to the observations for each year of the reconstruction and independent of the location (Fig. S1). The method does not make any assumptions on the prior and posterior distributions of the state of the climate system, and it can be used without a priori knowledge of the model, which is a clear advantage of the particle filter. Additionally, in this framework, the individual particles for a given time step are dynamically consistent climate states. However, it is important to keep in mind that even though atmospheric states are almost independent on an annual timescale (as in this study), the dynamics of the ocean and cryosphere components of the climate system involve very long timescales, and their states also influence the atmospheric variables. At the end of the DA procedure, the weighted mean of the particles for each year over the assimilation period provides an estimate of the climate variables (i.e., the ensemble mean). Since dynamical consistency in space and between variables involves nonlinear equations, the weighted ensemble mean is therefore not fully dynamically consistent. Finally, the weighted standard deviation of the particles is used to estimate the range of the reconstruction (i.e., the ensemble spread).

2.2 Climate model simulations

We use the ensemble of three simulations performed with the isotope-enabled Community Earth System Model version 1 (iCESM1; Brady et al., 2019; Stevenson et al., 2019; publicly available at https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.CESM_CAM5_LME.html; last access: 12 May 2021) as a basis for building the prior of our data assimilation experiments. The iCESM1 is a coupled atmosphere–ocean–sea ice model at 1.9^∘ latitude and 2.5^∘ longitude resolution for the atmosphere component and around 1^∘ for the ocean–sea ice component. The model simulations include anthropogenic (aerosol emissions and land use changes) and volcanic aerosol forcings, greenhouse gas emissions, solar irradiance variations, and Earth's orbital changes (Brady et al., 2019; Stevenson et al., 2019). The annually resolved iCESM1 simulations used in this study span the July 850–June 2005 CE time period. This corresponds to the time span of the simulations used to sample the particles during the data assimilation procedure by contrast with the time span of the reconstructions (i.e., the past 4 centuries). For our DA experiments (Sect. 3), all the particles are sampled in the three iCESM1 members for a total of 3465 particles.

2.3 Proxy data

In this study, we make use of three types of proxy data from two types of archives: the first is tree-ring width (Sect. 2.3.1) and the second is δ¹⁸O (i.e., the ratio of stable isotope oxygen in the ice core) and snow accumulation from ice cores (Sect. 2.3.2). They represent the best available continental proxies in the SH continents (South America, Australia, Tasmania, New Zealand and Antarctica) for reconstructing the past climate at high resolution (PAGES 2k Consortium, 2017).

2.3.1 Tree-ring width data

The PAGES 2k database (hereafter PAGES2k2017) includes 692 proxy records (PAGES 2k Consortium, 2017) selected for reconstructing the temperature of the last 2 millennia. The records are from a variety of archives (tree rings, ice cores, lake sediments, coral, speleothems, documentary evidence, etc.). TRW records represent a bit more than half of the database with 354 records, among which 12 are located below 30^∘ S, in New Zealand (four records), Tasmania (four records) and South America (four records). The previous version of the PAGES 2k TRW database (PAGES 2k Consortium, 2013; hereafter PAGES2k2013) contains 28 TRW records below 30^∘ S. Nine records have been excluded from the recent version of the PAGES 2k database because they were too short (less than 300 years). Those records are also included in this study and are all located in South America. The location of the 21 TRW records from the two databases can be seen on Fig. S2 (blue stars). The location of the 15 TRW records that are actually used in this study can be observed in Fig. 1 (see Sect. 2.4.1 and 2.4.2 for more details on the included sites).

Figure 1 Location of the 15 tree-ring width (TRW) records (blue stars) from the PAGES 2k databases (PAGES 2k Consortium, 2013, 2017) used in this study: the 49 ice core snow accumulation records (purple triangles) from Thomas et al. (2017) and Medley et al. (2018) as well as the 29 δ¹⁸O records (green diamonds) from Stenni et al. (2017) (Sect. 2.3). Background map from Hunter (2007).

Here, we use the already standardized tree-ring series directly derived from the PAGES 2k compilation without any modification. As the tree-ring series are from different studies (referred in the PAGES 2k database metadata; PAGES 2k Consortium, 2017), the standardization technique used can differ from one tree-ring series to another. Finally, the TRW records are normalized (i.e., with a null mean and a unity standard deviation) relative to the 1900–2000 CE time period for use with the MAIDEN model (Sect. 2.4.1), while for the regression-based model (Sect. 2.4.2), the anomalies are computed relative to the 1900–2000 CE time period.

2.3.2 Ice core data

In addition to tree-ring proxy records, which are the focus of this paper through the use of a process-based PSM (Sect. 2.4.1), we also assimilate other proxy data, i.e., ice core records, to cover a vast area of the Southern Hemisphere (i.e., Antarctica) where no tree-ring proxy data can be found in order to provide a comprehensive reconstruction of the Southern Hemisphere climate.

The ice core data selected here are the same as in Dalaiden et al. (2021), a study focusing on the mechanisms behind the West Antarctic climate changes. This includes the 48 annually resolved ice core snow accumulation records in Antarctica from the database of Thomas et al. (2017) spanning the last millennium (Medley and Thomas, 2019). The B40 ice core snow accumulation record of Medley et al. (2018) has also been added as this record was published after the publication of the Thomas et al. (2017) database. Finally, as in Dalaiden et al. (2021), we use the 29 annually resolved δ¹⁸O records in Antarctica from Stenni et al. (2017). Figure 1 represents the location of the snow accumulation and δ¹⁸O records (purple triangles and green diamonds, respectively). They are mostly situated in West Antarctica. The ice core records are averaged over a 500 km grid in order to reduce the nonclimatic noise in the observations. This represents a maximum of 20 snow accumulation and 13 δ¹⁸O time series over the last 4 centuries. The anomalies of the ice core data are computed relative to 1900–2000 CE.

2.4 Proxy system models

Proxy system models, also referred to as observation operators, are a crucial component of the DA framework that makes the link between the climate state variables of the prior, such as near-surface air temperature and precipitation, and the assimilated indirect observations from proxies, here tree-ring width. In this study, we make use of two types of tree-ring PSMs at opposite ends of the PSM spectrum: the first is a process-based ecophysiological model, MAIDEN (Sect. 2.4.1), and the second is a linear regression model (Sect. 2.4.2). If we neglect potentially important elements such as post-depositional effects, the ice core records can be directly compared with the climate model. In particular, the difference between precipitation and sublimation–evaporation (P–E) from iCESM1 is compared to the snow accumulation observations. It has already been shown in different studies that the main contributor to snow accumulation variability is indeed P–E, especially when working at large spatial scale (Agosta et al., 2019; Souverijns et al., 2018; van Wessem et al., 2018). Finally, the ratio of stable isotope oxygen δ¹⁸O is simulated in iCESM1 and can be directly compared to δ¹⁸O records.

2.4.1 MAIDEN

MAIDEN (Modelling and Analysis In DENdroecology; Misson, 2004; Gea-Izquierdo et al., 2015; Gennaretti et al., 2017) is a process-based dendroclimatic model which simulates the course of photosynthesis and the allocation of carbon to different pools (stem, roots, leaves, storage) at a daily time step as a response to climatic drivers (daily maximum and minimum temperature, precipitation, and atmospheric CO₂ concentration). The annual quantity of carbon which is allocated to the stem (hereafter Dstem) is assumed to be proportional to tree-ring growth and comparable to TRW data after normalization. Therefore, this variable is used for the comparison with the observed TRW. In practice, as previously discussed in Rezsöhazy et al. (2020, 2021), the comparison of the tree-growth observations with the forest carbon accumulation simulated by MAIDEN implies some limitations. For instance, Gennaretti et al. (2018) computed a wood biomass index directly comparable to what MAIDEN simulates. This involves using both tree-ring width and density (earlywood and latewood fractions) data, while the latter is less broadly available. Finally, we use a combined version of the model from Gea-Izquierdo et al. (2015) and Gennaretti et al. (2017), developed by Fabio Gennaretti (unpublished).

We must first determine where the MAIDEN model can be skillfully applied among the 21 available PAGES 2k sites (Sect. 2.3.1). MAIDEN was calibrated at the 12 SH sites from PAGES2k2017 in Rezsöhazy et al. (2021) driven by the daily maximum and minimum temperature as well as precipitation from the Global Meteorological Forcing Dataset for land surface modeling (v2) (http://hydrology.princeton.edu/data.php, last access: 5 May 2021; Sheffield et al., 2006) at 0.5^∘ resolution (hereafter GMF) over 1950–2000 CE. The model was validated using the same dataset over the 1901–1949 CE time period. The calibration method is based on a Bayesian procedure with Markov chain Monte Carlo sampling using the DREAMzs algorithm (Hartig et al., 2019). This Bayesian approach is used to calibrate the most sensitive parameters of the model from the photosynthesis and carbon allocation modules (18 parameters). MAIDEN was found to be robustly calibrated with significant (95 % confidence level) calibration and verification correlations greater than 0.3 at five sites south of 30^∘ S. The same methodology was used here to calibrate and validate MAIDEN for the nine records of PAGES2k2013. One site was found to be well calibrated and validated. The six sites for which MAIDEN calibration has successfully passed those tests will thus be incorporated into the DA procedure with MAIDEN. They can be seen in Fig. S3 (orange circles). In particular, there is no site in Tasmania.

Data assimilation in principle does not require such a selection based on correlation because the confidence we give to each record during the DA procedure is proportional to the observation error assigned to each record (Sect. 2.5). The potentially low correlation of some simulated tree growth with the observations is thus taken into account in the observation error and in the DA. In the Last Millennium Reanalysis project (Hakim et al., 2016; Tardif et al., 2019), for instance, such a selection is not made. In this study, we apply a selection using a low value of the correlation as a threshold, since tests have shown that including records with a very low correlation can increase the noise and reduce the skill of the reconstruction (Franke et al., 2020).

The iCESM1 data used in the DA framework are available at a monthly resolution. Yet, MAIDEN runs on a daily basis. To convert monthly data into daily data, we apply the same methodology to all the particles. We use a reference year of the daily GMF dataset (1950), which we modified to match the mean seasonal cycle over 1950–2000 CE, i.e., a climatological seasonal cycle. Then, for a given particle, the monthly anomalies of the iCESM1 data relative to 1950–2000 CE are added (near-surface air temperature) or multiplied (precipitation) to the reference year of the daily GMF data. Note that the iCESM1 does not include maximum and minimum temperature data. We use the anomalies of iCESM1 near-surface air temperature on both maximum and minimum daily GMF temperature. As for atmospheric CO₂ data, the annual data from Rubino et al. (2019a) (154–1996 CE; publicly available at https://doi.org/10.25919/5bfe29ff807fb, Rubino et al., 2019b, followed from 1996 by the annual data used as forcing in the Goddard Institute for Space Studies – GISS – climate model​​​​​​​, 1850–2011 CE; https://data.giss.nasa.gov/modelforce/ghgases/, last access: 10 January 2021) were linearly interpolated at a daily time step.

The temporal autocorrelation in the simulated Dstem can be important through the storage pool (Gea-Izquierdo et al., 2015; Gennaretti et al., 2017; Guiot et al., 2014) so that the performance of the model in reproducing the observed tree growth of 1 year highly depends on the performance of the previous years. As a consequence, we cannot run the MAIDEN model independently for each year driven by the iCESM1 data as is the case with the linear regression (Sect. 2.4.2). Instead, MAIDEN is sequentially run year by year using the reconstructed climate for the previous years and thus accounting for the effect of past climate history on tree-ring growth. For any particle, we use the CO₂ data corresponding to the assimilated year.

Dstem is not directly comparable with TRW data. The observations and simulation results are thus normalized (Sect. 2.3.1) to produce unitless indexes. The mean and standard deviation of a reference run of MAIDEN with iCESM1 inputs over 1900–2000 CE are used to normalize the annual Dstem output throughout the DA procedure.

2.4.2 Linear regression

A regression-based model was developed for each of the 21 TRW records (Sect. 2.3.1), similarly to Dalaiden et al. (2021) and to the Last Millennium Reanalysis (Tardif et al., 2019). Different climate variables were considered to model the relationship between climate and tree-ring growth. They are related to temperature, precipitation or both (univariate or bivariate model). For temperature, we considered the current year July to the next June temperature, the December–January–February (i.e., austral summer) temperature, and September to March (i.e., the growth period temperature). For precipitation, we considered the July to June precipitation, the July–August precipitation and the July–August–September precipitation. We built the different climate variables by averaging the monthly near-surface air temperature and precipitation from the GMF interpolated on the iCESM1 grid (Sect. 2.2) over the period of months of interest. The climate data were extracted from the continental grid cell closest to each individual TRW record.

Similarly to the MAIDEN calibration, the regression parameters of each statistical PSM were first calibrated over the 1950–2000 CE time period using the ordinary least squares method. For each record, we tested all the possible combinations of climatic variables. As in Dalaiden et al. (2021), the best model was selected for each TRW record using the Bayesian information criterion (BIC) (Schwarz, 1978), which has the advantage of penalizing the complexity (i.e., the number of explanatory variables) of a model.

The regression models were then run independently over the 1901–1949 CE verification time period. Among the 21 sites, only two locations satisfied the same selection criterion as MAIDEN (Sect. 2.4.1), i.e., a significant calibration and verification correlations greater than 0.3 (Table S1). As a consequence, the validity of the regression models out of the calibration interval is not ensured for most of them. Alternatively, we decided to use a less strict selection criterion similar to what is generally done in paleoclimate DA-based reconstruction using TRW proxies with a statistical model (e.g., the Last Millennium Reanalysis; Tardif et al., 2019). The regression parameters are thus calibrated using the same method over the entire period covered by climate data, i.e., 1901–2000 CE, in order to have the largest possible sample and thus the most robust estimate of the relationship between climate and tree growth. We then excluded nine sites (one in New Zealand and eight in South America) for which the explanatory climate variables were not significantly and positively correlated with observed TRW series. Following Cook et al. (2013) and as stated in PAGES 2k Consortium (2017), trees growing better under increasing temperature should be more reliable recorders of past temperature variability because they respond to temperature-limiting conditions. The same assumption can hold for precipitation. The 12 sites included in the DA procedure with the regression model (hereafter, selected sites) can be seen in Fig. S4 (blue circles). A similar criterion based on the calibration only could not be used to select the TRW records for MAIDEN. The complexity of the model and the number of calibrated parameters make it more prone to overfitting problems and biological inconsistency during the calibration (Rezsöhazy et al., 2020, 2021). As a consequence, the selection of the TRW records has to be more restrictive but may also lead to a better stability of the model. For comparison with the observed TRW series, anomalies of the simulated tree-ring indexes are computed relative to 1900–2000 CE.

2.5 Observation error

In the DA framework, the observation error determines the confidence that we can have in the climate information retrieved from the available proxy records. It directly influences the magnitude of the weight given to a particle when compared to an assimilated proxy observation, which in turn has direct consequences for the resulting reanalysis. We generally identify three types of observation error (e.g., Badgeley et al., 2020).

The first type of error is related to the inaccuracy of the PSM in reproducing the proxy dependency on climate as the PSMs cannot fully account for the complexity of the processes by which tree growth responds to climate and nonclimatic drivers. The error of MAIDEN is computed using the standard deviation of the residuals over the verification period (1901–1949 CE; Sect. 2.4.1). The error of the TRW regression-based PSM (Sect. 2.4.2) is computed as in Tardif et al. (2019) using the standard deviation of the residuals over the calibration period (1901–2000 CE). However, computing the error over the calibration period can have serious shortcomings. It does not account for the actual skill of the model to accurately estimate tree growth outside the calibration interval and is biased by possible overfitting problems associated with any statistical model. This could result in an underestimation of the error associated with the statistical PSM and, as a consequence, overconfidence in the PSM performance and finally in the climate reconstruction. The snow accumulation and δ¹⁸O records can be directly compared to the model outputs (Sect. 2.4) so that a PSM is not necessary, and therefore an error related to the performance of the PSM is not needed for these records.

Second, the representativeness error (e.g., Oke and Sakov, 2008) corresponds to the uncertainty related to the unresolved processes at the resolution of the climate model relative to the smaller scale potentially represented by proxy observations (tree-ring width, snow accumulation and δ¹⁸O in our study). For the regression-based TRW PSM, we consider the representativeness error to be implicitly taken into account by computing the error related to the PSM using observed data at the same resolution as the iCESM1 data (Sect. 2.4.2). For MAIDEN, the bias correction applied to the iCESM1 data compared with the high-resolution GMF data (Sect. 2.4.1) used for the calibration of the model is associated with an error, but we assume that it is weaker compared to the error related to the PSM.

Regarding the ice core records (snow accumulation and δ¹⁸O), the representativeness error was computed as in Dalaiden et al. (2021) on the basis of an Antarctic simulation performed with the latest version of the polar-oriented Regional Atmospheric Climate MOdel (RACMO2.3p2, hereafter RACMO) at 27 km horizontal resolution (van Wessem et al., 2018). They computed the standard deviation of the differences between the annual RACMO time series (1979–2016 CE) of snow accumulation averaged for all the ice core sites and the time series averaged over the 500 km grid, following the same approach as in Valler et al. (2020). As for the δ¹⁸O, the same methodology was applied but with the iCESM1 linearly interpolated on the RACMO grid over 1950–2005 CE, as RACMO does not simulate the water isotopes. The error was then corrected to account for the underestimation of spatial variability of snow accumulation and δ¹⁸O in models compared to point measurements as ice cores by multiplying the resulting standard deviation by a factor of 3 (see Dalaiden et al., 2021, for more detail). The last type of error is the instrumental error related to the measurement of the proxy data, which is often negligible compared to the other sources of error, particularly in paleoclimatology (Steiger et al., 2018; Tardif et al., 2019).

3 Experiment design and diagnostics

Different DA experiments are run in this study (Table 1), varying in terms of (1) the records assimilated, (2) the tree-ring PSM included and its associated period of calibration, (3) the PSM error, and (4) the period of the reconstruction.

Table 1 Data assimilation experiments in this study and experimental design: experiment name; types of proxies assimilated; proxy system model (PSM) used to assimilate the proxy data; error associated with the PSM; calibration period of the PSM; the proxy records assimilated. The notation n/a stands for not applicable. See Sect. 2.4.1 and 2.4.2 for details on the validated and selected sites, respectively. The number of tree-ring width records assimilated is in brackets.

A first set of experiments is conducted over the 1900–2000 CE time period and compared with gridded climate datasets in order to evaluate the performance of our DA-based reconstructions regarding the use of statistical (TIC-reg) and process-based (TIC-MDN) dendroclimatic models separately, assimilating both ice core and tree-ring proxies.

A second set of experiments is performed over the 1900–2000 CE time period to evaluate the contribution of tree-ring proxies in the data assimilation framework over the last century. Specifically, we perform an experiment assimilating tree-ring proxies only with the regression (TREES-reg) and MAIDEN (TREES-MDN) models and an experiment assimilating ice core proxies only (IC). These experiments are compared with the TIC-reg and TIC-MDN experiments incorporating both types of proxy data. The impact of the observation error associated with the tree-growth PSMs (Sect. 2.5) on the performance of the reconstruction is also evaluated by reducing the error by a factor of 2 in the experiments constrained by tree-ring proxies only (TREES-reg-05 and TREES-MDN-05).

Lastly, we perform a DA-based reconstruction over the period covered by a sufficient number of TRW proxy records, i.e., 1600–2000 CE, with the regression model (TIC-reg) and the MAIDEN model (TIC-MDN), with both tree-ring and ice core proxies assimilated.

Our analysis concentrates on four regions or subregions of the midlatitudes of the SH: South America, defined between 25 and 56^∘ S and between 55 and 85^∘ W (SA); the southwestern Andes, where the TRW proxy sites are located, defined between 35 and 42.5^∘ S and between 70 and 85^∘ W (SAmSW); New Zealand, defined between 29 and 51^∘ S and between 160 and 179^∘ E (NZ); and Tasmania, defined between 39.5 and 47^∘ S and between 140 and 154^∘ E (Tas). We did not include Antarctica in our study as tree-ring records have not been shown to significantly affect the reconstruction in this region (see Dalaiden et al., 2021, for more details on the reconstructions in Antarctica).

We evaluate our reconstructions against the most recent atmospheric reanalysis from the ECMWF, ERA5 (Hersbach et al., 2020), over the 1979–2000 CE time period and against the GMF dataset (see Sect. 2.4.1 for more details on the dataset) over 1901–2000 CE after linearly interpolating them on the iCESM1 grid (Sect. 2.2). We focus on the near-surface air temperature, the cumulative precipitation and the 500 hPa geopotential height (ERA5 only). In addition to these variables, we also evaluate the Southern Annular Mode (SAM) over 1957–2000 CE. The SAM is the principal mode of atmospheric variability in the extratropics of the SH (Marshall, 2003). It is associated with the north–south movement of the westerlies and often defined as the difference between the normalized monthly zonal mean sea level pressure at 40 and 65^∘ S (Gong and Wang, 1999; Marshall, 2003), as chosen here. The SAM index from Marshall (2003) is used as a reference for our evaluation.

For each of the experiments, the performance of our reconstruction is evaluated with the computation of two statistical indicators: the Pearson correlation coefficient (r; considered to be significant at the 95 % confidence level) and its associated p value as well as the coefficient of efficiency (hereafter CE). The coefficient of efficiency is computed as in Nash and Sutcliffe (1970):

$$\begin{array}{}\text{(1)}& \text{CE}=\mathrm{1}-\sum _{i=\mathrm{1}}^n{\displaystyle \frac{({\mathit{x}}_i-{\mathit{y}}_i{)}^{\mathrm{2}}}{({\mathit{x}}_i-\stackrel{\mathrm{‾}}{\mathit{x}}{)}^{\mathrm{2}}}},\end{array}$$

with n being the number of samples, x the observation vector and y the reconstruction vector. The overbar characterizes a time-averaged vector of values. The average is made over 1901–2000 CE in all experiments, except for the SAM index when compared to the observations from Marshall (2003), as they only cover the 1957–2000 CE period. While the Pearson correlation coefficient gives an estimation of the linear fit between the reconstructed and observed values, the CE compares the variance of the residuals with the variance of the predictand. A CE equal to 1 is associated with perfect predictions, a CE equal to 0 with a residual variance identical to predicting the observed mean, and a negative CE with larger residual variance than when predicting the mean. For regional near-surface air temperature and precipitation, the metrics are computed between the mean of the reconstruction and the mean of the ERA5 or GMF data over the continental mass that lies in the region delimited by the latitudes and longitudes above.

It is important to note at this stage that the numbers of proxy records, but also the tree-ring proxy records themselves, that are assimilated with the statistical or MAIDEN PSMs are different. The regression-based PSM experiment includes twice the number of tree-ring proxy records compared to the one using MAIDEN. This results from the different performance of the PSMs when applied to the same subset of SH tree-ring sites (Sect. 2.4.1 and 2.4.2). When working with actual proxy data, such a situation is difficult to avoid but we must be aware that this may have consequences for the skill of the resulting DA-based reconstructions.

4 Results and discussion

4.1 Comparison of the reconstructions using regression- and process-based dendroclimatic proxy system models over the last century

The spatial distribution of the correlation coefficients between different ERA5 and reconstructed climate fields over the recent period (1979–2000 CE) indicates generally good agreement of the reconstruction using MAIDEN as a PSM south of 30^∘ S (TIC-MDN; Fig. 2a). When compared with the DA experiment including the regression model to assimilate 12 TRW records (TIC-reg; Fig. 2b), the hemispheric 500 hPa geopotential height, near-surface air temperature and precipitation fields are reconstructed with a similar performance in both cases.

Figure 2 Pearson correlation coefficient between the reconstructed and ERA5 geopotential height at 500 hPa, near-surface air temperature, and cumulative precipitation over 1979–2000 CE for the TIC-MDN (a) and TIC-reg (b) experiments (see Sect. 3 and Table 1 for details on the experiments). The green dots indicate the localization of the assimilated tree-ring width records. The black dots indicate a significant correlation coefficient at the 95 % confidence level.

At the regional scale and over the recent period (1979–2000 CE), the correlation coefficients of the ERA5 and GMF near-surface air temperature with the reconstructed near-surface air temperature in the four regions (Table 2) are also generally similar for both TIC-reg and TIC-MDN. Except in Tasmania (Tas), where the TIC-MDN near-surface air temperature is highly and significantly correlated with both ERA5 and GMF, the correlation coefficients are low and not statistically significant for both experiments. This differs from the good general pattern previously noted (Fig. 2). Locally, the correlations can indeed be lower (see Figs. S5–S10 for the regional pattern of the correlations). It is also important to highlight that these correlations have been computed between the integrated fields and are not the average of the correlations represented on the map, which may result in slight differences. Regarding precipitation (Table 2), the correlation coefficients are particularly high and significant in South America (SA) and in the southwestern Andes (SAmSW) for TIC-reg and TIC-MDN.

Table 2 Pearson correlation coefficient between the reconstructed as well as GMF and ERA5 near-surface air temperature and precipitation over 1979–2000 CE for different DA experiments with both ice core and tree-ring proxies (Table 1) and for different regions (Sect. 3). SA stands for South America, SAmSW for the southwestern border of South America in the Andes where the TRW proxy sites are located, NZ for New Zealand and Tas for Tasmania. Asterisks stand for significant correlation coefficients at the 95 % confidence level.

Figure 3 Reconstructed (in color) and GMF (in black) near-surface air temperature anomalies (in ^∘C; relative to 1901–2000 CE) over 1901–2000 CE for different DA experiments with both ice core and tree-ring proxies (Table 1) in different regions (Sect. 3). The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

Figure 4 Reconstructed (in color) and GMF (in black) precipitation anomalies (in mm yr⁻¹; relative to 1901–2000 CE) over 1901–2000 CE for different DA experiments with both ice core and tree-ring proxies (Table 1) in different regions (Sect. 3). The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

Over the whole century (1901–2000 CE), both DA experiments have similar skills (ranging from −0.12 to 0.60 for the correlation coefficients and from −0.12 to 0.36 for the CE for the four regions; Fig. 3; Table 3) in reconstructing the near-surface air temperature when compared to the GMF climate dataset. Except in the southwestern Andes region (SAmSW in Table 3), the correlation coefficients and CE between reconstructed and GMF temperature are close for both experiments and particularly high in New Zealand (NZ). In Tasmania, no TRW records have been assimilated with MAIDEN, which may explain the lower CE (Table 3) for the TIC-MDN experiment, while the correlation coefficients are similar between the two experiments. As for precipitation (Table 3 and Fig. 4), both TIC-MDN and TIC-reg are more skillful in the vicinity of the South American TRW proxies (SAmSW), with significant correlation coefficients of 0.35 and 0.37, respectively, and positive CE.

Table 3 Pearson correlation coefficient and coefficient of efficiency (CE) between the reconstructed and GMF near-surface air temperature and precipitation over 1901–2000 CE for different DA experiments with both ice core and tree-ring proxies (Table 1) and for different regions (Sect. 3). SA stands for South America, SAmSW for the southwestern border of South America in the Andes where the TRW proxy sites are located, NZ for New Zealand and Tas for Tasmania. Asterisks stand for significant correlation coefficients at the 95 % confidence level.

Overall, the variance of the ensemble mean reconstruction is far smaller than the observed variance for both near-surface temperature and precipitation (Figs. 3 and 4), but observations are generally well in the range of the reconstruction (except for extreme years). Since this pattern is apparent in both experiments, i.e., when using a regression and when using a process-based model as a PSM, this is likely not predominantly related to the intrinsic characteristics of each PSM. Our DA experiments only include 6 and 12 TRW proxy sites in the midlatitudes of the SH for TIC-MDN and TIC-reg, respectively. Comparatively, in the high latitudes, the ice core records are far more numerous (Sect. 2.3.2). However, the ice core records can only have an influence on the reconstructed variance over a limited area. As a consequence, the number of assimilated TRW proxies in the midlatitudes is probably not sufficient to effectively constrain the climate model simulations during the DA procedure. In the DA framework, when no proxy is assimilated, the posterior as defined in Sect. 2.1 is equivalent to the prior, in which all the particles have the same weight. The posterior and the resulting ensemble mean are thus the same for each year of the reconstruction. As a consequence, the reconstruction has, by construction, a variance equal to zero. For a small number of proxies providing a weak constraint, many particles receive a similar weight, also resulting in small variance. By contrast, the uncertainty of the reconstruction (as measured by the standard deviation of the posterior or the ensemble spread; Sect. 2.1) remains large. Additionally, both the MAIDEN and regression models work on the basis of normalized observed and simulated tree-ring growth indexes. Consequently, in the DA procedure, the error in the simulated variance of the climate signal is likely not properly assessed when using such TRW PSMs to compare the climate model outputs with the proxy observations. It is in practice possible to transform tree-ring data in order to directly compare them with the output of the MAIDEN model without normalization (Gennaretti et al., 2018). This implies combining both tree-ring width and density measurements that are less widely available (PAGES 2k Consortium, 2017). This could improve the model–data comparison in the DA procedure and, consequently, the resulting reconstruction. This is not the case for the regression model.

The DA procedure seeks to identify the climate conditions (i.e., the inputs of the TRW PSM; near-surface air temperature and precipitation) that gave rise to the observed TRW through the use of a PSM. As a consequence, we must also establish whether the reconstruction of near-surface air temperature and precipitation translates into the local ability of the MAIDEN model to reproduce tree growth at the six assimilated TRW sites over the 20th century. First, we have retrieved the monthly climate fields from the ensemble mean of the TIC-MDN reconstruction. Then, the monthly climate inputs have been converted to daily inputs as in Sect. 2.4.1. For consistency, the tree-growth time series of the MAIDEN model have been normalized with the mean and the standard deviation of the reference run used in the DA procedure (Sect. 2.4.1). If we only look at the correlation coefficients between the observed and simulated tree-growth indexes (Fig. S11), MAIDEN has been successfully run over 1900–2000 CE with the TIC-MDN reanalysis of near-surface air temperature and precipitation as inputs at most assimilated sites (with a mean correlation of 0.59). The TRW sites displaying lower correlation coefficients between the observed and simulated tree-ring indexes are the ones associated with higher PSM errors (Sect. 2.5), i.e., Aus_002 and SAm_024 (i.e., the name of the tree-ring sites, as in the PAGES 2k database; Aus_* stands for a Tasmanian or New Zealand site, while SAm_* stands for a South American site). This result is logical as lower confidence is given to the latter in the particle filter (Sect. 2.1) so that the climate reconstructed at the local scale is less influenced by those records. However, the simulated variance is strongly underestimated compared to the observations, as seen in Fig. S11. This underestimation of the variance of tree growth comes directly from the underestimation of the variance in the reconstructed near-surface air temperature and precipitation, as highlighted above. The same pattern of small simulated variance is observed with the tree-ring indexes simulated with the linear regression.

Finally, the reconstructed SAM index is well and significantly correlated with the SAM index from Marshall (2003) for TIC-MDN (0.46) and TIC-reg (0.50) over the 1957–2000 CE period (Fig. 5). However, the TIC-reg experiment slightly better captures the significant observed trend of 0.28 standard deviation per decade over 1957–2000 CE, with a significant reconstructed trend of 0.34 and nonsignificant reconstructed trend of 0.20 standard deviation per decade for TIC-reg and TIC-MDN, respectively. However, both reconstructions display negative values for the CE, which indicates that they are unable to accurately reproduce the observed amplitude of the SAM variations.

Figure 5 Z-scored (relative to 1957–2000 CE) Southern Annular Mode (SAM) index of Marshall (2003) (in black) and reconstructed SAM index (in color) over 1957–2000 CE for different DA experiments with both ice core and tree-ring proxies (Table 1); Pearson correlation coefficient (r), coefficient of efficiency (CE), and decadal trend of the reconstructed SAM index (trend). Asterisks stand for a significant correlation coefficient or trend at the 95 % confidence level. The observed trend of the SAM index is 0.28^∗ standard deviation per decade over 1957–2000 CE. The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

Overall, the DA experiment using MAIDEN as a PSM for assimilating the TRW records does not outperform the experiment using the regression model. However, its skill is always as good as the TIC-reg experiment with potential better stability of MAIDEN through time, since it has been independently validated at the assimilated sites. It is also worth noting that the DA experiments with MAIDEN or the regression model have some characteristics that can influence their respective skills. Firstly, the DA experiment with the regression model includes twice the number of TRW records compared to the experiment with MAIDEN (12 and 6, respectively). The TRW records can also be located at different places, which may result in different local skills of the DA experiments (e.g., in Tasmania). Secondly, the regression model was calibrated over the whole century (1901–2000 CE), while MAIDEN was calibrated over the second half of the century. This may confer a comparative advantage to the regression-based DA experiment, particularly when evaluated against the GMF climate dataset over 1901–2000 CE.

4.2 Contribution of tree-ring proxies to the skill of the reconstructions over the last century

The experiments discussed in Sect. 4.1 include different types of proxies. For evaluating the specific contribution of tree-ring proxies using a regression- or process-based PSM, we focus on the regions and climate variables for which tree-ring data could likely provide a strong constraint (Sect. 4.1). This evaluation is essential to ensure that the few tree-ring proxies assimilated – compared to the ice core data – are efficiently taken into account in a realistic framework, including proxies in different regions using both MAIDEN and the regression-based PSM, and whether differences arise between the experiments with these PSMs. Complementary information can be obtained by analyzing experiments in which only a fraction of the proxies is assimilated.

Figure 6 Reconstructed (in color) and GMF (in black) temperature (in ^∘C; relative to 1901–2000 CE) or precipitation anomalies (in mm yr⁻¹; relative to 1901–2000 CE) over 1901–2000 CE for different DA experiments (Table 1, with ice core proxies only, with tree-ring proxies only or with both) in different regions (Sect. 3) and proxy sites (Fig. 1). Aus refers to a site in New Zealand and SAm in South America. The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

In South America (SA; Table 4 and Fig. 6), the correlations with observations are low for the experiments with tree-ring proxies alone, whatever the tree-ring PSM used (TREES-reg and TREES-MDN). Ice core proxies alone (IC) perform better than the tree-ring proxies alone (TREES-reg and TREES-MDN) or when combined with the tree-ring proxies (TIC-reg and TIC-MDN) in terms of correlation coefficients, but not in terms of reconstructed variance (standard deviation for TIC-reg, TIC-MDN, TREES-reg or TREES-MDN higher than IC). However, as previously noted, it is clear that, even though the reconstructions with tree-ring proxies display higher variance, none of the DA experiments accurately reconstruct the observed variance of precipitation (Fig. 6).

Table 4 Pearson correlation coefficient between the reconstructed and GMF near-surface air temperature and precipitation over 1901–2000 CE for different DA experiments (Table 1, with ice core proxies only, with tree-ring proxies only or with both) in different regions (Sect. 3) and proxy sites (Fig. 1). Standard deviation of the reconstruction in degrees Celsius (^∘C) (near-surface air temperature) or in millimeters per year (mm yr⁻¹) (precipitation). SA stands for South America, SAmSW for the southwestern border of South America in the Andes where the TRW proxy sites are located and NZ for New Zealand. Aus refers to a site in New Zealand and SAm in South America. Asterisks stand for significant correlation coefficients at the 95 % confidence level.

At the more local scale, in the southwestern Andes (SAmSW), the tree-ring proxies alone provide a better reconstruction skill than for South America (TREES-reg and TREES-MDN; Table 4). The correlation coefficients are, however, slightly lower than in the experiments using both tree-rings and ice core proxies or ice core proxies alone, in particular for TREES-MDN. Nevertheless, the TREES-reg and TREES-MDN reconstructions have higher reconstructed variance compared to IC, in particular for TREES-reg (Fig. 6). The combination of both tree-ring and ice core proxies (TIC-reg and TIC-MDN) provides the best reconstruction in the southwestern Andes in terms of correlation coefficients and reconstructed variance.

Only one proxy record is used in both TREES-reg and TREES-MDN in South America (i.e., SAm_11). For this site (i.e., the grid cell closest to the proxy site), the correlation coefficients are similar for all the experiments so that the different proxy records all bring some information, but when tree-ring records are included, an improvement of the reconstructed variance is obtained. More generally, at all spatial scales in South America the performance of TREES-reg in terms of correlation coefficients and reconstructed variance is slightly better than TREES-MDN.

In New Zealand, the correlations with observations are similar in all the experiments (NZ; Table 4). In other words, tree-ring proxies only are able to provide some skill at the scale of the New Zealand, even without assimilating ice core proxies and with slightly higher reconstructed variance (especially for TREES-reg). The contrast between the reconstructed variance in IC and TREES-reg or TREES-MDN experiments is stronger at the scale of the proxy site (Aus_005 and Aus_030, the tree-ring proxy sites assimilated with both MAIDEN and the regression model). At the Aus_030 proxy site the correlation coefficients are even higher for the experiments assimilating tree-ring proxies only compared to ice core proxies only. However, the reconstructed variance is still underestimated compared to observations (Fig. 6).

Figure 7 Reconstructed (in color) and GMF (in black) temperature (in ^∘C; relative to 1901–2000 CE) or precipitation anomalies (in mm yr⁻¹; relative to 1901–2000 CE) over 1901–2000 CE for different DA experiments (Table 1, with tree-ring proxies only, with the full observation error or half of the observation error; see Sect. 3 for more details) in different regions (Sect. 3) and proxy sites (Fig. 1). Aus refers to a site in New Zealand and SAm in South America. The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

The underestimation of the variance compared to the observations could arise from confidence that is too low given to the assimilated tree-ring records in the data assimilation procedure, which is directly related to the error associated with the PSMs (Sect. 2.5). In order to test this hypothesis, we have reduced this error by 2. This leads to an increase in the reconstructed variance at all spatial scales in both South America (precipitation) and New Zealand (near-surface air temperature; Table 5 and Fig. 7), in particular for TREES-MDN in New Zealand. This does not translate, however, into higher correlation coefficients, since they virtually remain unchanged (Table 5). In an ideal experiment in which the relationship between climate and tree-ring growth would be perfectly reproduced, we could thus expect the observed variance to be more skillfully reconstructed by the ensemble mean reconstruction.

Table 5 Pearson correlation coefficient between the reconstructed and GMF near-surface air temperature and precipitation over 1901–2000 CE for different DA experiments (Table 1, with tree-ring proxies only, with the full observation error or half of the observation error; see Sect. 3 for more details) in different regions (Sect. 3) and proxy sites (Fig. 1). Standard deviation of the reconstruction in degrees Celsius (^∘C) (near-surface air temperature) or in millimeters per year (mm yr⁻¹) (precipitation). SA stands for South America, SAmSW for the southwestern border of South America in the Andes where the TRW proxy sites are located and NZ for New Zealand. Aus refers to a site in New Zealand and SAm in South America. Asterisks stand for significant correlation coefficients at the 95 % confidence level.

More generally, this sensitivity analysis shows that tree-ring proxies have a positive effect on the reconstructed midlatitude climate changes in our framework with both the regression- and process-based dendroclimatic PSM, especially at the regional to local scale. The fact that the reconstructed variance can increase at the more local scale when assimilating tree-ring proxies alone indicates that the number of tree-ring proxies assimilated in our experiments likely has an impact on the reconstructed variance. If more proxies were available, we could expect that this local increase in the variance would lead to an increase at larger spatial scale. The definition of the error related to the PSMs likely plays a role as well. In particular, the computation of the error for MAIDEN is stricter than for the regression model (i.e., computed on the verification period and on the calibration period, respectively; Sect. 2.5) and the number of proxy records assimilated with MAIDEN is lower as well, which may advantage the regression-based PSM.

4.3 Data-assimilation-based reconstruction of the Southern Hemisphere climate over the last 400 years

In this section, we compare the DA-based reconstructions with the process- or regression-based PSM over the past 4 centuries to expand the conclusions obtained from the comparison of the two methodologies over the 20th century to longer timescales. If we first look at the regional reconstructions of near-surface air temperature (Fig. 8) and precipitation (Fig. 9) for the four regions, we observe that the ensemble spread is generally large in contrast to the interannual variability of the ensemble mean reconstructions, which is generally small. This could indicate that we overestimate the error associated with the PSM and therefore underestimate the confidence of TRW data in the DA procedure. As highlighted before, this is at least partly due to the small number of assimilated TRW proxies in the midlatitudes that prevents the prior distribution given by the climate model from being effectively constrained by the available observations. This pattern is stronger when the spatial scale is large (e.g., South America) and weaker when the spatial scale is small (e.g., Tasmania) and closer to the proxy sites. This is particularly true for MAIDEN, for which we assimilated half the number of TRW proxy records compared with the regression model.

Figure 8 Reconstructed TIC-reg (in blue) and TIC-MDN (in orange) near-surface air temperature anomalies (in ^∘C; relative to 1901–2000 CE) over 1600–2000 CE (Table 1) in different regions (Sect. 3). The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

Figure 9 Reconstructed TIC-reg (in blue) and TIC-MDN (in orange) precipitation anomalies (in mm yr⁻¹; relative to 1901–2000 CE) over 1600–2000 CE (Table 1) in different regions (Sect. 3). The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

In New Zealand, both TIC-reg and TIC-MDN display a smaller spread and larger variability of the reconstructed near-surface air temperature compared to the other regions (Fig. 8), which suggests that the DA has been more effective in updating the climate model simulations with the available proxy records. The reconstructed near-surface air temperature in New Zealand was also well evaluated against GMF observations over 1901–2000 CE in Sect. 4.1 (Table 3). The low-frequency variability of the two reconstructions over the last 400 years is relatively similar overall, with a very small trend before the industrial era (until around 1850–1900 CE) followed by an increase in temperature. TIC-reg and TIC-MDN display a significant trend of 0.047 and 0.023 ^∘C per decade during the 20th century (against 0.062 ^∘C per decade in the GMF observations) as well as a significant trend of 0.001 and nonsignificant trend of 0 ^∘C per decade before 1900 CE, respectively. The correlation coefficient between the two reconstructions is 0.63 (significant at the 95 % confidence level), which indicates that the trend and timing of the reconstructions are close to each other. However, TIC-MDN displays higher near-surface air temperature anomalies before 1850 CE than TIC-reg (−0.05 ^∘C on average for TIC-MDN against −0.14 ^∘C for TIC-reg relative to 1901–2000 CE). In Tasmania, there are no assimilated sites for MAIDEN (Figs. S3 and S7) and four assimilated sites for the regression over a very small area (Figs. S4 and S10), which may explain the difference between the reconstructed near-surface air temperature of TIC-reg and TIC-MDN (Fig. 8). This is consistent with the results obtained over the 20th century. Finally, we observe a difference in the reconstructed low-frequency variability between TIC-reg and TIC-MDN. Compared with TIC-reg, TIC-MDN tends to reconstruct a weaker temperature difference between the preindustrial and industrial era, particularly in South America, the southwestern Andes and Tasmania.

For precipitation (Fig. 9), the southwestern Andes is the region for which the ensemble spread is the smallest and the interannual variability the largest for both TIC-reg and TIC-MDN. Both experiments are also skillful for reconstructing precipitation in this region when compared with state-of-the-art climate datasets over the last century (Table 3; Sect. 4.1). The correlation coefficient between the two reconstructions is 0.58 and is statistically significant so that the trend and timing of the reconstructions are close to each other. Both reconstructions display a very small trend before 1900 CE (a significant trend of −0.586 and a nonsignificant trend of −0.042 mm yr⁻¹ per decade for TIC-reg and TIC-MDN, respectively) and a negative trend after 1900 CE (a significant trend of −12.79 and −8.89 mm yr⁻¹ per decade for TIC-reg and TIC-MDN, respectively, against a nonsignificant trend of −1.26 mm yr⁻¹ per decade in the GMF observations). However, TIC-MDN displays higher anomalies before 1900 CE than TIC-reg (28.23 mm yr⁻¹ on average for TIC-MDN against −5.70 mm yr⁻¹ for TIC-reg relative to 1901–2000 CE).

At the regional scale, we have highlighted from the reconstructed near-surface air temperature and precipitation that the experiments with the statistical or process-based PSM show different amplitudes of climate change over the past 400 years, even if close in terms of correlation (i.e., the trend and timing of change). It is not possible, of course, to determine which one is the most accurate. This difference in the reconstructions could be related to the number of proxies included in each of the DA experiments and their associated error (Sect. 2.5). This also points out that the commonly used linear regression may also not be fully able to always properly reflect the full spectrum of the climate variability over the past centuries. Alternatively, MAIDEN may underestimate low-frequency trends, for instance because it overestimates the effect of CO₂ concentration on tree growth (Rezsöhazy et al., 2021), which leads to an underestimation of the near-surface temperature low-frequency variability. It may also be related to the calibration of MAIDEN that may not be able to properly capture the low-frequency signal in tree-ring series.

Finally, we observe an increase in the reconstructed variability in the precipitation of the southwestern Andes (Fig. 9), particularly for TIC-MDN, from around 1800–1850 CE, i.e., when most of the ice core records start to become available (Fig. S12). Before 1830 CE, the standard deviation of the reconstruction is 30.7 and 21.0 mm yr⁻¹; after 1830 CE it is 56.3 and 44.2 mm yr⁻¹ for TIC-reg and TIC-MDN, respectively.

This increase in interannual variability is even more striking in the reconstructed SAM index for both reconstructions (Fig. 10). Before around 1830 CE particularly, too few TRW proxy sites have been assimilated in the midlatitudes to be able to reconstruct the large-scale variability adequately (Fig. S12). The standard deviation of the z-scored reconstructed SAM index is 0.52 and 0.51 before 1830 CE and 1.05 and 1.11 after 1830 CE for TIC-reg and TIC-MDN, respectively. However, the SAM reconstructions are still highly and significantly correlated with each other (r = 0.91), and they look essentially the same over the whole 1600–2000 CE period (Fig. 10). Before 1830 CE, when the number of ice core records is more limited (Fig. S12), the reconstructions are less correlated (r=0.77 and r=0.95 after 1830 CE).

Figure 10 Z-scored (relative to 1901–2000 CE) TIC-reg (in blue) and TIC-MDN (in orange) reconstructed SAM index over 1600–2000 CE (Table 1). The uncertainty (i.e., the ensemble spread; Sect. 2.1) associated with each reconstruction is represented by the shaded areas.

Lastly, over the last 400 years, the MAIDEN model generally reproduces the observed tree growth well in terms of correlation coefficients at the six assimilated sites (with a mean correlation of 0.53; Fig. S13). However, as previously highlighted in Sect. 4.1, the simulated variance is largely underestimated compared to the observations, which is not surprising regarding the reconstructions of near-surface air temperature and precipitation of the past 400 years.

5 Conclusions

In this study, we have included for the first time a dendroclimatic process-based model into a data assimilation procedure to reconstruct the past climate variability of the Southern Hemisphere. All the results from our DA experiments have been compared with the commonly used regression-based PSM. First, we have shown that including the ecophysiological model MAIDEN as a proxy system model in a data assimilation procedure is technically possible. The assimilation of actual TRW records with a tree-growth mechanistic model has never been published before. Knowing that it is at least achievable is promising.

The general pattern of near-surface air temperature, precipitation and the atmospheric circulation of the middle to high latitudes of the SH is reproduced well by the DA experiment with MAIDEN over the recent period. Over the last century, the DA experiment with MAIDEN is also able to produce reconstructions that are well correlated with different climate fields, such as near-surface air temperature, precipitation and atmospheric circulation in some regions of the Southern Hemisphere, with only a few TRW proxy records assimilated. The DA-based reconstructions with the regression have a similar skill as with MAIDEN.

Over the instrumental era, the skill of MAIDEN as a PSM in the present DA framework is not better but it is at least as good as the regression-based PSM. Given the advantages that a process-based PSM can still bring beyond the instrumental period in terms of nonlinearity of the relationship between climate and tree growth, but also in terms of climate variables included as inputs of the model (e.g., atmospheric CO₂ concentration), having skill comparable to a regression-based PSM is encouraging. There are still opportunities for improvement with the MAIDEN model, while for the regression, the perspectives are more limited. For instance, we could improve the model–data comparison by using an observed tree-growth index that is closer to the variable simulated by MAIDEN, i.e., forest carbon accumulation, by combining density and ring width measurements. Different shortcomings related to the calibration of the MAIDEN model have also been addressed in Rezsöhazy et al. (2020, 2021). If solved, this may likely help to increase the number of proxy sites where MAIDEN could be validated. It could also improve the performance of MAIDEN in simulating the relationship between climate and tree growth and therefore enhance the skill of the DA-based reconstruction using MAIDEN as a PSM. In particular, the calibration procedure could be more informative in the future, for instance in the definition of the prior ranges of the calibrated parameters that could be set by species or biome. Different data sources, such as satellite data and the associated vegetation indexes, could also be used to better inform the different biological processes in the model (e.g., related to photosynthesis or phenology). Moreover, in this study, the criteria for including a TRW proxy record with the regression had to be less strict than for MAIDEN. MAIDEN was indeed validated at all the assimilated sites, while the regression was not. This may provide better stability of the MAIDEN model.

It appears clearly that our DA framework does not provide reconstructions that properly reproduce the observed variance of near-surface air temperature and precipitation with both TRW PSMs, particularly at larger spatial scale (e.g., South America), likely due to the low number of TRW proxy records assimilated in the midlatitudes. The Southern Hemisphere indeed counts few TRW proxy records (PAGES 2k Consortium, 2017) but offers an interesting framework as few reconstructions of the Southern Hemisphere climate exist to date (Neukom et al., 2011, 2014). The underestimation of variance is a fundamental issue with the particle filter, but also potentially with other data assimilation methods (e.g., Tardif et al., 2019). In the ideal case, individual members of ensemble climate simulations have realistic temporal variability on all spatial scales. Any averaging of random, non-forced variability will reduce the variability. In other words, any ensemble mean will always have variability that is too low, regardless of how it is constructed. However, the extent of the reduction depends on whether members with more similar or more different variability are used for calculating the weighted ensemble mean. In contrast to unweighted ensemble means, the particle filter will give high weights to ensemble members that match the observations. If the empirical information strongly constrains the particle selection, the particles will be more similar than in a less constrained case, and there will be less underestimation of variance. In addition, the similarity measure is determined at the locations of the proxies, but the weights given to each particle are independent of location. Particles that have similar states at the proxy locations do not necessarily have similar states at other locations, and the reduction of variance in unconstrained locations is therefore likely to be larger than in constrained locations.

In our reconstructions, this pattern of variance becomes even more striking before 1830 CE, for which we observe smaller interannual variability of the reconstructed near-surface air temperature and precipitation in some regions of the midlatitudes as well as of the reconstructed SAM index compared to the last 150 years. However, the timing of the climate changes is well constrained in some localized areas of the midlatitudes, e.g., the southwestern Andes, New Zealand and Tasmania, as shown by the significant and relatively high correlations with the observations. As a consequence, the MAIDEN model simulates tree-ring growth series that are generally well correlated with the observations over the last 400 years but which display small variance compared to the observations.

Despite this, at the regional scale, we have highlighted a difference between the two DA experiments in the reconstructed amplitude of the climate change in the preindustrial era relative to the 20th century. Using MAIDEN as a PSM for TRW records tends to reconstruct an average warmer and wetter climate than the regression-based PSM. This questions the almost exclusive use of linear statistical PSMs for tree-ring proxies in data assimilation to reconstruct the past climate changes. However, this would need an in-depth investigation with more TRW proxy records in order to be confirmed, as well as additional analysis to ensure the validity and robustness of MAIDEN and its calibration in our framework, with generally more sophisticated tree-growth PSMs.

For now, our framework only offers limited skill for the reconstructions, and we are far from providing a complete hemispheric reconstruction of atmospheric fields in the Southern Hemisphere. However, our results are encouraging for MAIDEN and more generally for the use of dendroclimatic process-based models in DA. In the future, mechanistic models like MAIDEN will likely become an important tool of the DA framework that would help overcome the potential bias associated with regression models. However, this would require continuing to improve those models, including their calibration, so that they could be applied more efficiently and more skillfully at large spatial scale. Calibrating a model like MAIDEN at large spatial scale remains a challenge, and incorporating it into a data assimilation procedure is still more technically time-consuming than the traditional statistical approach. Finally, more investigation of the performance of the MAIDEN model as a PSM in the DA procedure compared with the commonly used regression model would also be needed using more TRW proxy records in order to fully document the added value of using a dendroclimatic process-based model in the data assimilation framework.