A method for generating coherent spatially explicit maps of seasonal palaeoclimates from site-based reconstructions

TL;DR

This paper introduces a 3D-Variational method to generate spatially continuous maps of seasonal palaeoclimates from site-based data, effectively integrating reconstructions with prior climate estimates.

Contribution

The novel approach combines site reconstructions and prior data using a 3D-Variational technique, allowing controlled smoothing and producing unbiased, continuous palaeoclimate maps.

Findings

Successfully applied to Last Glacial Maximum data for southern Europe.

Demonstrated how to select appropriate smoothing length scales.

Produced spatially explicit climate anomaly maps from pollen data.

Abstract

We describe a method for reconstructing spatially explicit maps of seasonal palaeoclimate variables from site-based reconstructions. Using a 3D-Variational technique, the method finds the best statistically unbiased, and spatially continuous, estimate of the palaeoclimate anomalies through combining the site-based reconstructions and a prior estimate of the palaeoclimate state. By assuming a set of correlations in the error of the prior, the resulting climate is smoothed both from month to month and from grid cell to grid cell. The amount of smoothing can be controlled through the choice of two length-scale values. The method is applied to a set of reconstructions of the climate of the Last Glacial Maximum (ca. 21,000 years ago, yr BP) for southern Europe derived from pollen data with a prior derived from results from the third phase of the Palaeoclimate Intercomparison Project (PMIP3). We demonstrate how to choose suitable values for the smoothing length scales from the datasets used in the reconstruction.

Equations

$$\mathbf{y}={\left({\mathbf{y}}_1^T|{\mathbf{y}}_2^T|\cdots |{\mathbf{y}}_N^T\right)}^T$$

$$\mathbf{x}={\left({\mathbf{x}}_1^T|{\mathbf{x}}_2^T|\cdots |{\mathbf{x}}_M^T\right)}^T.$$

$$\mathbf{T}={\left(T_1{T}_2\dots T_{12}\right)}^T$$

$$\mathbf{h}(\mathbf{x})=\mathbf{y}.$$

$${\mathbf{x}}_a=\underset{\mathbf{x}}{\min }J(\mathbf{x}).$$

$$J(\mathbf{x})=\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_b{)}^T{\mathbf{B}}^{-1}(\mathbf{x}-{\mathbf{x}}_b)+\frac{1}{2}(\mathbf{y}-\mathbf{h}(\mathbf{x}){)}^T{\mathbf{R}}^{-1}(\mathbf{y}-\mathbf{h}(\mathbf{x})).$$

$$\mathbf{R}=diag({\mathbf{v}}_y).$$

$$\mathbf{B}=\mathbf{\Sigma }\mathbf{C}\mathbf{\Sigma }$$

$$\mathbf{\Sigma }=diag({\mathbf{v}}_b^{\frac{1}{2}}),$$

$$c_L({\theta }_{ij})=\left(\frac{a}{L}\sin \left(\frac{{\theta }_{ij}}{2}\right)\right){\mathcal{K}}_1\left(\frac{a}{L}\sin \left(\frac{{\theta }_{ij}}{2}\right)\right)$$

$$({\mathbf{C}}_{L_s}{)}_{ij}=c_{L_s}({\theta }_{ij}).$$

$${\bf C}_{L_{t}}=\left(\begin{array}[]{c | ccc}1&0&\dots&0\\ \hline\cr 0&&&\\ \vdots&&\{c_{L_{t}}(\theta_{ij}\}_{ij}&\\ 0&&&\\ \end{array}\right).$$

$$\mathbf{C}={\mathbf{C}}_{L_s}\otimes {\mathbf{C}}_{L_t}$$

$$\kappa (\mathbf{M})=\frac{{\lambda }_{\text{max}(\mathbf{M})}}{{\lambda }_{\text{min}(\mathbf{M})}}$$

$${\mathbf{B}}^{\frac{1}{2}}\mathbf{w}=\mathbf{x}-{\mathbf{x}}_b$$

$$\mathbf{B}={\mathbf{B}}^{\frac{1}{2}}{\mathbf{B}}^{\frac{1}{2}}.$$

$${\mathbf{w}}_a=mi{n}_{\mathbf{w}}J(\mathbf{w}).$$

$$J(\mathbf{w})=\frac{1}{2}{\mathbf{w}}^T\mathbf{w}+\frac{1}{2}(\mathbf{y}-\mathbf{h}({\mathbf{x}}_b+{\mathbf{B}}^{\frac{1}{2}}\mathbf{w}){)}^T{\mathbf{R}}^{-1}(\mathbf{y}-\mathbf{h}({\mathbf{x}}_b+{\mathbf{B}}^{\frac{1}{2}}\mathbf{w})).$$

$$\nabla J(\mathbf{w})=\mathbf{w}-{\mathbf{B}}^{\frac{1}{2}}{\mathbf{H}}_{\mathbf{x}}^T{\mathbf{R}}^{-1}\left(\mathbf{y}-{\mathbf{h}}_u({\mathbf{x}}_b+{\mathbf{B}}^{\frac{1}{2}}\mathbf{w})\right)$$

$$\mathbf{A}=\left(\mathbf{I}-\mathbf{K}{\mathbf{H}}_{{\mathbf{x}}_b}\right)\mathbf{B}.$$

$$\mathbf{K}=\mathbf{B}{\mathbf{H}}_{{\mathbf{x}}_b}^T{\left({\mathbf{H}}_{{\mathbf{x}}_b}\mathbf{B}{\mathbf{H}}_{{\mathbf{x}}_b}^T+\mathbf{R}\right)}^{-1}.$$

$${\bf\hat{h}}({\bf x}_{i})={\bf\hat{h}}\left(\begin{array}[]{c}\bar{P}\\ \bar{{\bf T}}\\ \end{array}\right)=\left(\begin{array}[]{c}\mu({\bf x}_{i})\\ \bar{P}\\ mean(\bar{{\bf T}})\\ max(\bar{{\bf T}})\\ min(\bar{{\bf T}})\\ G(\bar{{\bf T}})\end{array}\right).$$

$$\mu ({\mathbf{x}}_i)=1+m({\mathbf{x}}_i)-{\left(1+m({\mathbf{x}}_i{)}^{\omega }\right)}^{\frac{1}{\omega }}$$

$$m({\bf x}_{i})=P\lambda\left[\sum_{k}^{12}{l}_{k}\frac{R(T_{k},S_{k})\frac{\partial e_{s}}{\partial T}\big{|}_{T_{k}}}{\frac{\partial e_{s}}{\partial T}\big{|}_{T_{k}}+\gamma}\right]^{-1}$$

$$\frac{\partial e_s}{\partial T}=\frac{10.5485}{{\left(237.3+T\right)}^2}exp\left(\frac{17.27T}{237.3+T}\right),$$

$$G(\bar{\mathbf{T}})=\frac{1}{N_y}\sum _k^{12}\{\begin{array}{ll}{l}_k\left({\bar{T}}_k-\frac{5}{T_s}\right)& {\bar{T}}_k>\frac{5}{T_s}\\ 0& \text{else}\end{array},$$

$$\mathbf{h}(\mathbf{x})=\left(\hat{\mathbf{h}}({\mathbf{x}}_1{)}^T|\hat{\mathbf{h}}({\mathbf{x}}_2{)}^T|\cdots \right)$$

$${\mathbf{H}}_{{\mathbf{x}}_b}{\mathbf{B}}^{\frac{1}{2}}{\mathbf{w}}_t\approx \mathbf{y}-\mathbf{h}({\mathbf{x}}_b)$$

$${\mathbf{H}}_{{\mathbf{x}}_b}(\mathbf{x}-{\mathbf{x}}_b)\approx \mathbf{h}(\mathbf{x}-{\mathbf{x}}_b).$$

$${\mathbf{x}}_a-{\mathbf{x}}_b\approx \mathbf{K}\left(\mathbf{y}-\mathbf{h}({\mathbf{x}}_b)\right),$$

Full text

\journalname

Submitted to Journal of Advances in Modeling Earth Systems (JAMES)

Department of Mathematics, University of Surrey, Guildford GU2 7XH, UK School of Archaeology, Geography and Environmental Sciences (SAGES), University of Reading, Whiteknights, Reading RG6 6AH, UK Department of Mathematics & Statistics, University of Reading, Whiteknights, Reading RG6 6AX, UK AXA Chair of Biosphere and Climate Impacts, Department of Life Sciences, Imperial College London, Silwood Park Campus, Buckhurst Road, Ascot SL5 7PY, UK

\correspondingauthor

Sean F. Cleators.cleator@surrey.ac.uk

{keypoints}

3D-Variational technique used to reconstruct seasonal palaeoclimates

Results reinforce spatial and seasonal coherence of temperature and moisture anomalies

Applied to southern Europe at the Last Glacial Maximum, ca 21000 years ago

1 Introduction

Past climates provide useful examples of how the climate system has responded to changes in external forcing, such as orbitally-induced changes in incoming solar radiation, and internal Earth system feedbacks, such as changes in atmospheric CO2 concentration ([CO2]) or ice sheet extent (Harrison \BBA Bartlein, \APACyear2012). Reconstructions of past climate states are now routinely used to evaluate the performance of the climate models that are used to project the trajectory of future climate changes (Harrison \BOthers., \APACyear2014; Schmidt \BOthers., \APACyear2014; Harrison \BOthers., \APACyear2015; Kageyama \BOthers., \APACyear2018). The Last Glacial Maximum (LGM, ca. 21,000 years ago) has been a major focus for these evaluations because the change in climate forcing was as large (albeit different in type) as ”high-end” changes projected for the end of the 21st century (Braconnot \BOthers., \APACyear2012; Kageyama \BOthers., \APACyear2018). These evaluations obviously depend on the availability of quantitative reconstructions of key climate variables and this has led to the creation of benchmark data sets documenting climate conditions over land (e.g. Bartlein \BOthers., \APACyear2011) and ocean (e.g. MARGO Project Members \BOthers., \APACyear2009).

Past climate conditions can be inferred from environmental records which respond to climate, including sedimentological, geomorphological, chemical, isotopic and biological records (Gornitz, \APACyear2008; Bradley, \APACyear1999). Quantitative reconstructions of climate variables can be obtained from these records either using statistical techniques based on modern day climate-response relationships (e.g. Ter Braak \BBA Juggins, \APACyear1993; see also discussion in Bartlein \BOthers., \APACyear2011) or by inversion of a model that simulates the response of a particular type of environmental record to climate (e.g. Garreta \BOthers., \APACyear2010; Steiger \BOthers., \APACyear2017). Pollen preserved in anoxic lake and bog sediments through time is the most widespread source of data for the reconstruction of terrestrial climates (Bartlein \BOthers., \APACyear2011; Marsicek \BOthers., \APACyear2018), because pollen abundance reflects the distribution of different plant taxa that have highly specific climatic requirements (Woodward, \APACyear1987; Harrison \BOthers., \APACyear2010) and the pollen-preserving sediments can be accurately dated using radiometric techniques. One important characteristic of all of the environmental records that are used for climate reconstruction, including pollen, is that both the primary data and the climate reconstructions are generated for individual sites. Geological and climatic factors mean that the distribution of potential sites is spatially non-uniform: speleothem records, for example, are confined to karst areas; pollen preservation requires anoxic conditions and thus pollen records are not common in arid regions. Furthermore, issues of accessibility and scientific interests means that the actual sampling of potential sites is non-uniform, so there are often large geographic gaps in the data coverage (Harrison \BOthers., \APACyear2016). While pollen records, for example, are well-sampled across Europe and North America, there are far fewer records from central Eurasia or the tropics. Furthermore, geological preservation issues mean that the number of sites available tends to decrease through time: there is an order of magnitude more pollen data available for the middle Holocene (ca 6000 yr BP) than for the LGM, for example Harrison \BOthers. (\APACyear2016). Ideally, a benchmark data set for model evaluation would provide continuous climate fields. However, while gridding the data sets at a scale comparable to that of the climate models (see e.g. Bartlein \BOthers., \APACyear2011) can improve the situation, this still does not solve the problem of significant gaps in site-based data coverage.

Alternative approaches to generating spatially continuous palaeoclimate reconstructions have been developed that involve combining observations with model simulations of palaeoclimates. Goosse \BOthers. (\APACyear2006), for example, used observations to select the most realistic member from an ensemble of climate-model simulations. They ran a relatively large ensemble of simulations using a range of different climate forcings and/or model parametrisations to encompass uncertainties, and then selected the members of the ensemble that best matched the observations at each time step before running these simulations for longer to gain an new estimate of the climate. In this approach, the most realistic climate is taken to be the simulated climate(s) that best matched observations after multiple simulations. Although this approach provides continuous and self-consistent fields of climate variables, the reconstructions cannot deviate fundamentally from the model predictions and thus could still be influenced by systematic errors inherent in the model construction. Annan \BBA Hargreaves (\APACyear2013) also used an ensemble of model simulations, but in this case they used multiple models. The ultimate climate reconstruction was assumed to be a weighted average of those climate models, where the weighting was determined by the goodness-of-fit to the observations. They applied a global weighting to each model rather than allowing the goodness-of-fit to vary regionally. As a result, there are regions where the reconstructed palaeoclimate is far from the observations, producing a palaeoclimate reanalysis that is highly influenced by systematic errors in the models.

Variational data assimilation techniques provide a way of combining observations and model outputs to produce climate reconstructions that are not explicitly constrained to a given source (Nichols, \APACyear2010; Lahoz \BBA Schneider, \APACyear2014). Variational techniques are widely used by the weather forecasting community (e.g. Daley, \APACyear1994) and have also been used to reconstruct palaeoclimate. Gebhardt \BOthers. (\APACyear2008), for example, used this approach to reconstruct European climates during the Last Interglacial. Simonis \BOthers. (\APACyear2012) applied the same basic approach to reconstruct January and July temperatures across European climate during the late Glacial (13,000 yr B.P.) and early Holocene (8,000 yr B.P.). The method involves applying a spatial constraint, based on a two-dimensional advection-diffusion equation of atmospheric dynamics, to upscale climate variables derived from statistical transfer functions relating the abundance of plant taxa with January and July temperature. In both examples, modern-day wind fields were used as the prior to determine the spatial scale and assumed to be the same in the past.

Tardif \BOthers. (\APACyear2018) also use variational techniques to create palaeoclimate reconstructions for the Last Millennium, using an ensemble of transient palaeoclimate simulations. They first determine the relationship between palaeoclimate reconstructions and the model-derived prior using linear regression, and then determine the best linear unbiased estimate (BLUE) using the Kalman formulation, to create the analytical reconstructions. Thus, temporal relationships are not based on an explicit analytical function designed to preserve structures (auto-correlations and/or discontinuities) in the prior and/or observations. Spatial correlation is generated from the prior ensemble with a covariance localisation applied to prevent spurious correlations. This, and the necessity to define scaling parameters, involves a number of arbitrary choices which influence the final reconstructions and make it difficult for these reconstructions to deviate substantially from the prior.

The 3D-Variational method finds the maximum a posteriori Bayesian estimate of the palaeoclimate given the site-based reconstructions and a prior estimate. While this could lead to the generation of reconstructions with sharp changes in time and/or space, it is possible to incorporate additional assumptions about the error of the prior estimate (the difference between the true climate and the prior) to prevent this by ensuring continuity of the solution. The degree of continuity in the change of the reconstructed climate field can be controlled by adjusting two length scales: a spatial length scale that determines how smooth the spatial correlation in the prior is between different geographical areas and a temporal length scale that determines how smooth it is through the seasonal cycle.

Here we apply this method to reconstruct six palaeoclimate variables across southern Europe at the LGM. The six climate variables are those provided in the Bartlein \BOthers. (\APACyear2011) dataset, namely mean annual temperature (MAT, ${}^{\circ}C$), mean temperature of the coldest month (MTCO, ${}^{\circ}C$), mean temperature of the warmest month (MTWA, ${}^{\circ}C$), growing degree days above a baseline of above $5^{\circ}$C (GDD5, $d^{\circ}C$), mean annual precipitation (MAP, mm) and an index of plant-available moisture (the ratio of actual to equilibrium evapotranspiration or $\alpha$ in Bartlein \BOthers. (\APACyear2011) re-expressed as a moisture index (MI, unitless) defined as the ratio of MAP to equilibrium evapotranspiration in our analyses. The conversion was made using the Zhang \BOthers. (\APACyear2004) formulation of the Budyko relationship). We use pollen-based reconstructions of climatic variables for the region of southern Europe (defined here as south of 50*∘N and extending eastward to 50∘*E) from Bartlein \BOthers. (\APACyear2011) as our observations. Although the sites from Europe were used to produce a gridded map in Bartlein \BOthers. (\APACyear2011), here we used the underlying individual site reconstructions. Some of the reconstructions used in Bartlein \BOthers. (\APACyear2011) were derived by model inversion, and these were excluded from our data set. Bartlein \BOthers. (\APACyear2011) gives mean values as anomalies from the modern climate, as well as standard errors. We use eight LGM climate simulations (CCSM4, CNRM-CM5, MPI-ESM-P, MRI-CGCM3, FGOALS-g2, COSMOS-ASO, IPSL-CM5A-LR, MIROC-ESM) from the 3rd phase of the Palaeoclimate Modelling Intercomparison Project (PMIP3: Braconnot \BOthers., \APACyear2012) to create a prior. These simulations were forced by changes in incoming solar radiation, changes in land-sea geography and the size and extent of ice sheets, and a reduction in atmospheric [CO2] (see Braconnot \BOthers., \APACyear2012 for details of the modelling protocol).

Our approach introduces features novel to the field of palaeoclimate data assimilation, explicitly designed to maximise the usefulness of the reconstructions for climate model evaluation. Specifically, by solving the full variational problem we take into account nonlinearities in the system. Furthermore we minimise the dependency of the final analytical reconstructions on the prior generated from the climate models by using a prescribed correlation function for the error of the prior and by using a resolution matrix (Menke, \APACyear2012; Delahaies \BOthers., \APACyear2017) to determine the temporal correlation length scale. The resolution matrix provides a particularly useful way to overcome problems caused by the sparsity of site-based palaeoclimate reconstructions at the LGM. In addition to investigating methods to determine appropriate spatial and temporal length scales, we provide a way of calculating the error in the final reconstructions.

2 Data Assimilation with Spatial and Temporal Correlations in the Prior

In this section we describe the underlying method used in this paper. Section 2.1 describes the inverse problem solved by the method and the types of data used. Section 2.2 shows how the different climate variables can be related to one another by specifying correlations from our prior estimate of the system. Finally section 2.3 describes how the problem is preconditioned in order to reduce the computation cost.

2.1 The Inverse Problem

Our problem is to determine the palaeoclimate that existed from a particular set of site-based reconstructions. We label the reconstructions as the column vector ${\bf y}_{i}\in\mathbb{R}^{6}$ for site $i$. For each reconstruction, ${\bf y}_{i}$, there are a total of $6$ variables that may have been reconstructed, namely; $\alpha$, MAP, MAT, MTCO, MTWA and GDD5. All these reconstructions together make the observations labelled ${\bf y}\in\mathbb{R}^{6N}$ such that

[TABLE]

where $N$ is the number of reconstructions. The reconstruction technique gives the variances for each reconstruction that we label as the column vector ${\bf v}_{y}\in\mathbb{R}^{6N}$ in the same order as ${\bf y}$. Not all variables are reconstructed at every site, for these variables the variance tends to infinity; this is achieved by setting their inverse to 0.

From these reconstructions we want to produce a gridded climate, the state vector, ${\bf x}\in\mathbb{R}^{13M}$ where there are $M$ grid cells. The $j$’th grid cell of the state is labelled ${\bf x}_{j}\in\mathbb{R}^{13}$ where

[TABLE]

For each grid cell the method determines a set of $13$ variables: the mean annual precipitation ($P$) and the $12$ average temperatures for each month, ${\bf T}$ where

[TABLE]

where $T_{m}$ is the temperature at month $m$.

For a general function ${\bf h}$ that maps a gridded climate ${\bf x}$ to the site-based observations we state the problem as trying to find an ${\bf x}$ such that

[TABLE]

Solving equation (3) for ${\bf x}$ is ill-posed as there are several ${\bf x}$ that are possible solutions. A prior estimate of the state called the background or prior (${\bf x}_{b}$) allows us to find the best ${\bf x}$ that solves equation (3) and remains close to the prior. The standard deviations of the prior are labelled as the vector ${\bf v}_{b}\in\mathbb{R}^{13M}$ in the same order as ${\bf x}_{b}$.

It can be shown (Nichols, \APACyear2010) that the optimal solution of equation (3) with a prior estimate of the state is defined as the analysis, ${\bf x}_{a}$, where

[TABLE]

with the cost function $J$ as

[TABLE]

Here ${\bf B}$ is the covariance of the uncertainties in the prior (conventionally denoted ${\bf B}$, for background) and ${\bf R}$ is the covariance of the uncertainties in the site-based reconstructions. Equations (4) and (5) ensure that the solution is the optimal distance from the observations subject to ensuring that the solution is not too far from the prior estimate, weighted by the error statistics in each. We assume that there are no correlations in the errors of the observations so we set

[TABLE]

The prior error covariance matrix can be represented as the product of the standard deviations of the prior and the correlations between the errors in the variables in the prior. Hence we write

[TABLE]

where

[TABLE]

is the diagonal matrix formed of the standard deviations of the prior error and ${\bf C}$ is the prior error correlation matrix.

2.2 Prior Error Correlation

The difference between the true ${\bf x}$ and the prior, the error in the prior, is expected to be smooth between adjacent grid cells and also from month to month since it would be unlikely that the observations would contain sharp jumps in climate that aren’t present in the prior. It would be unusual, for example, to have very high temperature in March if the temperatures in February and April are very low, if this behaviour isn’t seen in the prior. To achieve this we impose a structure on the prior error correlation matrix, ${\bf C}$, that weighs the cost function so that its minimum is smooth. This allows the prior error to be smooth, but still allows non-smooth areas if there is significant evidence to support it in the prior and/or the observations.

We assume there are two independent sets of correlations in the prior. The first correlation is spatially between the different grid cells. We also assume that the spatial correlation between the grid cells is homogeneous and valid on a sphere, so that for an angle $\theta_{ij}$ on a great circle of the Earth between each cell $i$ and $j$ the correlation is given by,

[TABLE]

where $c_{L}$ is a case of a Matérn function (Matérn, \APACyear1986; Handcock \BBA Wallis, \APACyear1994) with order 1 and $\mathcal{K}$ is the modified Bessel function of the second kind, evaluated using the boost C++ library (Maddock \BOthers., \APACyear2018). Here the correlation length scale is $L=L_{s}$ and $a=6371$km is the radius of the Earth. The correlation matrix between all grid cells, ${\bf C}_{L_{s}}$, is given as

[TABLE]

The choice of $L_{s}$ is dependent on the datasets used in ${\bf y}$ and ${\bf x}_{b}$ and so is specific to each problem. In section 3.2 a method of finding $L_{s}$ is shown for a particular experiment.

The second assumed correlation is between the error in the average temperatures of the prior. We assume that there is a correlation between the average temperatures of a month and the surrounding months given by equation (8). Here $\theta_{ij}=mod_{12}(\left|i-j\right|)$ between months $i$ and $j$. The correlation length scale is $L=L_{t}$ and $a=6/\pi$. The appropriate value of $L_{t}$ again depends on the datasets given and is shown for a particular experiment in section 3.2. For each grid cell the correlation between the different climate variables is given by ${\bf C}_{L_{t}}$ where

[TABLE]

Note how $\{c_{L_{t}}(\theta_{ij})\}_{ij}$ is offset by the first row and column due to the presence of the precipitation term which is uncorrelated to the temperature terms.

These two sets of correlations imply that all the variables in the error of the prior are correlated. For instance the grid cells $i$ and $j$ are correlated by $\left({\bf C}_{L_{s}}\right)_{ij}$ and the temperatures in month $l$ and $k$ are correlated by $\left({\bf C}_{L_{t}}\right)_{lk}$. This means that the temperatures in month $l$ in grid cell $i$ and month $k$ in grid cell $j$ are correlated by the product $\left({\bf C}_{L_{s}}\right)_{ij}\left({\bf C}_{L_{t}}\right)_{lk}$. Repeating this for every variable gives an overall correlation for the prior (${\bf C}$ from equation (6)) as

[TABLE]

where $\otimes$ is the Kronecker product of matrices.

The incorporation of correlations structures is due to the fact that the state space covers space and time. We introduce the ${\bf C}_{L_{s}}$ and ${\bf C}_{L_{t}}$ correlations to make the prior error smooth in space and time respectively. The presence of the scales $L_{s}$ and $L_{t}$ allows the adjustment of the smoothing in both dimensions and should depend, at least in part, on the spatial and temporal distribution of the prior and site-based reconstructions. In section 3.2 we discuss methods for choosing these scales.

2.3 Preconditioning and the Condition Number

The minimum of the cost function is sensitive to change in the input data of the problem and to computational errors. This sensitivity reflects the difficulty in solving the problem and is measured by the condition number of the Hessian of the cost function (Golub \BBA Loan, \APACyear1996). We define the condition number $\kappa$ of a symmetric positive definite matrix ${\bf M}$ to be

[TABLE]

where $\lambda_{\text{max}({\bf M})}$ and $\lambda_{\text{min}({\bf M})}$ are the maximum and minimum eigenvalues of ${\bf M}$. Here, ${\bf M}$ is the Hessian of the cost function, given by its (first order) second derivative ${\bf S}={\bf H}{\bf B}{\bf H}^{T}+{\bf R}$. This condition number indicates the computational work needed to minimise the cost function. Equation (11) shows how the condition number of ${\bf S}$ represents the disparity in scales of the problem. As the eigenvalues represent the sizes of the scales of ${\bf S}$, their ratio represents the largest scale that will be encountered when inverting ${\bf S}$. Since large scale differences create more numerical inaccuracy, a large condition number will increase the computational cost and lead to an inaccurate solution.

Haben \BOthers. (\APACyear2011) shows that the bounds on the condition number can be reduced by minimising the cost function around ${\bf w}$ instead of ${\bf x}$ where

[TABLE]

where ${\bf B}^{\frac{1}{2}}$ is the symmetric square root of the matrix ${\bf B}$ such that

[TABLE]

The use of this linear transformation can be thought of as a z-score to work with uncorrelated states.

Equation (12) transforms the inverse problem from equation (4) into finding

[TABLE]

where $J({\bf w})$ is

[TABLE]

We use the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method to find the state, ${\bf w_{a}}$, which has the minimum $J$, L-BFGS is a quasi-Newton method that maintains a limited memory version of an approximated Hessian as described in Liu \BBA Nocedal (\APACyear1989). At each evaluation step we calculate the gradient of $J$ as

[TABLE]

where ${\bf H}_{{\bf x}}$ is the Jacobian of ${\bf h}$ evaluated at ${\bf x}$. Once ${\bf w}_{a}$ is found we use equation (12) to find ${\bf x}_{a}$, the solution.

The error covariance of the analysis, ${\bf x}_{a}$, is denoted by ${\bf A}$ and is given (to first order) following Nichols (\APACyear2010) as

[TABLE]

where the gain matrix ${\bf K}$ is

[TABLE]

3 Experimental Design

We use our method to reconstruct the palaeoclimate of southern Europe during the Last Glacial Maximum (LGM). The LGM had insolation forcing relatively similar to the present day but northern hemisphere ice sheets were more extensive, sea-level was lower and the area of the continents therefore larger, and the atmospheric [CO2] was less than half of the concentration today. In this section we describe the choices of ${\bf h}$, ${\bf y}$ and ${\bf x}_{b}$ used to make this reconstruction and our choices for $L_{t}$ and $L_{s}$, the correlation length scales.

3.1 Experiment Setup

We use pollen-based reconstructions of climatic variables from Bartlein \BOthers. (\APACyear2011) as our observations. Bartlein \BOthers. (\APACyear2011) gives means as anomalies from the modern climate as well as standard errors. We add the anomalies to the CRU CL v2.0 dataset (New \BOthers., \APACyear2002) to derive absolute climate reconstructions. We non-dimensionalise the climate variables in order to avoid computational issues because they are on different scales in the calculation of the cost function. After solving for the non-dimensionalised case, we re-dimensionalise each of the variables to be on the original scale of the observations and the prior. Details of the dimensionalisation and non-dimensionalisation of variables can be found in Appendix A. We use the non-dimensionalised variables as our ${\bf y}$ and their non-dimensional standard errors, formed from the product of the standard errors and the derivative of $D_{y}$ (equation 22), as ${\bf v}^{\frac{1}{2}}_{{\bf y}}$.

We use the LGM outputs from PMIP3 as our prior. We use the variables of monthly precipitation (that are summed to annual precipitation), monthly temperature and monthly total cloud fraction. For each of the selected PMIP models that ran an LGM experiment we interpolate the output to a $2^{\circ}\times 2^{\circ}$ grid producing a set of maps all at the same resolution. In order to minimise the impact of potential individual systematic model biases the simulated climate at the LGM, experiments are generally expressed relative to that specific model’s pre-industrial control (PI) experiment. We therefore interpolate each of the PI experiments to the same grid and take the difference between the LGM and PI experiments of each model as the anomaly to the modern day for each model. We then sum each model’s anomalous values with values from the modern day (from CRU CL v2.0, as above, bilinearly interpolated to the $2^{\circ}\times 2^{\circ}$ grid) in order to produce absolute values for each model. For each variable in the set we take the mean and variance across the set of all models to produce a gridded map. As for the observation space, we non-dimensionalise the state space to remove any dimensional effects using equation (24). The non-dimensional variables form the prior ${\bf x}_{b}$ and their non-dimensional variances, formed from taking the product of the variances and the derivative of equation (24), form ${\bf v}_{{\bf x}_{b}}$.

The observation function, ${\bf h}$, links together the variables from both datasets. At each site, $i$, we define the observation function as

[TABLE]

The derivatives, $\frac{\partial max(T)}{\partial T_{m}}$ and $\frac{\partial max(T)}{\partial T_{m}}$ are taken to be $1$ if $T_{m}$ is the maximum or minimum of $T$ and [math] elsewhere. The moisture index function $\mu$ is

[TABLE]

as given by the Budyko curve with $\omega=3$ as described in Zhang \BOthers. (\APACyear2004). The moisture index $m$ is calculated as

[TABLE]

where $\gamma$ ($0.067kPaK^{-1}$) is the psychrometer constant at sea level, $l_{j}$ is the length of month $j$ in days and where

[TABLE]

is the differentiated Roche-Magnus formula from Allen \BOthers. (\APACyear1998). The function $R(T_{k},S_{k})$ is the daily net radiation at the vegetated surface defined in Davis \BOthers. (\APACyear2017) for the middle day in month $k$. The variable $S_{k}$ is the total cloud fraction for month $j$ which is taken from the PMIP3 average described above. We define

[TABLE]

and the mean function to be $mean(\bar{{\bf T}})=\frac{1}{N_{y}}\sum_{k}^{12}l_{k}\bar{T}_{k}$ and $max(\bar{{\bf T}})$ and $min(\bar{{\bf T}})$ to be the maximum and minimum temperature in $\bar{{\bf T}}$ respectively. The full observation function, ${\bf h}$, is formed by applying $\hat{{\bf h}}$ at each grid cell where there is an observation and defining

[TABLE]

and so ${\bf h}$ will have the dimension $6N$, and hence the Jacobian of $h$, $H$, will have dimension $(6N)^{2}$.

3.2 Determining $L_{t}$ and $L_{s}$

The two correlation length scales, $L_{t}$ and $L_{s}$, in ${\bf C}$ (section 2.2) determine the strength of the correlation in the errors of the prior. By varying the length scales we can vary how smooth the error of the prior is and hence how smooth the solution is. If the length scale is too large then the error will be over-smoothed and the solution will miss smaller scale features such as inter-annual temperature changes or spatially small features such as topography. A length scale too small will mean the solution will be too coarse and contain unrealistic jumps.

In order to determine a suitable value for $L_{t}$ we consider a single grid cell with a single simulated observation at $37.50^{\circ}$N and E$33.73^{\circ}$, which allows us to ignore the effects of ${\bf C}_{L_{s}}$. The example only has observations of MTCO and MTWA ($-15^{\circ}$C and $30^{\circ}$C respectively), allowing us to ignore the non-linear effects of calculating $\alpha$. Fig. 1 shows the prior and observations for the sample as well as the estimated states after assimilation for different values of $L_{t}$. For all values of $L_{t}$ the analysis doesn’t match the observed MTCO since the prior temperature for January has low uncertainty. Low values of $L_{t}$ create an analysis that swaps between the prior and the observations. Although the solution always matches either the reconstructions or the prior, the jumps between them are unrealistic. On the other hand high values of $L_{t}$ create an analysis that follows the prior too closely and is unable to create high and low temperatures. The value of $L_{t}=1$ produces an assimilation that follows the shape of the prior but lies between the values of the prior and the observations.

We can further understand $L_{t}$ by seeing how information is changed by the method. If we consider the hypothetical, true solution to the inverse problem, ${\bf w}_{t}$, then by equation (3) we have that

[TABLE]

since, up to first order,

[TABLE]

Further Nichols (\APACyear2010) shows how

[TABLE]

where ${\bf K}$ is the gain matrix defined in equation (17). Hence we can consider the change from true solution to our computed one (${\bf w}_{a}$) as being given by

[TABLE]

where

[TABLE]

is the resolution matrix as described in Menke (\APACyear2012); Delahaies \BOthers. (\APACyear2017).

Resolution matrices where the diagonal elements are close to 0 describe a situation where, if perfect information is input, then the solution would only contain part of this information. In situations where the resolution matrix has large off-diagonal terms, the solution is degraded by interference between variables. If the opposite is true, the resolution matrix is close to the identity matrix. The best method will have a resolution matrix that resolves as many variables as possible whilst having few variables interfering with each other.

Fig. 2 shows how the resolution matrix changes with respect to $L_{t}$ for the same test grid cell as in Fig. 1. The simulated prior temperatures are closest to the observations in January and July such that for small values of $L_{t}$, the method resolves temperatures in these months well. However, for large $L_{t}$ the method improves the patterns away from these months whilst degrading reconstructions of January and July. Values of $L_{t}$ in between the large and small values show a mixture of both high resolution and low interference. These results together with the results from Fig. 1 suggest a value of $L_{t}=1$ is suitable for this problem.

The choice of the other scale, $L_{s}$, is especially relevant for the relatively sparse dataset used here. A higher $L_{s}$ represents errors in the prior being correlated even though they are far away, whereas a low $L_{s}$ represents errors not being highly correlated even though they are close together. A large $L_{s}$ means that information from the reconstructions could be be propagated over a large distance. While this is useful in maximizing the use of a geographically sparse data set, it could be unrealistic if this propagation extends too far beyond the source area for the pollen on which the site reconstructions are based (which is generally, though not always, of the order of $20-100$km around the site). In order to obtain a realistic solution whilst maximising the use of the data we choose $L_{s}$ such that the assumed average source area of the different sites does not overlap.

$L_{s}$ corresponds to the area that each observation impacts, so an increase in $L_{s}$ gives higher utilisation of observations. Haben \BOthers. (\APACyear2011) show that the condition number of the inverse problem is proportional to the distance between the reconstruction sites which, in this case, is proportional to $L_{s}$. However, the condition number corresponds to the sensitivity of inverting the Hessian to inputs and so is inversely proportional to the computational accuracy of the problem, up to first order. Hence, it is important to check that a choice of large $L_{s}$ doesn’t lead to a condition number for the problem that is too large to give an accurate result.

Fig. 3 plots $\kappa({\bf S})$ against $L_{s}$ and shows how $\kappa({\bf S})$ begins to increase with higher $L_{s}$. Also Fig. 3 shows several inflection points which could indicate values of $L_{s}$ that allow multiples of observations to interact. For this paper we pick a value of 400km for $L_{s}$ as this is large enough to propagate information sufficiently far from the different reconstructions. As seen in Fig. 3, $L_{s}=400$km still has a relatively low condition number and hence the solution will be relatively accurate.

4 Results

The solution using scaling values of $L_{t}=1$ and $L_{a}=400$ (Fig. 4) produces climates at 50 sites and surrounding grid cells that are close to the reconstructions, as expected, over much of the region. However, this is not the case for the MI values of the 3 sites at the eastern tip of the Black sea (Apiancha, Kobuleti, Sukhumi). These discrepant cases occur either where there is significant disagreement between different reconstructions and/or disagreement between the reconstructions and the prior with at least one of the reconstructions having relatively low variance. This reconstruction is weighted highly in the cost function and the solution does not meet the other reconstructed variables or the prior. This creates a situation in which the best possible solution differs from both the reconstructions and prior.

The difference between the solution and the prior, transformed by equation (18) at each grid cell and dimensionalised via equation (22), shows that the climate is much drier than the prior in the western part of the area, as shown by MI and precipitation (Fig. 5). MAT has increased in some regions but decreased in others; this suggests that the inclusion of ${\bf C}_{L_{s}}$ is working as intended, since although there are varied changes in MAT, the changes occur in a spatially smooth way. Furthermore there has been an increase in temperature seasonality as MTCO has become colder at all sites and MTWA has become warmer at most sites. This, together with the changes to MAT and GDD5 suggests that ${\bf C}_{L_{t}}$ is having the desired effect; as the changes to MTCO and MTWA are impacting the whole of the seasonal cycle of the climate and giving reasonable and smooth changes to both MAT and GDD5.

In general (Fig. 6) grid cells near reconstruction sites have less error, because the solution is using information from both the prior and the reconstructions, while grid cells further away from reconstruction sites have higher error by defaulting to the error in the prior. There are some areas near reconstruction sites with high errors in MTCO, particularly in the northeast. This could reflect the fact that vegetation towards the cold and dry end of the winter temperature gradient is less sensitive to temperature change than in the Mediterranean region. However, the high median error for MTCO overall shows that there need to be large changes in MTCO from the prior to match the observations.

5 Discussion

Our final temperature reconstructions show good coherence spatially, plausible seasonal relationships, and no systematic discrepancies from pollen-based reconstructions at individual sites. However, the reconstructions of moisture variables, MAP and MI, are wetter than indicated by the pollen-based reconstructions. This was expected and is realistic. The atmospheric CO2 concentration, [CO2], was considerably lower during the LGM than it is today (180 ppm compared to 280 ppm in the PI simulations, and ca. 400 ppm today). Low [CO2] decreases the water-use efficiency of plants and favours drought-adapted plants at the expense of trees, even without a change in climate (Jolly \BBA Haxeltine, \APACyear1997; Prentice \BBA Harrison, \APACyear2009). Although there are methods of accounting for this direct [CO2] effect (Prentice \BOthers., \APACyear2017), statistical techniques based strictly on the application of modern analogues do not account for this impact. All of the pollen-based reconstructions for southern Europe from the Bartlein \BOthers. (\APACyear2011) data set are based on statistical reconstruction techniques. Application of the theoretically-based correction factor derived by Prentice \BOthers. (\APACyear2017) to the reconstructed moisture variables would be a useful next step to improve their realism.

Sites suitable for obtaining pollen records are not uniformly distributed geographically, and in any case the actual sampling of potential environments is extremely uneven in many regions of the world (Figure 4; Bartlein \BOthers., \APACyear2011; Harrison \BOthers., \APACyear2016). We have shown that the condition number can be used to identify an appropriate scale for interpolating the site-based data spatially, and that a scale of 400-500km appears to be appropriate for southern Europe at the LGM given the data currently available. This spatial scale is not uniformly appropriate, however. The standard deviation of the reconstructions (Fig. 6) provides a measure of how reliable the interpolation is. More importantly, the standard deviation of the reconstruction could be used to determine when the interpolated values provide a realistic measure of the actual climate and when they do not. Establishing an acceptable threshold value for reliability would be a useful step in the creation of the kind of palaeoclimate reanalysis we are proposing here.

Whilst the values of both scales, $L_{s}$ and $L_{t}$, have been shown to be appropriate for the example shown in this paper, they are somewhat subjective. The spatial scale, $L_{s}$, is chosen to give high utilisation of sparse observation data and is shown, by the condition number in Fig. 3, not to lead to a numerically inaccurate solution. A value for $L_{t}$ is determined by plotting the resolution matrix for multiple $L_{t}$, as shown in Fig. 2; however, this only provides a range of possible values. A more objective method for selecting $L_{t}$ could be developed by selecting the $L_{t}$ which gives the resolution matrix closest to the identity.

6 Conclusions

In this paper we have demonstrated a novel method for reconstructing spatially explicit palaeoclimate reconstructions from site-based data. The method allows the effects of each site in the dataset to be tuned by imposing a structure on the error of the prior that creates reconstructions that are spatially smooth and hence more realistic. By assuming that the error in the prior with respect to temperature has a given correlation month by month, it also allows the generation of a solution that is temporally smooth. We show that a length scale $L_{t}$ of 1 provides a smooth solution for the seasonal cycle, both using single sites and over multiple grid cells. Our analyses suggest that a spatial length scale ($L_{s}$) of 400km is reasonable for southern Europe at the LGM; although this is larger than the assumed source area of most of the reconstruction sites, it reflects the large-scale coherence of the regional climate change between LGM and present. Additional work could help to determine a more objective way to determine these length scales, but nevertheless the final climate maps appear plausible and suggest that the application of this new method should yield more robust data sets for climate-model evaluation.

Appendix A Non-dimensionalisation

Most of the variables from the site-based reconstructions and PMIP3 have a dimension. This can cause a problem when computing the cost function as different variables can be at different scales and it is difficult to compare different scales together computationally. To avoid this problem we non-dimensionalise all the variables involved before computing the cost function and then re-dimensionalise the variables when the analysis has been found.

We non-dimensionalise the observation space using

[TABLE]

where $N_{y}$ is the number of days in a year, $T_{s}$ is a temperature scaling value ($5^{\circ}C$). The function $D_{P}$ is defined as

[TABLE]

where $I_{sc}$ is the solar constant ($1360.8Wm^{-2}$) and $\lambda$ is the latent heat of vaporisation of water ($2.45MJkg^{-1}$). $D_{P}$ ensures that the method never creates a situation where $P<0$. Similar to the observation space, we also non-dimensionalise the state space using

[TABLE]

Acknowledgements.

SFC was supported by a UK Natural Environment Research Programme (NERC) scholarship as part of the SCENARIO Doctoral Training Partnership at the University of Reading. SPH acknowledges support from the ERC-funded project GC 2.0 (Global Change 2.0: Unlocking the past for a clearer future, grant number 694481). ICP acknowledges support from the ERC under the European Union’s Horizon 2020 research and innovation programme (grant agreement No: 787203 REALM). This research is a contribution to the AXA Chair Programme in Biosphere and Climate Impacts and the Imperial College initiative on Grand Challenges in Ecosystems and the Environment (ICP). NKN is supported in part by the NERC National Center for Earth Observation (NCEO). We thank PMIP colleagues who contributed to the production of the palaeoclimate reconstructions. We also acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and the climate modeling groups in the Paleoclimate Modelling Intercomparison Project (PMIP) for producing and making available their model output. For CMIP, the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. The analyses and figures are based on data archived at CMIP on 12/09/18. We thank the Next-Generation Vegetation Modelling group for providing model code for the calculation of bioclimatic variables and for discussion of the results.

Code and data availability

The reconstructions of southern Europe in this paper, as well as the code used to perform the data assimilation, are currently being archived at the University of Reading Research Data Archive (https://researchdata.reading.ac.uk/). The CRU CL v2.0 dataset was downloaded from the University of East Anglia Climatic Research Unit and was published as New \BOthers. (\APACyear2002). The PMIP3 LGM simulations (Braconnot \BOthers., \APACyear2012) are available at CMIP5 archives and, for this paper, were downloaded from the Earth System Grid Federation at the Pierre Simon Laplace Institute (https://esgf-node.ipsl.upmc.fr/projects/esgf-ipsl/). The pollen reconstructions used are available from Bartlein \BOthers. (\APACyear2011).