2.1 Model simulation: CCC400

We start from an existing DA system, which is described in Bhend et al. (2012) and Franke et al. (2017a). We use the same atmospheric model simulation as in the previous studies. The model simulation, termed as Chemical Climate Change over the Past 400 years (CCC400), has 30 ensemble members that are used as background to reconstruct monthly climate states between 1600 and 2005. Simulations were performed with the ECHAM5.4 climate model (Roeckner et al., 2003) at a resolution of T63 with 31 levels in the vertical. The 30 ensemble members were forced and driven with the same external forcings and with the same boundary conditions. For sea-surface temperatures (SSTs), which have a particularly large effect on the simulations, the reconstruction by Mann et al. (2009) was used. At the time when the model simulation was run, this was the only available global gridded SST reconstruction that dated back until 1600. The surface temperature reconstruction by Mann et al. (2009) is based on a multiproxy network and was produced by a climate field reconstruction method. The SST reconstruction by design captures interdecadal variations (Mann et al., 2009); hence, intra-annual variability dependent on a El Niño–Southern Oscillation reconstructions (Cook et al., 2008) was added to the SST fields. Further forcings include solar irradiance, volcanic activity and greenhouse gas concentrations (for more details, see Bhend et al., 2012; Franke et al., 2017a). The land-use reconstruction by Pongratz et al. (2008) was used to derive the land-surface parameters. The 6-hourly output fields provided by the model were transformed to monthly means. To reduce the computational burden, only every second grid point in latitude and longitude was selected. We limit the analysis in this study to 2 m temperature, precipitation and sea-level pressure.

2.2 Observational network

In this study, we use the same observational network of tree-ring proxies, documentary data and early instrumental measurements as described in Franke et al. (2017a) (Fig. 1), but we only assimilate tree-ring proxies and instrumental data. The temporal resolution of the instrumental air temperature and sea-level pressure measurements is monthly. The tree-ring proxy records have annual resolution. Trees respond to a locally varying growing seasons. We consider temperature from May to August and precipitation from April to June to possibly affect tree-ring width data. The maximum latewood density proxies were considered to be affected by temperature over May until August. The observations were quality checked before the assimilation, and outliers which were more than 5 standard deviations away from the calculated 71-year running mean were discarded, both for instrumental and proxy data.

Figure 1The observational network in 1904, before the quality check.

2.3 Assimilation method

In our paleoclimate reconstruction, we combine the CCC400 model simulation with the observations as described above by implementing a modified version of the ensemble square root filter (EnSRF; Whitaker and Hamill, 2002). This ensemble-based DA method is called ensemble Kalman fitting (EKF; Franke et al., 2017a). In fact, the EKF is an offline version of the EnSRF, and the update step of the EKF remains the same as of the EnSRF. EKF is described in more detail in Bhend et al. (2012) and Franke et al. (2017a). Here, we shortly highlight the most important aspects of the EKF. The update step in the EnSRF scheme has two parts: updating the mean ($\stackrel{\mathrm{‾}}{\mathit{x}}$), and for each member, the deviation from the mean (x^′). They are calculated as

where K and $\stackrel{\mathrm{̃}}{\mathbf{K}}$ are

The background state vector (x^b) contains the variables of interest from CCC400 (Table 1). In the EKF, the length of the assimilation window is 6 months (October–March and April–September), which was adapted to the Southern and Northern Hemisphere growing seasons to effectively incorporate the proxy records stored in trees. Due to the 6-monthly assimilation window, x^b contains the variables of 6 months. x^a stands for the analysis state vector. H is the forward operator that maps the model state to the observation space (here, it is linear). H differs depending on the type of observation being assimilated. In the case of tree-ring width data, H extracts temperature between May and August and precipitation between April and June from the model; these fields are transformed to observational space by using a multiple regression approach (for more details, see Franke et al., 2017a). y represents the observations. K is the Kalman gain matrix, and $\stackrel{\mathrm{̃}}{\mathbf{K}}$ is the reduced Kalman gain matrix. P^b is the background-error covariance matrix, estimated from the 30 ensemble members. A common assumption is to treat the observation-error covariance matrix (R) as a diagonal matrix: it is presumed that the elements of R are uncorrelated. Therefore, the observations can be processed serially. We set the error variances of instrumental temperature observations to 0.9 K² and of instrumental pressure data to 10 hPa². The error variances are rough estimates that include, for instance, measurement uncertainties, temporal inhomogeneities and the fact that a station is not representative for a grid cell (see Frei, 2014; Franke et al., 2017a). The errors of tree-ring proxy data are calculated as the variance of the multiple regression residuals of the H operator. The assimilation is conducted on the anomaly level: we subtract both from model and from observational data their 71-year running mean in order to deal with the biases related to systematic model errors and inconsistent low-frequency variability in the paleoclimate data.

Table 1Defined localization length scale parameters.

Download Print Version | Download XLSX

The use of DA in an offline manner is typical in paleoclimate reconstructions (e.g., Dee et al., 2016). Bhend et al. (2012) argue that the assimilation step is too long for initial conditions to matter, whereas there is some predictability from the boundary conditions. In addition, Matsikaris et al. (2015) found that both online and offline DA methods perform similarly in their paleoclimate reconstruction setup. Furthermore, the offline DA is advantageous as it allows using the precomputed simulations. In our case, we can use CCC400 (Bhend et al., 2012) and test the method without having to repeat the simulations.

2.4 Spatial localization

As R is a diagonal matrix, the EKF can be used to assimilate the observations one by one; that is, after the first observation is assimilated and the analysis is obtained, this analysis field becomes the background state for the next observation (see the arrow pointing from x^a to x^b on Fig. 2). This serial implementation makes the calculation of P^b simpler. H becomes a vector (not a matrix) of the same length as x^b. It is zero everywhere except for a few elements (those required to model the observation). This translates to only a few columns of P^b that are actually required. HP^bH^T and R are then scalars (Whitaker and Hamill, 2002). This procedure also makes the localization simpler, as it needs to be applied only to those columns. In the original setup, the elements of P^b were multiplied (Schur product) with a distance-dependent function (see Eq. 7 in Franke et al., 2017a). For all the variables in the state vector, the same Gaussian function was used but with different localization length scale parameters (Table 1). The localization length scale parameters are defined based on the spatial correlation of the variables in the monthly CCC400 model simulation fields. For the cross-covariances between two variables, the smaller localization length scale of the two variables is applied. With the serial implementation, the calculation and localization of P^b is significantly simplified.

Figure 2The main steps of the blending experiment in one assimilation window. The blended covariance matrix P^blend is calculated as a linear combination from the year specific and climatological covariance matrices. The calculation of the Kalman gain (K) and reduced Kalman gain ($\stackrel{\mathrm{̃}}{\mathbf{K}}$) matrices is the same as in Eqs. (3) and (4) except the covariance matrix is replaced with P^blend. The observation is assimilated to both state vectors and these analysis become to the starting point for assimilating the next observation.

Download

3.1 Spatial localization

In most of the studies, the localization function is implemented in an isotropic manner. In the original setup, the same horizontally isotropic localization function was used with different localization parameters. However, such spatial symmetries may not be realistic. In the real atmosphere, correlation lengths might be longer in the zonal than in the meridional direction, due to the prevailing winds and the weaker large-scale temperature gradients in this direction. On multi-annual to multi-decadal timescales, multiple processes act in the meridional direction, e.g., a widening/shrinking of the Hadley cell, shifts of the Intertropical Convergence Zone or changes in atmospheric modes like the Atlantic Multi-Decadal Oscillation or the North Atlantic Oscillation. These can shift the zonal circulation northward or southward, but the zonal coherence will be less effected. Hence, instead of using a circular Gaussian function, we conducted an experiment with a spatially anisotropic localization function

where d_z and d_m are the distances from the selected grid box in the zonal and meridional directions, respectively. L_z and L_m are the length scale parameters used for localization in the zonal and meridional directions, respectively. As a first experiment, we tested a 2:1 ratio for L_z:L_m. We used the values from Table 1 in the meridional direction and doubled them in the zonal direction. Thus, the resulting localization function has an elliptical shape.

3.2 Inflation techniques

Covariance inflation techniques are another possible method to compensate for errors in the DA system (Whitaker et al., 2008). The multiplicative inflation technique uses a small factor γ(γ>1) with which the x^′b is multiplied (Anderson and Anderson, 1999). This type of covariance inflation accounts for filter divergence due to sampling error (Whitaker and Hamill, 2002) but can be also applied to take into account system errors (Whitaker et al., 2008). We conducted some experiments using multiplicative inflation, although in our offline approach, filter divergence is not the main concern as ensemble members are not propagated in time.

The other methodology that we adapt shows similarities with the additive inflation technique (e.g., Houtekamer and Mitchell, 2005) and with the hybrid DA scheme (e.g., Clayton et al., 2013). In both methods, P^b is modified by either adding model error (Houtekamer and Mitchell, 2005) or a so-called climatological covariance matrix (Clayton et al., 2013) to P^b. This has given rise to the idea of generating a climatological ensemble in order to alleviate the effect of the small ensemble size. In the original setup, P^b is approximated from only 30 members. Here, we additionally build a climatological state vector (x^clim) from randomly selected ensemble members from our 400-year long model simulation. The number of ensemble members should be higher than the original ensemble size, but still computationally affordable. The climatological state vector is created as follows: (1) define the ensemble size (n) of x^clim; (2) select n random years between 1601 and 2005; (3) every year has 30 members from which one member is randomly selected and kept; (4) the chosen members are combined in x^clim. x^clim is randomly resampled after every second assimilation cycle. Using x^clim in the assimilation leads to increased computational cost, which partly comes from the creation of x^clim. The other time-consuming part comes from the updating of the climatological part after each observation is assimilated (the standard way when observations are assimilated serially). We tested numbers between 100 and 500. From x^clim, a climatological background-error covariance matrix (P^clim) can be obtained by using the ensemble perturbations. The background-error covariance matrix used in the blending experiments (P^blend) is built as a linear combination of the sample covariance matrix (P^b) and the climatological covariance matrix (P^clim):

where β₁, β₂ are the weights given to the covariance matrices. The sum of the weights is unity.

Figure 2 shows the main steps of the blending assimilation process. First, the covariance matrices were localized separately; then, we blended them according to the given weights. We conducted several experiments to tune the ratio between the two covariance matrices while using different localization length scale parameters (L) (Table 2). We used the same L values for localizing P^b in most of our experiments to evaluate improvements in comparison with the original setup. For this study, we calculated the latitudinal dependency of correlation of the state variables from a bigger ensemble of the model than in Franke et al. (2017a). The result suggested that longer L values can be applied in the tropics and the L of precipitation is probably too strict. Based on the rather strict L values in the previous study and the assumption that the covariances can be better estimated from a bigger ensemble, we used doubled length scale parameters (2L) in some of the experiments for localizing the climatological covariances. In this case, the L for temperature is 3000 km, which means that the correlation is decreased close to zero approximately 6000 km away from the observation.

Since observations are assimilated serially, we also update x^clim after an observation is assimilated with the same Kalman gain matrices as x^b. Thus, in the assimilation process, we update 30+n ensemble members.

3.3 Temporal localization

Localizing observations in time is a special feature of the EKF due to its 6-month assimilation window. Having the state vector in half-year format, every month within the October–March or April–September time window is updated by each single observation. To test whether the covariances between a single observation and the multivariate climate fields are correctly captured, we ran an instrumental-only experiment with temporal localization. We set covariances between different months to zero.

3.4 Skill scores

The EKF method is tested with different localization functions and with a set of mixed background-error covariance matrices as described above. We have performed the experiments by assimilating either only proxy records (proxy-only experiment) or only instrumental data (instrumental-only experiment). The proxy-only experiments were carried out between 1902 and 1959, because many proxy records already end in the 1960s, while the instrumental-only experiments were tested over the 1902–2002 period. We separated the different observation types to see whether different settings perform better depending on the type of data being assimilated. We do not compare proxy-only results with instrumental-only results; hence, the difference in time periods used does not matter; we simply use the longest possible time period. To evaluate the reconstructions, we examined two verification measures: correlation coefficient and reduction of error (RE) skill score (Cook et al., 1994). We use the Climate Research Unit (CRU) TS 3.10 dataset (Harris et al., 2014) for reference in the validation process. The presented verification measures are functions of time. Correlation is calculated between the absolute values of the ensemble mean of the analysis and the reference series at each grid point. The RE compares, in our case, the reconstruction with the model simulation, both expressed as deviations from a reference.

where x^u is the ensemble mean of the analysis, x^f is the ensemble mean of the model background state, x^ref is the reference dataset, and i refers to the time step. The RE skill scores are computed based on anomalies with respect to the 71-year running climatologies. Note that x^f comes from a forced model simulation; therefore, it already has skill compared with a climatological state vector. The RE is 1 if the x^u is equal to x^ref. Negative RE values indicate that the background state is closer to the reference series than the analysis.

To test which experiments have significantly different skill compared with the original skill, we carried out a permutation test following the method described in DelSole and Tippett (2014). Permutation was performed 10 000 times. If the difference between the median of the experiment and the median of the original data falls outside of the 95 % confidence level of the interval calculated from permuted data, then the experiment is significantly different from the original data.

In the next section, we will focus on analyzing the result of the experiments mainly over the extratropical Northern Hemisphere (ENH), because most of the data are located in this region. The skill scores refer to seasonal averages of the ensemble mean.

4.1 Localization function

4.1.1 Results

We compared the original setup applying an isotropic localization function and the experiment in which an anisotropic localization function was used, to test whether we can obtain a more skillful reconstruction by implementing anisotropic localization method. As an example of the spatial reconstruction skill, we show the RE values of temperature (Fig. 3). The figures reveal that the type of localization function only resulted in small differences in both experiments. Nonetheless, there are larger areas of negative RE values (Greenland, Siberia) with the anisotropic localization function in the proxy-only experiment (Fig. 3). In the instrumental-only experiment, the decrease of RE values occur in the northern high latitudes and in the Tibetan Plateau in both seasons (Fig. 3). To have a better overview of how the skill scores changed, we summarize their distributions with the help of box plots. Figure 4 shows the differences of skill scores between the aniso experiment and the original skill for the three variables (temperature, precipitation and sea-level pressure) in the ENH region. In the instrumental-only experiment, correlation values of temperature and sea-level pressure decreased in both seasons, while for precipitation they remained mostly unchanged. The RE values show that the experiments with anisotropic localization function reduced the skill of the reconstructions, but the extent of the reduction varies with the variables and with the seasons (Fig. 4). In general, the same holds for the proxy-only experiment (Fig. 4).

Figure 3Spatial skill of temperature reconstruction presented by RE values, assimilating only instrumental data (a, b, d, e) and only proxy records (c, f). Comparing the skill of the reconstruction using isotropic localization function (a–c) versus an anisotropic localization function (d–f). Skill in the winter season (a, d) and in the summer season (b, e, c, f) are shown.

Figure 4Difference between the aniso experiment and the original setup in terms of skill scores over the ENH region. Distributions of correlation values and of RE values are on the left and right figures, respectively. Distribution of temperature (a, b), precipitation (c, d) and sea-level pressure (e, f) are shown. Blue color indicates the instrumental-only experiment and yellow indicates the proxy-only experiment. The midline of the box is the median. The lower (upper) border of the box is the first (third) quartile. The whiskers extend up to 1.5 times the interquartile range; beyond these distances, the number of outliers is given under the box plots. The grid boxes were not area weighted. The asterisk above the box indicates significant differences between the median of the experiment and the original setup.

Download

4.2 Inflation experiments

4.2.1 Results

The main problem of ensemble-based DA techniques is the computationally affordable limited ensemble size. Due to the finite ensemble size, the estimation of P^b suffers from sampling error. Applying inflation techniques is one method to mitigate its effect (see Sect. 3.2).

Figure 5Distribution of correlation coefficients differences between the mixed background-error covariance matrix experiments and the original setup over the ENH region. Panels (a), (d) and (g) show the skill of the reconstruction in the winter seasons, while panels (b), (c), (e), (f), (h) and (i) show it for the summer season. The labels on the x axis indicate the experiments. Box plot, color, number and asterisk are the same as in Fig. 4.

Download

Using the multiplicative inflation method, the deviations from the ensemble mean are multiplied with a small factor (γ). To find the optimal γ, a set of experiment runs is required. We used γ=1.02 and γ=1.12 in our experiments, where only instrumental data were assimilated. We chose γ from a range that was previously tested by Whitaker and Hamill (2002). Multiplying the deviations from the ensemble mean with γ=1.02 in the assimilation process hardly affected the skill of the reconstruction over the ENH region (not shown). When we increased the value of γ to 1.12, the RE values slightly decreased (not shown). We did not carry out further experiments since, based on the results, randomly increasing the error in background field did not lead to improvement.

In the other set of experiments, we used P^blend in the update equation (Eq. 6). The experiments were run with using β₂ equal to 0.25, 0.50, 0.75 and 1 to estimate the P^blend (denoted 25c, 50c, etc.). Besides the varying weight given to P^clim, the applied L values on P^b and P^clim differed as well. Three L values were used: no localization (termed “no”), applying L values as in Table 1 (L) and doubling these numbers (2L). Different combinations of the fraction of P^clim and L values were termed accordingly (e.g., 50c_PbL_Pc2L).

Figure 6Distribution of RE value differences between the mixed background-error covariance matrix experiments and the original setup over the ENH region; otherwise, it is the same as in Fig. 5.

Download

We expect that estimating the covariances from a bigger ensemble size (n = 100–500) instead of 30 members leads to a more accurate background matrix. In most of our experiment, n is 250. Hence, P^clim is likely less affected by the sampling error implying that long-range spurious correlations are less prominent, which makes localization less needed. We presume that using P^blend helps to better reconstruct areas which were characterized with lower skill score values in the original setup and to improve the estimation of unobserved climate variables. The reconstruction skill of the blending experiments is always calculated from x^a (Fig. 2).

For the ENH region, we present how the verification measures changed by replacing P^b with P^blend in the assimilation process. We conducted an experiment without localizing P^clim and using L values from Table 1 on P^b in the construction of P^blend. However, the skill of the reconstruction was largely reduced, implying that 250 members are not enough to avoid localization altogether (not shown).

Figures 5 and 6 show the distribution of the differences of the skill scores between the experiments and the original analysis for correlation coefficients and RE values, respectively. Depending on the variables and the data type being assimilated, different setups perform better. In the case of assimilating only instrumental data, one of the largest increases of median for temperature reconstruction was obtained from the 100c_PcL experiment in both seasons (Figs. 5, 6). Precipitation records were not assimilated; thus, a reasonable estimation of the cross-variable covariances is essential. The skill of the precipitation reconstruction with the original setup, in terms of correlation, is better than the forced simulation (not shown); however, the RE values are negative over large regions in the ENH (Fig. 7). Using P^blend in the assimilation, with, e.g., the settings of 50c_PbL_Pc2L experiment, led to more positive RE skill (Fig. 7). The biggest improvement, in terms of RE skill score, was found in Europe (Fig. 7). The 50c_PbL_Pc2L analysis also has higher skill in North-America, especially in the summer season (Fig. 7). The skill of the sea-level pressure reconstruction also improved in the 50c_PbL_Pc2L experiment (Figs. 5, 6). In the proxy-only experiments, 75c_PbL_Pc2L is among the best performing experiments for all the variables (Figs. 5, 6).

Figure 7Spatial reconstruction skill of precipitation in terms of RE values, assimilating only instrumental data. Panels (a) and (b) show the skill of the original setup, and (c) and (d) show the result of the 50c_PbL_Pc2L experiment. The skill in the winter season is presented in panels (a) and (c), and for the summer season in panels (b) and (d).

We also investigated the effect of the ensemble size in the estimation of P^clim. To test whether further improvements can be achieved by doubling the ensemble size of x^clim, we ran an experiment with the following setup: β₁ and β₂ are equally weighted, and L and 2L are applied on P^b and P^clim, respectively (Table 2). In the experiment, we assimilated only instrumental data. The skill scores of x^a (corr, RE) from the 500-ensemble-member experiment showed no marked improvement compared with the same experiment with 250 ensemble members. An additional experiment was carried out with the same setup but using only 100 ensemble members in the construction of x^clim. The verification measures of the 50c_PbL_Pc2L_100m experiment are higher than the original one, and the distribution of the skill scores over the ENH region is very similar to what we obtain by using 250 members in P^clim for temperature and precipitation. However, the sea-level pressure fields from the 50c_PbL_Pc2L experiment have higher skill than those in the 50c_PbL_Pc2L_100m experiment (not shown).

Furthermore, we conducted two experiments in which only x^b was updated after an observation was assimilated, and x^clim was kept constant in the assimilation window. However, the ensemble members of x^clim were randomly reselected for each year (October–September). The advantage of this setup compared to the setup described in Sect. 3.2 is that it is computationally less demanding since only the original 30 members keep being updated with the observations. In the first test, we give β₂=0.75 weight to P^clim with 2L values. In the second test, β₂=1; that is, only P^clim used for updating x^b and the L values in Table 1 were applied for the localization. By comparing the skill of the reconstructions without and with updating the climatological part, we see that the skill scores are higher when the climatological part is also updated with the information from the observations (Fig. 8). The only exception is the correlation values of sea-level pressure: when keeping the climatological part constant, they are slightly higher in both seasons (Fig. 8). Nonetheless, by keeping the climatological part static in one assimilation window, the experiments still outperform the original reconstruction (Fig. 8).

Figure 8Distribution of skill scores over the ENH region. The skill of the original setup is compared with experiment 75c_PbL_constPc2L, 75c_PbL_Pc2L, 100c_constPcL and 100c_PcL. Distribution of correlation coefficients in the winter (left column) and in the summer (right column) seasons. Distribution of RE values in the winter (left column) and in the summer (right column) seasons.

Download

4.2.2 Discussion

We have tested a number of configurations of the mixed covariance matrix P^blend to evaluate the effect of the sampling error. In numerical weather predication (NWP) applications, various methods have been designed to better estimate the errors of the background state. In hybrid DA systems, the advantages of variational and ensemble Kalman filter techniques are combined (Hamill and Snyder, 2000; Lorenc, 2003). In another method, the background-error covariances are obtained from an ensemble of assimilation experiments performed by a variational assimilation system (Pereira and Berre, 2006). In an additive inflation experiment, a term is added to the x^a to account for the errors of the DA system (Whitaker et al., 2008).

Figure 9Difference of the RE skill between the temporally localized experiment and the original setup, when only instrumental data are assimilated. Temperature (a, b) and precipitation (c, d) differences are shown in the winter (a, c) and in the summer (b, d) seasons. The black dots indicate the negative RE values in the temporally localized experiment.

In our implementation, P^blend is calculated from x^b and x^clim. Using P^blend in the assimilation process improved on the reconstruction performed with the original setup. The skill scores show the largest improvement in the sea-level pressure reconstruction. Moreover, the skill of the precipitation reconstruction also improved, indicating that P^clim helps to better estimate the cross-covariances of the background errors between the variables. In general, increasing the weight of P^clim in forming P^blend positively affected the skill of the analysis. The 100c_PcL experiment, in which P^blend is equal to P^clim, is similar to the DA technique used in the last millennium climate reanalysis (LMR) project (Hakim et al., 2016). In the LMR, 100 randomly chosen ensemble members form a climatological state vector, which is used in each assimilation window and is updated with the observations. In this study, x^clim is randomly resampled every year and primarily used in the estimation of P^blend. The settings used in the 100c_PcL experiment led to one of the largest increases in the median for temperature reconstruction when only instrumental measurements are assimilated. However, other settings resulted in larger increase of median for different variables and observation types. By applying no localization on P^clim in the 50c_PbL_PcnoL experiment, we obtained a less skillful reconstruction than by using the other two localization schemes. The skills reduced especially over the areas where no local observations were assimilated. Using 2L values for localizing the covariances of P^clim in the instrumental-only experiments resulted in higher correlation values of sea-level pressure (50c_PbL_Pc2L) and helped to obtain higher correlation scores of precipitation in summer. Among the proxy-only experiments, 75c_PbL_Pc2L shows the largest increase of median for pressure reconstruction. Here, pressure data are not assimilated, and the result suggests that by applying longer L values, the cross-variable covariances are treated better. We tested whether the skill of the experiments performed with various settings is significantly different from the skill of the original analysis. We compared the median value of the skill scores from the experiments and the original data, and with most of the settings a significant difference was obtained for all the variables. The results of the experiments show that with a mixed covariance matrix implementation, a major drawback of the ensemble-based DA system, due to the limited ensemble size, can be improved.

4.3 Localization in time