Assimilating monthly precipitation data in a paleoclimate data assimilation framework

Data assimilation approaches such as the ensemble Kalman filter method have become an important technique for paleoclimatological reconstructions and reanalysis. Different sources of information, from proxy records and documentary data to instrumental measurements, were assimilated in previous studies to reconstruct past climate fields. However, precipitation reconstructions are often based on indirect sources (e.g., proxy records). Assimilating precipitation measurements is a challenging task because they have high uncertainties, often represent only a small region, and generally do not follow a Gaussian distribution. In this paper, experiments are conducted to test the possibility of using information about precipitation in climate reconstruction with monthly resolution by assimilating monthly instrumental precipitation amounts or the number of wet days per month, solely or in addition to other climate variables such as temperature and sea-level pressure, into an ensemble of climate model simulations. The skill of all variables (temperature, precipitation, sea-level pressure) improved over the pure model simulations when only monthly precipitation amounts were assimilated. Assimilating the number of wet days resulted in similar or better skill compared to assimilating the precipitation amount. The experiments with different types of instrumental observations being assimilated indicate that precipitation data can be useful, particularly if no other variable is available from a given region. Overall the experiments show promising results because with the assimilation of precipitation information a new data source can be exploited for climate reconstructions. The wet day records can become an especially important data source in future climate reconstructions because many existing records date several centuries back in time and are not limited by the availability of meteorological instruments.

1 Introduction

Precipitation is one of the key components of the climate system. Understanding its variability is fundamental due to its impact on the ecosystem and on human society. Instrumental observations are the main data source for studying its spatiotemporal variability. However, instrumental measurements are often insufficient because long time series are rarely available and their spatial distribution is rather sparse, especially further back in time. To examine the decadal variability of precipitation, longer time series are needed. This is often done by analyzing proxy records, documentary data, or model simulations.

Simulations suggest, for instance, that the tropical monsoon regions are characterized by the largest interannual variability of annual precipitation, and the interannual variability exhibits significant changes on the multi-decadal scale (Yang and Jiang, 2015). Precipitation reconstructions are needed to validate such results. Their number increased in recent years on all scales from local to large scale. To reconstruct local millennia-long hydroclimate variability, tree-ring series were used, for example, in southern–central England (Wilson et al., 2013) and in southern Scandinavia (Seftigen et al., 2017). Pauling et al. (2006) reconstructed a 500-year-long seasonal precipitation field over Europe back to the 16th century by using instrumental measurements, documentary data, and proxy records. A reconstruction of Northern Hemisphere hydroclimate variability from multi-proxy records and documentary data is available between the 9th and 20th century (Ljungqvist et al., 2016). A similar reconstruction was also produced for southern South America for the last 500 years (Neukom et al., 2010). Centuries-long tree-ring drought atlases are available for North America (Cook et al., 2010b), Asia (Cook et al., 2010a), Europe (Cook et al., 2015), and eastern Australia and New Zealand (Palmer et al., 2015) to study long-term hydroclimate variability. Steiger et al. (2018) produced the first global hydroclimatic reconstructions at annual and seasonal resolutions by combining multi-proxy data with the Community Earth System Model Last Millennium Ensemble model simulations (Otto-Bliesner et al., 2016) over the last 2 millennia. A multi-century global reconstruction making use of observational precipitation data is still missing.

One challenge in global reconstructions is the required observation network density. To test how dense a network for a climate field reconstruction needs to be, pseudo-proxy experiments can be performed. A set of such experiments was conducted by Gomez-Navarro et al. (2015) to reconstruct precipitation over Europe with three climate field reconstruction methods (canonical correlation analysis, analog method, Bayesian hierarchical method). They found that the skill of the precipitation reconstruction increases close to the proxy locations (correlations >0.8) independently of the method used. However, within a few hundred kilometers correlation values drop below 0.2 (with seasonal dependency), implying that precipitation is no longer accurately reconstructed at these distances. Therefore, they conclude that a dense network of proxy records, in accordance with the high spatiotemporal variability of precipitation, is essential for successful reconstructions.

Thanks to the introduction of new techniques in the field of paleoclimatology, nowadays spatially complete and physically consistent reconstruction can be derived. Paleoclimate data assimilation provides a framework in which observational data and model simulations are optimally combined to obtain global, three-dimensional climate fields. Using the data assimilation approach is advantageous because by assimilating only one type of climatic variable (e.g., temperature), we can gain information from other climatic fields present in the model simulation based on their covariances between the observed and unobserved variables. Monthly instrumental temperature and sea-level pressure data, documentary temperature data, and tree-ring proxy records have been successfully used with the ensemble Kalman fitting (EKF) method to reconstruct past climate fields (Franke et al., 2017a). In the EKF400 atmospheric paleo-reanalysis, skillful temperature and sea-level pressure reconstruction can be expected a few thousand kilometers away from observations, whereas precipitation fields show limited skill (Franke et al., 2017a). The reasons are the high spatial heterogeneity of precipitation and the fact that no precipitation data have been assimilated.

In this paper, monthly precipitation information in the form of precipitation amounts or the number of days with precipitation in each month (wet days) is assimilated with the EKF method (Franke et al., 2017a) to reconstruct monthly climate fields. Early instrumental precipitation measurements are available since the 17th century (New et al., 2001). Descriptive weather journals were kept by scholars even before the instrumental era, usually including the number of wet days. Examples of such weather diaries exist, for example, for Kassel (Germany) by Uranophilus Cyriandrus covering 1621–1650 (Lenke, 1960), for Zürich (Switzerland) by Wolfgang Haller between 1545 and 1576 (Pfister et al., 1999), and for Savanna-la-Mar (Jamaica) by Thomas Thistlewood between 1750 and 1786 (Chenoweth, 2003). In future reconstructions, these documentary records and numerous others could add valuable information to the limited instrumental measurement data to assess the natural variability of precipitation. Here, we investigate the possibility of using these new observation sources in a climate context by conducting and evaluating a set of experiments over the 1950–2004 time period. The skill of the reconstructions with the assimilation of precipitation amount versus wet day records is compared. The precipitation amount and wet day records are assimilated together with other instrumental measurements to study their combined effect.

This paper is structured as follows: in Sect. 2 an overview of the EKF techniques is given and the model simulation and the observations are introduced. Section 3 describes the experimental framework and the skill metrics used for evaluation. In Sect. 4 the results are presented, and they are discussed in Sect. 5. We conclude how monthly precipitation information can be assimilated best in Sect. 6.

2 Methods

2.1 Ensemble Kalman fitting

To reconstruct climate fields by assimilating monthly precipitation amounts and the number of wet days per month, the ensemble Kalman fitting (EKF) data assimilation technique was used (described in detail in Bhend et al., 2012, and Franke et al., 2017a). EKF is used in an offline manner; that is, only the update step of the Kalman filter is implemented. In the update step of the EKF precomputed model states (see Sect. 2.2) are updated with observational information as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{b}}+\mathbf{K}\left(\mathit{y}-\mathbf{H}{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{b}}\right),\text{(2)}& {\displaystyle }{{\mathit{x}}^{\prime }}^{\mathrm{a}}={{\mathit{x}}^{\prime }}^{\mathrm{b}}+\stackrel{\mathrm{̃}}{\mathbf{K}}\left(\mathbf{H}{{\mathit{x}}^{\prime }}^{\mathrm{b}}\right),\end{array}$$

where the analysis ensemble mean (${\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}}$) and the deviation from it (${{\mathit{x}}^{\prime }}^{\mathrm{a}}$) are updated separately. x^b is the background state vector that is provided by the model simulation. H is the forward operator, by which the model value is transformed to observation value. In our setup H is linear. K and $\stackrel{\mathrm{̃}}{\mathbf{K}}$ denote the Kalman gain matrix and the reduced Kalman gain matrix, which weight the model and observation information based on their error estimates:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathbf{K}={\mathbf{P}}^{\mathrm{b}}{\mathbf{H}}^T{\left({\mathbf{HP}}^{\mathrm{b}}{\mathbf{H}}^T+\mathbf{R}\right)}^{-\mathrm{1}}\text{(4)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{̃}}{\mathbf{K}}={\mathbf{P}}^{\mathrm{b}}{\mathbf{H}}^T& {\left({\left(\sqrt{{\mathbf{HP}}^{\mathrm{b}}{\mathbf{H}}^T+\mathbf{R}}\right)}^{-\mathrm{1}}\right)}^T,\\ & \times {\left(\sqrt{{\mathbf{HP}}^{\mathrm{b}}{\mathbf{H}}^T+\mathbf{R}}+\sqrt{\mathbf{R}}\right)}^{-\mathrm{1}}\end{array},\end{array}$$

where P^b is the background-error covariance matrix and R is the observational-error covariance matrix. In ensemble-based Kalman filter methods, P^b is calculated from the ensemble members (Evensen, 2003). The estimation of the observational error is discussed later in Sect. 2.3.2. Gaussian statistics in the errors are assumed.

The EKF is implemented serially; that is, observations are assimilated one by one, which makes the assimilation process computationally more efficient. The localization of P^b is necessary if the covariances are calculated from a much smaller ensemble size than the model dimension (discussed among others in Kepert, 2009). Therefore, P^b was multiplied element-wise with a correlation function that decreases as distance increases. The correlation function was represented with the isotropic Gaussian function:

$$\begin{array}{}\text{(5)}& G=\exp \left(-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}{L}^{\mathrm{2}}}}\right),\end{array}$$

where z denotes the distance between two grid boxes, and L is the length scale parameter (Gaspari and Cohn, 1999). In most of the experiments of this study, we use the localization length scale parameters determined by Franke et al. (2017a). The localization length scale parameter for temperature is 1500 km, for pressure 2700 km, and for precipitation 450 km. In the case of the number of wet days, the localization length scale parameter is set to 450  km. Between two different variables the smaller L is always used for localization.

2.2 Model simulation

The Chemical Climate Change over the past 400 years (CCC400) model simulations serve as the model background. They were produced with the ECHAM4.5 model (Roeckner et al., 2003) and have 30 ensemble members (Bhend et al., 2012). CCC400 is a forced model simulation of the atmosphere that was performed at T63 horizontal resolution with 31 layers in the vertical (for more details see Bhend et al., 2012, and Franke et al., 2017a). CCC400 was used in previous climate studies (Franke et al., 2017a; Valler et al., 2019), but here an additional climate field (wet days) was added to the already existing monthly model fields. Hence, in this study the background state vector (x^b) has the following monthly climate fields: temperature, precipitation, sea-level pressure, and wet days. Previous climate model evaluation studies showed that precipitation characteristics such as frequency and intensity are often over- or underestimated in the model simulations (e.g., Sun et al., 2006; Koutroulis et al., 2016). The fraction of light precipitation is overestimated by the models; that is, it drizzles too frequently (Sun et al., 2006). Therefore, in this study the threshold for wet day was set to ≥1 mm d⁻¹. The precomputed model states cover the period between 1601 and 2004. However, in this study we mainly focus on the period between 1950 and 2004.

Figure 1 The location of precipitation observations (a) used in the exp_R, exp_W, exp_R_2L, and exp_W_2L experiments. The location of instrumental precipitation (blue), temperature (red), and sea-level pressure (green) observations (b) assimilated in the exp_TPR and exp_TPW experiments.

2.3 Observation

2.3.1 Observation dataset

Precipitation observation data are obtained from the Global Historical Climatology Network – Daily database (GHCN-Daily v3.22; Menne et al., 2012a). Before 1950, the number of available station data in GHCN-Daily drastically decreases in large parts of the world (Menne et al., 2012a). Therefore, we only use time series from the 1950–2005 period. Furthermore, to ensure homogeneous data completeness and spatial coverage, we kept only those stations that fulfill the following criteria: (1) at least 50 years of data are available; and (2) the distance between two stations is at least 200 km. After the selection process, 432 stations remained (Fig. 1a). In the extratropical Northern Hemisphere the station network is dense and equally distributed. However, in the Southern Hemisphere outside Australia, very few station data are available.

The daily precipitation sums were converted to monthly totals. For the conversion of daily precipitation sums into the number of wet days, a threshold of 1 mm was used (i. e., days with precipitation <1 mm were not considered wet days).

2.3.2 Representativity error

Rain gauge measurements are subject to systematic, random, and representativity errors (e.g., Lopez, 2013). Representativity error, in particular, arises when comparing point measurements to area averages such as grid points in model simulations. For precipitation, it is arguably the largest error source due to the high spatial variability of this variable.

To estimate the representativity error, we analyzed all GHCN-Daily station data over the contiguous United States and adjacent territories (lat: 24–48^∘ N, long: 126–68^∘ W) in the 1961–1990 period. This region has the required network density and covers several climatic zones, which makes it a good subset to estimate global errors in the monthly precipitation amount and number of wet days.

Figure 2 The distribution of the representativity error of precipitation (a) and the number of wet days (b). The red line indicates the median.

The representativity error was estimated using the following procedure for both precipitation amounts and wet days: (1) an average value for all months between 1961 and 1990 was calculated at each single grid box from the stations located within one grid box according to the resolution of the CCC400 model simulation. (2) The spatially averaged monthly time series of each grid box were subtracted from all station series within the grid box. (3) The standard deviation of the resulting series was calculated. (4) The median value of the standard deviation was taken as the representativity error. Figure 2 shows the distribution of the standard deviations. The medians of the standard deviation of precipitation and the number of wet days are 24.63 mm (equivalent to 28 % error) and 1.49 d (equivalent to 21 % error), respectively.

3 Experiments

3.1 Experiment design

An important feature of our assimilation process is that only anomalies are assimilated. Both model and observation data are transformed into anomalies using a moving 70-year time window around the current year (this window is shorter at the edges of the available period, from which the anomalies are calculated). Working with anomalies alleviates the problem of model biases (see Franke et al., 2017a).

If more than one station is available within the same grid box, then the average of the observation anomalies was assimilated. As in previous studies (e.g., Franke et al., 2017a) the length of the assimilation window is 6 months. Therefore, in x^b six monthly fields of different variables are combined. The assimilation windows cover the period between October–March and April–September, which were chosen due to the growing seasons of trees if tree-ring proxy data are assimilated. We eventually aim to assimilate precipitation data in combination with proxy data. Hence, the half-yearly assimilation scheme was kept the same.

When data for 6 months are combined together in x^b for one assimilation step, the assimilated observations can influence all the variables of the 6 months. However, Valler et al. (2019) found that in the case of multiple months combined into one state vector, the analysis improved when time localization was applied; that is, observations can only affect the current month. Hence, time localization was applied in all experiments described in this paper.

Table 1 Summary of experiments and their setups.

A set of experiments was conducted, which are summarized in Table 1. In our first experiment we analyze the skill of the reconstruction when assimilating monthly instrumental precipitation amount in the period 1950–2004 (exp_R). In the second experiment we assimilate the monthly number of wet days instead of the precipitation amount (exp_W). Both experiments were repeated with an extended spatial localization distance in order to make better use of the available observations (exp_R_2L and exp_W_2L). The length scale parameter for both precipitation amount and wet days was doubled from 450 to 900 km based on the calculation of the decorrelation distance of precipitation in the model space. In two further experiments, we added instrumental temperature and pressure data to the assimilated precipitation amount or to the wet day records (exp_TPR and exp_TPW, respectively) to study the combined influences of various variables (Fig. 1b). Due to the localization of P^b, the sequence of assimilated observations can affect the analysis (x^a) (Greybush et al., 2011). We kept the same assimilation order used in a previous reconstruction (Franke et al., 2017): assimilating temperature data first then pressure measurements. Precipitation data were assimilated last due to the bigger uncertainties characterizing precipitation measurements. The assimilated temperature and pressure time series are described in Franke et al. (2017a). In all experiments, we used an observational error of 30 % or at least 10 mm for precipitation amount and of 2 d for wet days. These values were chosen in agreement with our estimation of the representativity error (Sect. 2.3.2). For temperature and pressure we applied the observational errors of $\sqrt{\mathrm{0.9}}$ K and $\sqrt{\mathrm{10}}$ hPa, respectively (Valler et al., 2019).

Finally, an additional experiment was conducted to reconstruct a severe drought event in Europe in order to demonstrate the potential of assimilating precipitation data. Six available stations from GHCN-Daily were used to reconstruct the extreme drought year of 1842 (Brázdil et al., 2019, and references therein). Brázdil et al. (2019) collected several documentary data sources, from diaries and newspapers to scientific papers, to analyze the severity of the drought. Besides documentary data, instrumental measurements were also examined. The year 1842 was characterized by low water levels in many countries from the Netherlands to western Ukraine and from Italy to Sweden, affecting shipping, the availability of drinking water, and agriculture. We looked at how this drought event is captured in an independent reconstruction using the Twentieth Century Reanalysis version 3 (20CRv3; Slivinski et al., 2019a). In 20CRv3 only pressure measurements are assimilated into a model with prescribed sea surface temperatures, sea ice concentrations, and radiative forcings. The 20CRv3 reanalysis has 80 ensemble members that provide an idea of the confidence in the reanalysis at any time and location. For 1842 data are available from 39 distinct locations, where the reanalysis is more reliable compared to regions with no available observations. The monthly precipitation anomaly fields of 20CRv3 are calculated from the 1836–1877 time period. The presented results compare the 20CRv3 ensemble mean with the ensemble mean of the experiments described above. In addition, we compared our exp_R experiment with a seasonal high-resolution precipitation reconstruction over Europe (Pauling et al., 2006). Pauling et al. (2006) used various data types (instrumental, documentary, proxy) and principal component regression techniques to reconstruct precipitation between 1500 and 2000. From their dataset, we calculated the seasonal relative precipitation anomalies of 1842 by dividing the absolute anomaly with the 71-year climatology centered on 1842. Relative precipitation anomalies were used for monthly field comparison as well.

3.2 Evaluation

In order to evaluate the skill of the sensitivity experiments, two commonly used skill metrics – correlation and the root mean squared error skill score (RMSESS) – are calculated over the 1950–2004 time period. The CRU TS3.10 dataset (Harris et al., 2014) was chosen as a reference dataset to evaluate the precipitation, temperature, and sea-level pressure reconstructions. CRU TS3.10 is a gridded dataset, which is relaxed to the 1961–1990 climatology in the case of no available observations. Correlation is calculated between the absolute values of the ensemble mean of the analysis and the reference series; it expresses the covariability between the two series but does not measure biases. In Sect. 4.1, we show the differences between the correlations calculated from the analysis–reference data and from the model–reference data. We show the correlation differences because the CCC400 model simulations are transient, forcing-dependent simulations and already have skill. The RMSESS metric is calculated as

$$\begin{array}{}\text{(6)}& \mathrm{RMSESS}=\mathrm{1}-{\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n({\mathit{x}}_i^u-{\mathit{x}}_i^{\mathrm{ref}}{)}^{\mathrm{2}}}}{\sqrt{\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n({\mathit{x}}_i^f-{\mathit{x}}_i^{\mathrm{ref}}{)}^{\mathrm{2}}}}},\end{array}$$

where x^u is the reconstruction, x^f is the model simulation, x^ref is the reference dataset, and i refers to the time step. In experiments when only the precipitation amount or the number of wet days is assimilated, the RMSESS is calculated on the anomalies from the ensemble mean of the analysis and the ensemble mean of the model simulation and compared to the reference dataset. In the case of the exp_TPR and exp_TPW experiments, x^f is replaced with the ensemble mean of the analysis from the exp_TP experiment. The two skill metrics measure different properties, and therefore the skill of the reconstruction also depends on the analyzed skill metrics (Franke et al., 2020). When, for example, only the number of wet days is assimilated, then the temperature reconstruction is fully independent of the reference dataset. However, in other cases the skill of the reconstruction can be overestimated (see also Franke et al., 2020). Monthly skill metrics are discussed for the exp_R experiment, while for the other experiments the seasonal skill metrics – calculated from seasonal averages (winter: October–March, summer: April–September) – are analyzed. Both seasons are discussed later, but only the results of summer are shown in the paper (figures for winter are shown in the Supplement).

Figure 3 Summer season correlation differences of precipitation (a, b, c, d), temperature (e, f, g, h), and sea-level pressure (i, j k, l) between the analyses and the CCC400 model simulation ensemble means by assimilating only precipitation amounts (a, e, i), only precipitation amounts with a doubled localization length scale parameter (b, f, j), only the number of wet days (c, g, k), and only the number of wet days with a doubled localization length scale parameter (d, h, l).

4 Results

4.1 Skill score analysis

Assimilating monthly precipitation amount data (exp_R) led to improved precipitation skill compared to the existing correlation between the CCC400 model simulation ensemble mean and the reference dataset. The monthly correlation differences of precipitation show clear improvement (winter: Fig. S1 in the Supplement, summer: Fig. S2). The correlation differences between the analysis and the model simulation are almost always positive in the case of temperature and sea-level pressure in all months (Figs. S1, S2). In terms of RMSESS values, during boreal winter the monthly precipitation fields show reduced skill in the high northern latitudes and in Siberia (Fig. S1). The skill of the precipitation reconstruction gradually decreases from October to March over Siberia (Fig. S1). The negative RMSESS values remain present in April and to a lesser extent in May over northern North America and Siberia, while the RMSESS values are mainly positive throughout June and September in the Northern Hemisphere (Fig. S2). The precipitation reconstruction has mixed skill over Australia. Mostly negative RMSESS values characterize the northern and northwestern regions, while positive skill is seen in the southern and eastern parts (Fig. S2). The RMSESS values of the winter monthly temperature fields are in general positive, except Australia (Fig. S1). However, the skill of the temperature reconstruction is negative in large parts of North America and Europe, especially from June to August (Fig. S2). In the winter sea-level pressure fields, an improvement can be seen mainly over Europe and Asia (Fig. S1), while in the summer months the effect of precipitation on the pressure fields is mixed (Fig. S2). The winter seasonal skills of exp_R are shown in the Supplement (Figs. S3 and S4), while summer seasonal skills are shown in Figs. 3 and 4.

Figure 4 Summer season RMSESS values of precipitation (a, b, c, d), temperature (e, f, g, h), and sea-level pressure (i, j k, l) by assimilating only precipitation amounts (a, e, i), only precipitation amounts with a doubled localization length scale parameter (b, f, j), only the number of wet days (c, g, k), and only the number of wet days with a doubled localization length scale parameter (d, h, l).

The assimilation of wet days (exp_W) resulted in mainly positive correlation differences of all three variables (Figs. 3, S3). The correlation skill of the exp_W analysis is very similar to the skill of exp_R, in which precipitation amounts were assimilated. However, the skill of the reconstructions shows a different picture when the RMSESS values are analyzed. The precipitation field shows improved skill almost everywhere in summer (Fig. 4c), and in winter the RMSESS values are negative only in the high northern latitudes and over the west and central part of northern Asia (Fig. S4). The negative skill over Europe seen in the summer temperature field when assimilating precipitation amounts (Fig. 4e) is diminished by using the number of wet days (Fig. 4g). However, the negative RMSESS values over the central Siberian region seen in summer in the exp_R experiment did not change with the assimilation of wet days (Fig. 4d). In general, the assimilation of wet days (exp_W) has higher skill than the assimilation of precipitation amounts (exp_R) when temperature fields are compared. The positive effect of assimilating the number of wet days instead of the precipitation amount is also evident in the RMSESS values of the sea-level pressure fields (Figs. 4k, S4). The sea-level pressure analysis of exp_W has positive skill over Europe in both seasons.

Doubling the correlation length scale parameter of precipitation resulted in a very similar correlation skill of precipitation (Figs. 3b, S3) as in the exp_R experiment. Increasing the correlation length scale parameter mainly positively affected the sea-level pressure reconstruction (Fig. 3j). The same holds for the temperature correlation (Fig. 3f). While the increased localization distance positively affected the correlation of the reconstructed fields, it had a more mixed impact on the other skill metric. The RMSESS values of precipitation decreased in the high northern latitudes, over North America, and in the Mediterranean region (Figs. 4b, S4). This negative effect of the increased localization length scale is also seen in the sea-level pressure, primarily in summer (Fig. 4j). The temperature reconstruction is affected differently in the two seasons. In winter the RMSESS mainly increased in the Northern Hemisphere (Fig. S4), while in summer larger areas are affected negatively (Fig. 4f).

The same experiment with a doubled localization length scale parameter was conducted with the assimilation of the number of wet days. Correlation coefficients of precipitation in the exp_W_2L experiment remained high (Figs. 3d, S3), similar to what was obtained with a more strict localization scheme in the exp_W experiment. Both sea-level pressure and temperature reconstructions benefited from the larger localization length scale, except a few grid boxes (Figs. 3h, l; S3). The precipitation reconstruction shows similar skill to the exp_R experiment in terms of RMSESS values. However, the increase in the localization length scale had a positive impact on the sea-level pressure fields over Europe in both seasons (Figs. 4l, S4).The regions over central South America and southern Africa that were negatively affected in the exp_W experiment show worsening skill in the exp_W_2L experiment (Fig. S4). Similarly to the sea-level pressure field, in some regions the temperature reconstruction became more skillful, while in others the skill decreased (Figs. 4h, S4).

The distributions of the skill matrices of the presented experiments over the extratropical Northern Hemisphere are summarized as box plots (Fig. S5). Increasing the localization length scale parameter positively affected the reconstructed temperature and sea-level pressure fields in terms of correlation, especially in boreal winter. The sea-level pressure and temperature RMSESS values are less affected by the applied localization. However, the RMSESS of the precipitation reconstruction decreased in the experiments with a doubled localization length if precipitation amounts were assimilated.

Figure 5 Summer season correlation differences between the analysis of exp_TP and the CCC400 model simulation ensemble means: precipitation (a), temperature (d), sea-level pressure (g). In the middle column (b, e, h) the correlation differences between the ensemble mean of exp_TPR and exp_TP are shown, while in the right column (c, f, i) the correlation differences between the exp_TPR and exp_TP ensemble mean are depicted.

Figure 6 Summer season RMSESS values of exp_TP (a, d, g), exp_TPR (b, e, h), and exp_TPW (c, f, i). x^f is the ensemble mean of the CCC400 model simulations when the RMSESS values are calculated for the exp_TP experiment, while x^f was replaced with the analysis mean of exp_TP for the other two experiments.

In the next experiment different observation types such as temperature, sea-level pressure, and precipitation amounts were combined (Fig. 1b). In previous studies, monthly instrumental temperature and sea-level pressure data were already successfully assimilated (Franke et al., 2017a; Valler et al., 2019). Here we added precipitation measurements to the assimilated variables. To see the effect of assimilating precipitation data an experiment using only temperature and sea-level pressure data was carried out over the 1950–2004 time period (exp_TP), which resembles the original setup in Franke et al. (2017a) except that it does not assimilate any proxy information. In the exp_TPR experiment, first the temperature, then the sea-level pressure data, and finally the precipitation amounts were assimilated. As mentioned above, in this experiment the skill of the exp_TPR analysis is compared to the exp_TP analysis mean. Precipitation correlations show a marked improvement over exp_TP in both seasons (Figs. 5b, S6), as expected. Temperature and sea-level pressure correlations are not changed in most parts of the world (Figs. 5e, h; S6). The positive influence of assimilating the precipitation amount on the precipitation reconstruction is not only seen in the correlation but also in the RMSESS values. The exp_TP experiment has positive RMSESS values, mainly over Europe, in both seasons (Figs. 6a, S7). Adding a precipitation amount to the assimilated observations (exp_TPR) increased the precipitation reconstruction RMSESS in winter, for example, across Europe to the Urals (Fig. S6), while in summer the skill tends to be more positive from the Urals to northern and eastern Asia (Fig. 6b). The temperature reconstruction in terms of RMSESS values mainly shows slight improvement in winter (Fig. S7), but in summer the skill is negatively affected by the assimilated precipitation amount in North America and Europe (Fig. 6e). The impact of the precipitation amount on the sea-level pressure reconstruction is mixed. In Europe many pressure records are available, and it is over Europe where the skill decreases the most in both seasons (Figs. 6h, S7). Besides Europe, North America is also negatively affected, especially in summer (Fig. 6h). The experiment was repeated with the number of wet days added instead of precipitation amounts (exp_TPW), and similarly as a reference the exp_TP experiment was used. The correlation differences of all three variables are very similar to the differences seen between exp_TPR and exp_TP (Figs. 5, S6). The precipitation RMSESS values also show a similar pattern to the exp_TPR experiment, but, for instance, over Australia in boreal winter the skill tends to be higher (Fig. S7). The summer season temperature reconstruction of exp_TPW is in some regions more skillful than the exp_TP experiment (Fig. 6). The sea-level pressure fields of exp_TPW depict negative RMSESS over Europe (Figs. 6i, S7). To further assess the impacts of assimilating precipitation data in regions with dense observations such as Europe, the distribution of skill metrics was analyzed. The reconstructed temperature fields of the three experiments have very similar skill over Europe (Fig. S8). Over Europe a loss of skill in sea-level pressure and temperature is seen for the RMSESS in both seasons, especially in summer, when the skill of exp_TPR and exp_TPW are compared to exp_TP (Fig. 6). If the RMSESS is calculated for all three experiments using the model simulation as x^f, the distributions over Europe indicate that the skill loss in sea-level pressure is smaller than the gain in precipitation for exp_TPR and exp_TPW (Fig. S8).

Figure 7 Relative precipitation anomalies over Europe in 1842. Monthly relative precipitation anomaly fields between April and September from exp_TP (a), exp_R (b), exp_W (c), and 20CRv3 (d) are depicted. Note that the anomalies are calculated from the 1807–1877 period for exp_TP, exp_R, and exp_W, as well as for 1836–1877 for 20CRv3. The red circles indicate the location of the time series used as input in the assimilations.

4.2 Case study: 1842 drought

Six stations in the GHCN-Daily dataset in Europe fulfilled the requirement to have continuous data in 1842. These stations, with their starting year in parentheses, are the following: Prague (1804), Jena (1826), Bologna (1813), Genoa (1833), Mantova (1849), and Armagh (1835). Using all available data in the 71-year time window centered on 1842, the precipitation amounts and the number of wet days were converted to anomalies. We calculated the relative precipitation anomalies for the months between April and September from exp_TP, exp_R, exp_W, and 20CRv3.

In the exp_R experiment (when only the precipitation amount was assimilated) the period between April and August in central Europe is characterized by negative relative precipitation anomalies (Fig. 7b). The experiment was repeated by assimilating the number of wet days (exp_W) instead of the precipitation amounts. The monthly relative precipitation anomaly fields are similar to the exp_R precipitation fields, but the deviations from the climatology are smaller (Fig. 7c). However, in May the south of France and northern Italy are wetter than the climatology in the wet day experiment (Fig. 7c), while in the exp_R experiment these regions are drier and closer to the climatology (Fig. 7b). We also investigated how this drought event is represented in the exp_TP experiment when no direct precipitation data were assimilated. The relative anomalies of precipitation indicate a dry period over central Europe from June to August, with the largest precipitation deficit in August (Fig. 7a).

These monthly relative anomaly fields were compared to the 20CRv3 reanalysis. As mentioned above, no precipitation data were used in 20CRv3. The relative precipitation anomaly fields of 20CRv3 and exp_R show similar precipitation patterns over the central region from France to western Ukraine. A precipitation deficit from April to August and a wet September are reproduced in 20CRv3 (Fig. 7d). The 20CRv3 reanalysis ensemble mean was remapped to the resolution of EKF400 in order to present the differences of relative precipitation anomalies between the experiments in this study and 20CRv3 (Fig. S9). The differences between the experiments and 20CRv3 are all negative in August and September over central Europe, while in the other months generally wetter conditions are present in exp_TP; exp_R and exp_W are in closer agreement with the reanalysis.

The drought over Europe in exp_R stands out even better on a seasonal timescale (June–July–August) (Fig. S10). Except northern Europe, Spain, and the southern part of France and Italy, negative relative precipitation anomalies define Europe (Fig. S10). Despite assimilating precipitation amount data from only six stations the obtained precipitation pattern is very similar to the reconstruction made by Pauling et al. (2006), which is based on several precipitation series, documentary data, and proxy records (Fig. S10). However, there are regional differences. In their reconstruction, a larger area with more negative anomalies is present over central Europe, while we find the strongest relative anomaly in exp_R over the northern part of central Europe. The sign of the precipitation anomaly differs over the eastern Mediterranean region, but in exp_R no data are used from this area.

Figure 8 Outcome of the Shapiro–Wilk test for normality for precipitation (a, b) and the number of wet days (c, d) in January (a, c) and July (b, d). The grey area shows where the null hypothesis of normality is rejected (p<0.05).

5 Discussion

5.1 Skill score analysis

In weather forecasting, there have been many attempts to make use of precipitation measurements from radar (e.g., Lopez, 2011), satellite (e.g., Lien et al., 2016), and rain gauges (e.g., Lopez, 2013). However, the usage of variables with a non-Gaussian distribution – such as precipitation – is a challenging task for many data assimilation methods; therefore, different transformation techniques have been tested (e.g., Lien et al., 2016). Monthly data have the advantage of being closer to a Gaussian distribution as a consequence of the central limit theorem. In the case of precipitation, however, this only applies to sufficiently wet climates in which several precipitation events occur each month.

Hence, the question remains as to whether the errors of precipitation amounts and the number of wet days are normally distributed, a fundamental assumption of EKF. In Fig. 8, we show the outcome of a Shapiro–Wilk test for normality (Wilks, 2011) applied to the deviations of the model ensemble members from the model ensemble mean, which represents our best estimation of the background error. On the one hand, we see that the distributions in arid climates and seasons are significantly different from a normal distribution: in these conditions the assimilation of precipitation by means of EKF cannot provide the optimal solution. On the other hand, large parts of the midlatitudes show errors close to a normal distribution, particularly in winter. The hypothesis of normality is less often rejected in the number of wet day data than in precipitation amount data, indicating that the former are more suitable for assimilation. This is probably one of the reasons why the skill of the reconstruction is much better in Europe when using wet days, although the climate of southern Europe is too dry in summer for either variable to be assimilated successfully.

Another advantage of the number of wet days over precipitation amounts is lower representativity error (see Sect. 2.3.2). The observation error of the number of wet days is also arguably lower because it is less affected by undercatch. Yang et al. (1999) showed that unshielded manual rain gauges underestimate solid precipitation by 10 % to 50 % on average, depending on location. The impact of undercatch on the number of wet days is much smaller, since it only affects days with precipitation between 1 and 2 mm (assuming a 50 % undercatch).

Assimilating monthly precipitation information such as precipitation amounts and the number of wet days, in general, shows positive improvements in all variables compared to the CCC400 model simulations when correlation is analyzed. The RMSESS of precipitation from both the exp_R and exp_W experiments is negative at the high northern latitudes and over large parts of Asia, especially in the winter season. As discussed earlier, precipitation observations are not error-free. Prein and Gobiet (2017) compared different gridded precipitation datasets over Europe and found large differences between them, comparable to the uncertainties seen between regional climate model simulations. They defined several possible sources of observational uncertainties, such as precipitation undercatch correction or station densities. The northern Asia region prior to 1957 shows the largest differences between different datasets in the Northern Hemisphere (Harris et al., 2014). Hence, it is difficult to evaluate the skill of the precipitation reconstruction from the exp_R and exp_W experiments in this region. Using an ensemble of observational datasets for the evaluation of the reconstruction could provide a better strategy (Prein and Gobiet, 2017). The RMSESS of temperature and sea-level pressure from the exp_R and exp_W experiments is not always positive, which suggests that the variance of the reconstructions is different from the chosen reference series. However, assimilating the number of wet days instead of the precipitation amounts greatly improved the RMSESS of both temperature and sea-level pressure reconstructions over Europe, especially in summer.

In addition to assimilating only precipitation amounts or only the number of wet days, the effect of assimilating them in combination with other types of instrumental measurements was also tested. Correlation and RMSESS metrics were calculated using the reconstruction based on temperature and sea-level pressure data (exp_TP) as a baseline. The exp_TP experiment shows a clear positive impact on the correlation values of all three variables in both seasons. Due to the high skill of temperature and sea-level pressure reconstructions, further improvements with the assimilation of precipitation amounts or the number of wet days are mainly seen in the precipitation field. Assimilating the precipitation amount or the number of wet days has a small impact on the temperature and sea-level pressure correlation skills. If only precipitation information is available from a given region, its affects the reconstruction of the other fields more; for example, over northern and western Australia, where no temperature observations are assimilated, precipitation information improved the skill of the temperature field in boreal winter (Fig. S4). The RMSESS values of all variables from exp_TPR and exp_TPW are very similar in winter, indicating that both types of precipitation data perform similarly. Moreover, assimilating the precipitation amount or wet day records has a mainly positive effect in the regions where temperature observations are absent. However, sea-level pressure fields suggest that with the assimilation of precipitation observations the skill of the pressure reconstruction decreased in the regions with available sea-level pressure measurements. In general, assimilating precipitation amounts (exp_TPR) performs worse than assimilating wet days (exp_TPW) in summer. The better performance of wet days in summer is expected due to lower spatial variability. However, even wet day records have a negative effect on the sea-level pressure RMSESS, especially in Europe. To minimize this negative effect, ignoring or limiting the cross-covariance updates between precipitation and other variables will be tested in future experiments.

5.2 Case study: 1842 drought

A case study was conducted to test how well precipitation amounts and wet day records are able to reproduce a severe drought event in Europe. Only a sparse network with six stations provided data from 1842. In the model simulations (CCC400) the precipitation anomalies are much smaller compared to the reconstructed precipitation fields in the exp_R and exp_W experiments, indicating that with the assimilation of six stations this drought event is reconstructed. The reconstructed precipitation anomaly fields are very similar in the case of both precipitation amounts and the number of wet days (Fig. 7), which is promising, since in the past more wet day records are available. In the exp_TP experiment the dry event is also present, especially in August; however, the number of temperature and pressure observations that were assimilated is much higher than the six available precipitation stations.

Based on the historical sources gathered by Brázdil et al. (2019) a significant precipitation deficit was reported in Bohemia from April to December, causing a forest fire in June. Similarly, a very dry period was recorded in Prague between April and August, with extreme dryness in August. Besides the Czech lands, several documentary data between April and August describe similar conditions in Germany, like fires in Hamburg, a lack of drinking water, a low water level in the Danube, and praying for rain. In the GHCN-Daily dataset there are observation measurements from these regions, and in the exp_R and exp_W experiments notable negative relative precipitation anomalies are apparent from April to August (Fig. 7), with August being the most severe in accordance with the documentary data. In 20CRv3 the relative precipitation anomaly in August does not appear to be as large over this part of Europe; however, no pressure data are included from eastern–central Europe. A larger precipitation deficit is depicted over northern Europe in August in 20CRv3, which can be supported with documentary data about forest fires in Norway (Brázdil et al., 2019). Several other documentary records across Europe from France to Transylvania describe an extraordinary drought, which is reflected with varying intensity in the precipitation anomaly fields of the exp_R and exp_W experiments as well as in the 20CRv3 reanalysis.

6 Conclusions

As the application of data assimilation techniques has become more widespread in the field of paleoclimatology, more and more different observational sources, such as early instrumental measurements, documentary records, and various types of proxy records, have been used in the assimilation process. In this paper, new observation data sources – precipitation amounts and the number of wet days – were tested in an offline ensemble-based Kalman filter framework. The experiments in which only precipitation amounts (exp_R) and only wet day records (exp_W) are assimilated performed similarly in winter, but in summer exp_W has better skill in the case off all three examined variables. Moreover, the results of the exp_W_2L experiment suggest that the localization used for wet days could be increased, by which a better use of observation data can be achieved. In the exp_TPR and exp_TPW experiments, the skill of the two reconstructions was compared to exp_TP to examine the effect of adding precipitation data to the assimilated observations. In general, both precipitation amounts and the number of wet days had a rather positive impact on the temperature reconstructions in winter, while in summer only the number of wet days had an overall positive effect on temperature. The skill metrics of sea-level pressure clearly indicate that precipitation data should be used if pressure measurements are not available from a given region. Therefore, it might be better to limit how precipitation data can update the other fields of the state vector, for example with the implementation of an asymmetric localization function, in forthcoming experiments. The reconstructed monthly precipitation fields of the severe 1842 drought in Europe are mostly in agreement with documentary data, showing that precipitation amounts or wet day records can be useful sources in a paleoclimate data assimilation framework.