Assimilation of Pseudo-Tree-Ring-Width observations into an Atmospheric General Circulation Model

Paleoclimate Data Assimilation (DA) is a promising technique to systematically combine the information from the climate model simulations and the proxy records. Here we investigate the assimilation of Tree-Ring-Width (TRW) chronologies into an atmospheric global climate model using Ensemble Kalman Filter (EnKF) techniques and a process-based tree-growth forward model as observation operator. Our results, within a perfect-model experiment setting, indicate that the “on-line DA” approach did not outperform the “off-line” one, despite its considerable additional implementation complexity. On the other 5 hand, it was observed that the nonlinear response of tree-growth to surface temperature and soil moisture does deteriorate the operation of the time-averaged EnKF methodology. Moreover, for the first time we show that this skill loss appears significantly sensitive to the structure of growth rate function, used to represent the Principle of Limiting Factor (PLF)s within the forward model. Additionally, our experiments showed that the error reduction achieved by assimilating a pseudo-TRW chronologies is modulated by the strength of the yearly internal variability of the model. This result might help the dendrochronology 10 community to optimize their sampling efforts.


Introduction
The low-frequency temporal variability in the climate system cannot be estimated from the available time span of instrumental climate records.Accordingly, paleoclimate reconstruction must necessarily rely on the use of the paleoclimate proxy records.These natural archives exhibit several problematic features, e.g., low time resolution, sparse and irregular spatial distribution, complex nonlinear response to climate, and high noise levels.Therefore the proper extraction of the climate signal contained therein can often remain opaque (Evans et al., 2013).To date, many different ideas have been proposed in order to link proxy records to the paleoclimate conditions where they were created, e.g., data-driven statistical techniques, climate model hindcasts and Bayesian probabilistic methods (see Crucifix, 2012, for a recent review).Among this plethora of approaches, data assimilation (DA) methodologies are today particularly appealing as they deliver estimates of paleoclimate quantities by systematically combining the information of paleoclimate records with the dynamical consistence of climate simulations (Brönnimann, 2011;Hakim et al., 2016).
An important difference between paleo-DA and traditional meteorological DA is that the assimilation period might be very long compared to the timescales of the dynamical model.Under these conditions, the randomizing action of the chaotic model dynamics becomes dominant and consequently the forecast appears completely de-correlated from the previous analysis state.This phenomenon, currently referred to as an "off-line regime", has been observed in several paleo-DA studies (Huntley and Hakim, 2010;Bhend et al., 2012;Pendergrass et al., 2012;Matsikaris et al., 2015).Furthermore, some recent studies have assumed the presence of the off-line condition and accordingly have removed the reinitialization step after assimilation altogether (Steiger et al., 2014;Dee et al., 2016;Hakim et al., 2016).These types of DA methodologies will be referred to in this paper as "offline DA" techniques, in order to contrast them with traditional "online DA" techniques, where the state of the model is updated after the assimilation of observations.Note that despite their lack of accumulation of observational information over time, off-line DA methodologies have already been shown to be more robust than traditional climate field reconstruction (CFR) techniques based on orthogonal empirical functions and stationarity assumptions (Steiger et al., 2014;Hakim et al., 2016).
A typical assumption in most of the paleo-DA studies so far conducted is that the climate-proxy relation is linear.Nonetheless, currently it is widely recognized that climate proxies are the result of complex recording processes, which can be of a physical, chemical and biological nature.More realistic methodologies have been recently sculpted by the paleoclimate community in order to investigate the climate-proxy relationship, considering the distinct processes whereby the climate signal is recorded in proxy archives.Proxy forward modeling (Hughes et al., 2010;Evans et al., 2013) appears to be one of the most promising methodologies in this area.In a proxy forward model the climate forcing is used as input data for producing the artificial proxy records which can be directly compared with the actual ones.One application of proxy forward models is the prediction of the evolution of proxy archives (Vaganov et al., 2006).They can also be applied within climate reconstruction strategies by using probabilistic inversion methods like Bayesian hierarchical modeling (Tolwinski-Ward, 2012), Markov Chain Monte Carlo (MCMC) (Boucher et al., 2014) and DA (Hughes et al., 2010).
Several recent studies have investigated the applicability of process-based forward models in a paleo-DA setting (Hughes et al., 2010;Acevedo et al., 2015;Dee et al., 2016;Hakim et al., 2016).Acevedo et al. (2015) (AC15, hereafter) utilized the process-based tree-ring-width (TRW) forward model Vaganov-Shashkin Lite (VSL) (Tolwinski-Ward et al., 2011) and an online EnKF scheme to assimilate pseudo-TRW records into a chaotic two-scale dynamical system.They found that the nonlinearities of the forward model may deteriorate the performance of the EnKF.Furthermore, they observed that this loss of skill may be ameliorated by means of a fuzzy logic (FL)-based extension of the VSL model.Matsikaris et al. (2015) compared the offline and online implementations of a "degenerate particle filter" applied to a low-resolution Earth system model.They found similar skill for both methods on the continental and hemispheric scales.Nonetheless, they concluded that in the off-line method the temporal consistency of the model is lost.On the other hand, Dee et al. (2016) used three different nonlinear proxy forward models (including VSL) and an off-line EnKF scheme to assimilate TRW, coral and ice core records into two different isotope-enabled atmospheric general circulation model (AGCMs).They demonstrated that the linear-univariate models for tree-ring width may not capture the AGCM's climate, especially for regions where the tree's growth is dominated by moisture.
This paper follows the rationale of AC15 but within a more realistic scenario, where an AGCM is used as a dynamical system and the observational network resembles the currently available TRW chronologies.The purpose of this study is then to contribute to the present knowledge of paleo-DA techniques by addressing the following two questions: i.Does the off-line regime naturally appear for the assimilation of TRW records into an AGCM?
ii.Is the FL-based extension of the VSL model still useful to improve the performance of a time-averaged EnKF technique when a climate model is used?
This study is structured as follows.In Sect. 2 we describe the DA technique, the TRW forward model and the climate model as well as the experimental setting used.Our numerical results are shown in Sect.3, followed by a discussion in Sect. 4.

Data assimilation basics
In this paper, the term DA designates the process of estimating the state of a system using observations and the physical laws governing the evolution of the system as represented in a numerical model (Talagrand, 1997).In a typical DA scheme, a dynamical model is integrated in time until observations become available.Afterwards, the predicted state, also known as forecast, is "updated" using the observational information in order to obtain a corrected state, also known as the analysis.Finally, the model is re-initialized from the analysis state and propagated in time until the next assimilation time, completing the "analysis" cycle.DA methods have evolved from very empirical approaches, such as Newtonian relaxation, to probabilistic ones that attempt to estimate the probability density function (PDF) of the model state conditional on the observations (see Kalnay, 2003, Lahoz et al., 2010and Reich and Cotter, 2015, for reviews).
Among the currently available DA techniques, EnKF (Burgers et al., 1998) occupies an outstanding position for several reasons.It offers an appealing trade-off between accuracy and computational expenses.It works robustly for sparse observational networks and a moderate number of ensemble members (Whitaker et al., 2009).Furthermore, EnKF's implementation does not require any modification of the model's code and uncertainty estimates can be directly obtained from the ensemble spread (Hamill, 2006).The main disadvantage of EnKF, within a paleoclimate setting, is its inability to handle non-Gaussian PDFs, which can easily arise from the nonlinearities of climate models and observation operators.Recently, there have been several developments in the field of nonlinear DA for high-dimension systems (Van Leeuwen et al., 2015); however, at the present fully non-Gaussian DA techniques are still prohibitively expensive to run for general circulation climate models.

Kalman filter
Within the Kalman filter (KF) (Kalman, 1960), the PDF of forecast state p(x) is assumed to be given by a Gaussian function with mean x f and covariance P f : The observations y(t j ) are also assumed to have Gaussian errors and therefore the conditional probability of the observation vector y given the state x is where H is the observation operator and R is the observation error covariance matrix.Following the Bayes theorem, the conditional probability of the state given the observations, i.e., the analysis PDF, is Finally, assuming that H is a linear function, p (x | y) is also a Gaussian function whose mean and covariance can be calculated by the Kalman update equations (Lorenc, 1986): where the Kalman gain matrix K is given by:

Ensemble Kalman filter
For realistic geophysical models, the dimensionality of the model state can be very high and then the calculation and storage of the covariance matrices can be prohibitively expensive.A solution to this problem is provided by the EnKF (Evensen, 1994), which uses an ensemble of model states (X(t) = (x 1 , . .., x m )) to approximate KF equations.Following this approach, the mean and covariance of the forecast take the following form: Here X f ∈ R n×m denotes the forecast ensemble deviation matrix: where e = (1, . .., 1) ∈ R m .An analysis ensemble whose covariance satisfies Eq. ( 5) can be generated in different ways, which can be classified into two main families: stochastic and deterministic filters (Hamill, 2006).
In the stochastic approach an observational ensemble Y is generated by adding realizations of the observational noise to the observation vector y.The analysis ensemble is then created by the following update equation: In the deterministic approach, instead of creating an ensemble of observations, the analysis mean and deviations are calculated using update formulae which do not involve random numbers (see Tippett et al., 2003 for further references).
A practical problem of EnKF schemes is that due to the limited ensemble size, the forecast uncertainty is usually underestimated.This leads to an excessive confidence in the forecast, and after several assimilation cycles the observations may be completely ignored.This situation is normally avoided by means of an ad hoc procedure known as "covariance inflation", where the forecast covariance matrix is multiplied by a constant greater than 1.Another undesired consequence of the limited ensemble size is that the ensemble state at any grid point will present non-negligible spurious correlations with observations located far apart in space.This difficulty is solved using another ad hoc procedure known as "covariance localization".Here, we utilize the R localization (Hunt et al., 2007), where the entries of the observation error covariance matrix are multiplied by a function that increases exponentially with distance.

Time-averaged ensemble Kalman filter
The EnKF algorithm was initially designed to estimate the instantaneous state of a model given instantaneous observations.As a consequence, EnKF cannot be directly applied to paleoclimate data given that the observational information present in proxy records is typically the average of a function of the state over long time periods.A solution to this conflict is provided by the time-averaged EnKF (Dirren and Hakim, 2005), where the instantaneous forecast is decomposed into its time-averaged part and the anomalies around it.Afterwards, the original EnKF update formula is used to assimilate the time-averaged observations into the time-averaged forecast, obtaining the time-averaged analysis.Finally, the instantaneous analysis is form by adding the unaltered timeaveraged forecast anomalies to the time-averaged analysis.This approach is based on the assumption that the observations can only contain time-averaged information (Dirren and Hakim, 2005).If the time-averaged period is considerably longer than the timescales of the dynamical model, which may easily be the case for paleoclimate records, after assimilation the ensemble spread can reach climatological levels, leading to a complete lack of estimation skill for the forecast quantities.This behavior, currently referred to as the off-line regime, was first observed for the time-averaged EnKF applied to a quasi-geostrophic atmospheric jet model (Huntley and Hakim, 2010;Pendergrass et al., 2012).Afterwards, several studies used the simplified off-line time-averaged EnKF approach with global climate models (Bhend et al., 2012;Steiger et al., 2014;Dee et al., 2016) assuming the presence of the off-line regime.However, to our knowledge, there had not been numerical evidence of the onset of off-line conditions for a full time-averaged EnKF algorithm applied to an AGCM.As mentioned in the introduction, filling this knowledge gap is one the objectives of this paper.

TRW forward model
The VSL model for TRW chronologies offers an intermediate-complexity approach between ecophysiological and completely data-driven models (Tolwinski-Ward et al., 2011;Tolwinski-Ward, 2012), where the climatedriven component of tree-ring growth is parameterized by way of a simple representation of the principle of limiting factors (PLF) (Fritts, 1976).The biological concept of PLF states that the pace at which a plant develops is controlled by the single basic growth resource, typically either energy or water, that is in shortest supply.Within VSL the limiting factors considered are near-surface air temperature (T ) and soil moisture (M).These variables influence tree-ring growth via the "growth response" functions and where is the piece-wise linear "standard ramp" function (Tolwinski-Ward et al., 2014 T L and M L denote minimum thresholds for temperature and moisture below which there is no growth, and T U and M U are upper thresholds above which tree growth is optimal.Afterwards, the growth rate G MIN is determined by the smallest growth response, i.e., Finally, the yearly TRW values W are obtained by integrating in time the growth rate modulated by the relative insolation I : Here n stands for the number of the year, t n is the time at the end of the year n and τ is the length of the year (12 months).
Regarding DA, VSL presents two challenging nonlinear aspects: (i) a "thresholded response" (Evans et al., 2013), arising from the insensitivity of trees to climate variability during dormancy and optimal growth, and (ii) "abrupt shifting of recorded variable" (Acevedo et al., 2015), coming from the structure of Eq. ( 12).There, the use of the minimum function implies that tree growth at a particular time can be limited by either temperature or moisture.Accordingly, transitions between growth limitation regimes necessarily happen in a sudden manner.

VSL from the fuzzy logic viewpoint
The term fuzzy logic was coined by Zadeh (1975) and refers to a mathematical theory which has been very successful in modeling complex systems involving imprecise data and vague knowledge of the underlying mechanisms.Since its introduction, FL has greatly influenced many disciplines, most notably control theory (Nguyen et al., 2002).Within the environmental sciences, FL has been applied in ecological and hydrological modeling (Marchini, 2011;Salski, 2006;Se, 2009).Regarding climate proxy forward modeling, AC15 recently showed that the VSL model can be completely embedded into the framework of FL.Within this reinterpretation, the growth response function g T (g M ) corresponds to the membership function to the set S T (S M ) of optimal temperature (moisture) conditions for tree growth.Temperature (moisture) values lying below T L (M L ) present null values for g T (g M ) and accordingly do not belong to S T (S M ).On the other hand, temperature (moisture) values lying above T U (M U ) lead to g T (g M ) values equal to 1, meaning they belong completely to S T (S M ).All the other temperature (moisture) conditions present growth responses between 0 and 1, and consequently they are considered to belong partially to S T (S M ).This idea of partial membership is the basis of fuzzy logic and the sets defined this way are called fuzzy sets.Furthermore, the intersection of the fuzzy sets S T and S M is again a fuzzy set S T∧M , whose membership function can be calculated by evaluating the minimum between G T and G M : Equation ( 14) is completely equivalent to the Eq. ( 13).Then VSL's growth rate function can be interpreted as the membership function for the fuzzy intersection set S T∧M .In FL theory, the minimum function (Eq.14) is one of the most popular representations of the intersection operation; however, it is not the only one as a whole family of appropriate functions actually exist referred to as t-norms (see Nguyen et al., 2002).In AC15 a number of t-norms was tested as replacement for VSL's growth rate function within a highly simplified paleo-DA setting.In particular, it was found that the product t-norm g T∧M = g T • g M might significantly improve the performance of the time-averaged EnKF technique.Accordingly, beside the minimum t-norm, we also consider the product growth response VSL-Prod in this paper: An important aspect of the product growth response VSL-Prod is the presence of an additional growth limitation regime where T and M concurrently limit tree-ring growth.This "co-limitation" regime allows a gradual transition between temperature-and moisture-limited growth limitation regimes and accordingly a progressive alternation of the recorded variable.Growth co-limitation was initially recognized by Singh and Lal (1935), who called attention to the limitations of the original formulation of the PLF due to Blackman (1905).Within the ecological research community, co-limitation has been widely acknowledged as an crucial resource limitation phenomenon (Harpole et al., 2011), with abundant observational support both in terrestrial (Niinemets and Kull, 2005) and aquatic (Saito et al., 2008) environments.Interestingly, within the vegetation modeling community the original PLF formulation is still the predominant approach to model photosynthesis (Yin and Struik, 2009).

Experimental design
Following the rationale used in the experiments of AC15, we conducted a set of DA experiments using the Simplified Parameterizations, primitivE-Equation DYnamics (SPEEDY) model (Molteni, 2003) as a dynamical system and the VSL forward model as an observation operator.The time-averaged state of the atmosphere is estimated via the EnKF approach of Dirren and Hakim (2005).In the following, we describe in detail each of the components of our experimental setting.

Atmospheric general circulation model
The SPEEDY model (Molteni, 2003) is an intermediatecomplexity atmospheric general circulation model (AGCM) comprising a spectral dynamical core and a set of simplified physical parameterizations, based on the same principles as state-of-the-art AGCM but tailored to work with just a few vertical levels.Regarding the interaction with the ocean, SPEEDY offers two possible configurations: (i) prescribed ocean -sea surface temperature is directly imposed as forcing; (ii) slab ocean -the atmospheric model is coupled to a slab ocean model (q flux adjusted mixed layer model) forced by climatological ocean dynamics.Despite its low resolution and the relatively low complexity of its parameterizations, SPEEDY still captures many observed global climate features in a realistic way, while its computational cost is at least 1 order of magnitude lower than the one of sophisticated state-of-the-art AGCMs at the same horizontal resolution (Molteni, 2003).The latter makes SPEEDY especially suitable for studies involving long ensemble runs, like the ones necessary for this study.

DA technique
The SPEEDY model was embedded by Miyoshi (2005) into the SPEEDY-LETKF framework, which offers a parallel implementation of a local ensemble transform Kalman filter (LETKF) (Hunt et al., 2007).Among the different types of EnKF, LETKF is particularly promising for high-resolution models given that the calculation of the analysis for a particular grid point requires only the information of the neighboring observations.Therefore, LETKF offers outstanding scalability properties.SPEEDY-LETKF is an open-source software which has already been used for several DA studies (Li et al., 2009;Miyoshi, 2010;Lien et al., 2013;Ruiz et al., 2013;Amezcua et al., 2014).Here, SPEEDY-LETKF was extended to allow the assimilation of pseudo-TRW observations by means of either online or off-line time-averaged EnKF methods.
For the experiments presented in this paper, we employed ensembles of 24 members due to computational constraints.We used a constant multiplicative inflation of 1 % after the ensemble update and R localization via the following formula: where r h and r v stand for the horizontal and vertical distances, respectively.Their corresponding scaling parameters were set to the values of λ h = 500 km and λ v = 0.4 ln p.

Numerical experiments
We performed a set of Observation System Simulation Experiments (OSSEs) (see Fig. 1), consisting of (i) a single model trajectory x Nature , referred to as "true" run or "nature" run, that is used as a prediction target, (ii) pseudoobservations created by applying the observation operator to x Nature and adding simulated observational noise, (iii) an observationally constrained ensemble run X DA , where the pseudo-observations are assimilated, and (iv) a free ensemble run X Free , where no observations are assimilated and then the ensemble just freely evolves under the action of the model  dynamics.X Free is intended to provide a benchmark of performance, against which it is possible to assess the added value of the DA scheme.
Initially, a 1-year-long spin-up run is performed starting from 1 January 1860.Afterwards, the final state of this model trajectory is used as the initial condition for a 150-year-long nature ("true") run.The ensemble runs with and without DA are identically initialized from a set of states gathered from the last 2 months of the spin-up run (lagged 2-day initialization).Note that the nature run and the different ensemble runs are generated with the same time-varying forcing fields.Regarding the atmosphere-ocean coupling, we used SPEEDY's slab ocean configuration, motivated by the fact that the slow variability in the slab ocean may lend predictability to the atmosphere.In these conditions, the online DA technique should have higher chances to outperform the off-line technique.

Observation generation
Pseudo-TRW observations are produced following VSL's formulation plus a final white noise addition step, in which random draws from a Gaussian distribution are imposed on the time-averaged observations.Noise levels are assessed by means of the signal-to-noise ratio (SNR), given by the ratio of the standard deviation of the unpolluted pseudo-TRW observations to the standard deviation of the additive white noise.Most of the results reported in Sect. 3 correspond to the optimistic value SNR = 10, with the exception of the last sensitivity study, which analyzes the dependency on observational noise levels.Regarding the geographical distribution of the observations, we place a pseudo-TRW chronology at every grid box where at least one actual TRW chronology from the database of Breitenmoser et al. (2014) is present.This strategy yields a realistic observational network comprising 257 stations (see Fig. 2).
Concerning the configuration of the observation operator, we focus our study on the role of the growth rate function by configuring VSL in such a way that no thresholded response takes place.This is done by setting the upper and lower response thresholds to the maximum (minimum) values during the nature (true) run so that the response functions reduce to linear rescaling operators.We consider three different growth rate functions leading to three VSL configurations: (i) VSL-T, where the growth rate is directly given by the growth response to temperature (G T = g T ); (ii) VSL-Min, where the original "minimum" t-norm G MIN is used; and (iii) VSL-Prod, where the FL-based "product" t-norm G PROD is utilized.Note that within the setup described above, VSL-T generates purely temperature-limited pseudo-chronologies, while VSL-Min and VSL-Prod generate mixed ones, where temperature and soil moisture determine tree-ring growth pace.
Finally, respecting the soil moisture fields used to drive the VSL model, we consider two options: (i) extracting soil moisture time series from the climatological surface boundary conditions of the SPEEDY model and (ii) using the precipitation and temperature output of SPEEDY as input for the leaky bucket model (LBM)1 (Huang et al., 1996), as it is normally done with the VSL model.Note that with the first option, the time series obtained present only intra-annual variability (yearly periodicity), while with the second option moisture time series do exhibit interannual variations.

Diagnostics
Thanks to the availability of the truth model evolution for our OSSEs, the forecast and analysis skill of the ensemble runs can be directly assessed.Given the annual resolution of TRW chronologies, we study the filter performance for yearly averaged values of near-surface temperatures.We focus our analysis on near-surface temperature due to the larger error reduction in this field as compared to other variables (e.g., humidity, u wind, v wind) when DA is applied.The behavior of ensemble runs is monitored by means of the root mean square error (RMSE) for the ensemble means.The results are shown as (1) time series of globally averaged temperature RMSE, (2) histograms of these time series and (3) maps of time-averaged (150 years) temperature RMSEs.

Unconstrained ensemble run
An AGCM is an example of a nonautonomous system, and accordingly the evolution of its state is determined by both the atmospheric dynamics and the external forcing.The influences of these two distinct factors can be disentangled to some extent by considering atmospheric variability to be a superposition of an internal component, caused by the intrinsic dynamics, and an external one, resulting from the variations in the boundary conditions (Deza et al., 2014).Under this assumption, internal and external variability can be separated by way of a free ensemble run, using the ensemble mean as an estimate of the forced variability and the ensemble spread as an estimate of the internal variability.Following this train of thought, Fig. 3a shows the magnitude of the yearly internal variability in near-surface temperature.Geographical regions around the Equator present negligible yearly internal variability, which can be explained by the fact that these areas are dominated by the daily cycle, and accordingly the 1-year-long averaging strongly attenuates their internal variability.On the other hand, maximum spread values occur at high latitudes around ±70 • .These yearly internalvariability maxima can be related to planetary-scale patterns fluctuating over timescales longer than 1 year, such as the "annular modes" (Thompson and Wallace, 2000) and displacements of the jet stream (Woollings et al., 2011).Regarding the error of the free ensemble run, Fig. 3b shows that the RMSE maxima clearly coincide with the internal-variability maxima, given that quantities without internal variability can be well described by an ensemble run without DA.From that, it becomes clear that in our setting, DA can only be beneficial to estimate quantities displaying significant yearly internal variability.

Assimilating linear univariate time-averaged observations
In order to focus on the comparison between online and off-line DA techniques, we first consider the assimilation of temperature-limited pseudo-chronologies produced with VSL-T.Note that this observation operator generates linear univariate time-averaged observations, and accordingly the time-averaged EnKF must be in good conditions to operate, given that no nonlinearities are present in the observation operator.
Figure 4b and c show how the assimilation effectively reduces the error of the analysis for both online and off-line methods.All areas adjacent to the observational network show low RMSE values; however, it is important to highlight that stations located in areas of strong yearly internal variability are more efficient than the others at reducing the error of the analysis quantities.An example of this is the chronologies placed in Alaska, where RMSE reductions of up to 1.5 • take place.Conversely, stations situated in regions with weak yearly internal variability do no lead to significant error reductions, as is the case of the chronologies located in the Himalaya area, where SPEEDY presents no significant yearly temperature internal variability to be constrained via DA.Therefore, we claim that the error reduction for analysis quantities is modulated by the magnitude of the yearly internal variability at a specific site.
On the other hand, online forecast quantities do not present significant error reductions, and consequently the online time-averaged EnKF appears to work under the off-line regime.This situation can be seen in Fig. 4a, which presents a strong resemblance with Fig. 3b, corresponding to the RMSE of the free ensemble run.Furthermore, looking at global RMSE time series (Fig. 5), it is evident that the errors for the online technique increase with time, as opposed to the corresponding off-line quantities, which presented stable behavior all along the running time interval.Overall the off-line DA methodology leads to lower error levels on the global scale.The existence of this upward trend for online errors might be attributed to an undesired effect arising from the re-initialization of the ensemble after each assimilation step.For all our online DA experiments we observed the aforementioned lack of forecast skill.Consequently, for the rest of this section we will focus on the off-line DA technique.

Assimilating mixed TRW chronologies
Here we analyze the role of the structure of the growth rate function in the performance of the off-line DA scheme.Given that both G Min and G Prod are nonlinear functions, their use within the observation operator degrades the operation of the time-averaged EnKF, as it is illustrated by Fig. 6.Some areas displaying very low RMSE values for VSL-T, such as Alaska, Siberia and central Asia, present considerably higher error levels when VSL-Min or VSL-Prod are utilized.This loss of skill to nonlinearities in the observation operator is also apparent in the global RMSE time series, as can be seen in Fig. 7, where VSL-T clearly outperforms the other two observation operators.
Regarding the dependency on the representation of the PLF, Fig. 6 shows that the use of VSL-Prod leads to considerably lower RMSE values in central Europe and the east coast of the United States.This edge of VSL-Prod over VSL-Min is also observed on the global scale, as can be seen in Fig. 7.These results concur with the ones of AC15, where in a highly simplified paleo-DA setting the use of VSL-Prod instead of VSL-Min appeared also beneficial to the performance of the time-averaged EnKF technique.Concerning the soil moisture fields, Fig. 8 shows how using the soil moisture calculated by the LBM reduces the error levels of both VSL-Min and VSL-Prod DA runs.However, these improvements in filter performance are more significant for VSL-Min than for VSL-Prod.
Finally, respecting the dependency of the filter performance on the observational noise levels, Fig. 9 shows how the off-line DA scheme is considerably effective for SNR values higher than 1.Increasing further the observational noise (reducing SNR) leads to a fast increase for the aver-

Conclusions
Using the time-averaged EnKF methodology and a processbased proxy forward model (VSL), we assimilated pseudo-TRW chronologies in an AGCM (SPEEDY).Using a set of perfect model experiments we studied two different aspects of the paleo-DA problem: (i) the onset of the off-line regime in the assimilation of observations averaged during long time periods and (ii) the impact of a nonlinear observation operator on the performance of EnKF-based time-averaged DA approaches.
Our online DA experiments in general showed no forecasting skill, and accordingly they appear to operate under the off-line regime.Moreover, they exhibited a detrimental increasing error trend not present for our experiments with off-  In these conditions, the off-line time-averaged EnKF appears to outperform its online counterpart.This result complements the studies of Huntley and Hakim (2010) and Pendergrass et al. (2012), where the off-line regime was observed to emerge within a DA setting comprising the time-averaged EnKF and a quasi-geostrophic atmospheric jet model.From our point of view, the off-line regime can arise either from the dynamical model or from the DA scheme.
Regarding the dynamical model, the main reason for losing forecast skill is that the between consecutive observations considerably exceeds the predictability horizon of the model, and accordingly the ensemble forecast completely forgets the observational information assimilated at the analysis steps.For SPEEDY, due to its purely atmospheric nature, the 1-year averaging period of TRW chronologies appears to be considerably longer than its dynamical timescales.Therefore, the appearance of the off-line regime is to be expected.
On the other for more realistic climate models, there is already evidence of processes with timescales longer than 1 year (Smith et al., 2012).Examples of these sources of slow variability are (i) the "annular modes" (Thompson and Wallace, 2000), which strongly affect climatic conditions in high latitudes, (ii) the latitudinal oscillation of the cell structure, which influences the position of the jet streams, and (iii) the oscillations of the Intertropical Convergence Zone (ITCZ), which greatly impacts humidity in several parts of the globe (Holton and Hakim, 2013), as well as the El Niño-Southern Oscillation (ENSO), which affects the climate of a large portion of the tropical and subtropical areas.Consequently, we consider that for more realistic paleo-DA settings using more comprehensive Earth system models, the off-line regime is less likely to arise.Regarding the DA scheme, a possible reason for the appearance of the off-line regime is the timeaveraged update strategy (Dirren and Hakim, 2005), which assumes complete decorrelation between time-averaged and instantaneous variables.This condition might not be satisfied by SPEEDY, and accordingly the estimation of instantaneous quantities might be badly affected.Another limitation of the q q q q q q q q q q q q q q q q q q Free DA_Min DA_Min_LBM DA_Prod DA_Prod_LBM 0.25 0.30 0.  time-averaged EnKF is its reliance on Gaussianity.This assumption might be easily violated in a climate model, for example by definite positive quantities such as humidity.Consequently non-Gaussianity might be also one of the culprits of the onset of the off-line regime in our experiments.
Concerning the influence of nonlinearities in the observation operator, the performance of the off-line time-averaged EnKF appeared to be significantly sensitive to the selection of the t-norm used to represent the PLF.In our experiments, the product t-norm outperformed the original minimum tnorm, as previously observed for a two-scale Lorenz (1996) model in AC15.Tolwinski-Ward et al. (2014) described trees as fundamentally lossy 2 recorders of climate, due to the integrated nature of the information contained in them and the standardization process used to minimize the non-climatic effects on growth.In the same vein, we argue that the "abrupt shifting" of the recorded variable (temperature or moisture) -implied by the minimum function used in VSL's original formulation -might constitute an additional source of lossiness, which can be reduced by resorting to an FL-based representation of the PLF.In particular, for the product t-norm, the existence of an additional co-limitation regime makes a smoother shifting of the recorded variable possible.As a cautionary remark, we want to highlight that the pseudoobservations assimilated in our experiments present several important limitations: (i) the thresholded response of trees to temperature and moisture was not considered in order to focus on the role of the growth rate function; (ii) VSL's parameters were set in a completely homogeneous fashion for all the observational stations, whereas actual TRW networks are strongly heterogeneous, comprising chronologies generated under highly dissimilar growth limitation regimes.More realistic TRW assimilation experiments will probably have to address these issues as well as the necessity of considering model errors by conducting imperfect model OSSEs.
Finally, we want to mention that the translation of VSL into the FL language suggests other possible extensions for VSL.(i) Growth response functions can be generalized using the extensive knowledge on membership functions gathered in the FL research community (Nguyen et al., 2002).In particular, it might be possible to tailor the shape of the growth response functions so as to optimize the performance of VSL regarding the particular application at hand, e.g., Gaussian DA, as was done in this paper.(ii) Additional limiting factors, e.g., CO 2 concentration, can be easily incorporated into the fuzzy inference system by adding new rules designed to express expert knowledge about the influence of these factors on tree growth.(iii) The intrinsic uncertainty regarding VSLs parameters might be taken into account via the emerging theory of stochastic FL (Luhandjula and Gupta, 1996).Furthermore, this ability of the FL approach to efficiently simulate complex processes involving vaguely understood mechanisms can be used in the development of new proxy forward models as well as in the extension of the existing ones.

Figure 1 .
Figure 1.Schematic of a typical observation system simulation experiment (OSSE) with ensemble online (with cycling) and offline (no-cycling) DA methods.t designates the time axis and X(Y ) denotes the model state (observation) space.Boxes with sharp (rounded) corners represent data (processes).Red (green) vertical bars indicate the forecast (analysis) spread.Vertical dotted lines represent the assimilation steps.

Figure 3 .
Figure 3. Yearly near-surface temperature spread (a) and RMSE (b) for the free ensemble run.

Figure 4 .
Figure 4. Yearly near-surface temperature RMSE for the ensemble run constrained by VSL-T pseudo-TRW observations.(a) Online forecast quantities, (b) online analysis quantities and (c) off-line analysis quantities

Figure 5 .
Figure 5. Global near-surface temperature RMSE for the forecast ensemble run constrained by VSL-T pseudo-TRW observation (red) and the free ensemble run (black) for the analysis of online (green) and off-line (blue) DA.Panel (a) shows time series and panel (b) their corresponding histograms.Horizontal lines represent the mean values.

Figure 6 .
Figure 6.Yearly near-surface temperature RMSE for the analysis of the ensemble runs using off-line DA and climatological soil moisture.(a) Assimilating VSL-T pseudo-TRW observations, (b) assimilating VSL-Min pseudo-TRW observations and (c) assimilating VSL-Prod pseudo-TRW observations

Figure 7 .
Figure 7. Global near-surface temperature RMSE for the off-line analysis of the ensemble runs constrained by VSL-Min (red), VSL-Prod (green) and VSL-T (blue) pseudo-TRW observations.Climatological soil moisture is used to drive the VSL model.Horizontal lines represent the mean values.Panel (b) shows the histograms of the time series.

Figure 8 .
Figure 8. Global yearly near-surface temperature RMSE box plots for the free ensemble run forecast (free) and the analysis of the offline DA runs using VSL-Min observation operator and climatological soil moisture (DA_MIN), VSL-Min observation operator and leaky bucket model (DA_MIN_LBM), VSL-Prod observation operator and climatological soil moisture (DA_PROD), and VSL-Prod observation operator and leaky bucket model (DA_PROD_LBM).

Figure 9 .
Figure 9. Averaged global yearly near-surface temperature RMSE for the analysis of the ensemble run with off-line DA, VSL-Min observation operator and different signal-to-noise ratios.The green star shows the corresponding value for the free run.