3.1 Functional form of the linear response function

In general, the identification of linear response functions from data may be performed either pointwise (e.g., Ragone et al., 2016) or assuming a functional form (e.g., Maier-Reimer and Hasselmann, 1987). Both approaches usually lead to an ill-posed problem and therefore to similar difficulties in finding the solution (see more details in Sect. 3.3.1). Although the RFI method may be applied in either case, here we assume that the response function consists of an overlay of exponential modes. By this ansatz we guarantee from the outset that the response relaxes to zero for t → ∞, which is consistent with the expectation that real systems have finite memory. Besides constraining the function space for the derivation of the response function, another added value of this approach is that in principle it also gives the spectrum of internal timescales of the response.

Assuming this ansatz, the question on the functional form of χ(t) arises. In climate science, it is typically assumed that the response function can be described by only a few exponents (Maier-Reimer and Hasselmann, 1987; Enting and Mansbridge, 1987; Hasselmann et al., 1993, 1997; Grieser and Schönwiese, 2001; Li and Jarvis, 2009; Joos et al., 2013; Colbourn et al., 2015; Lord et al., 2016), i.e.,

where the τ_i values are interpreted as characteristic timescales and the g_i values are their respective weights. τ_i and g_i are then obtained by applying some fitting technique taking a fixed number of terms M. Thus, an important step in this type of approach is to determine a suitable value for M. A common practice is to initially take only a small number of terms M, solve the problem, and then add terms progressively until the addition of a new term does not improve the fit anymore according to some quality-of-fit criterion (e.g., Kumaresan et al., 1984; Maier-Reimer and Hasselmann, 1987; Hasselmann et al., 1993; Pongratz et al., 2011; Colbourn et al., 2015; Lord et al., 2016). Thereby it is assumed that once results stabilize, the information in the data has already been fully exploited, so that fitting of additional terms would be artificial. Nevertheless, finding the parameters τ_i and g_i either from a given χ(t) by Eq. (7) or from ΔY(t) by inserting Eq. (7) into Eq. (1) means solving a special case of a Fredholm equation of the first kind (see Appendix A), which is an ill-posed problem (Groetsch, 1984). This implies that even though the obtained solution may give a very good fit to the data, it may significantly differ from the exact solution (see, e.g., the famous example from Lanczos, 1956, p. 272).

Therefore, to avoid the complication of determining M, we assume instead that the response function is characterized by a continuous spectrum g(τ) (Forney and Rothman, 2012):

Accordingly, we assume that the response is dominated by relaxing exponentials, meaning that potential contributions from oscillatory modes are not distinguishable from noise. By this approach the timescale τ is not an unknown anymore but is given after discretization by a prescribed distribution with M terms covering a wide range of τ_i values. Thus, only a discrete approximation to the spectrum g(τ) needs to be found. In this way the functional representation is made independently of the question of information content as long as the spectrum of discrete timescales is chosen to be sufficiently large and dense to widely cover the spectrum of internal timescales of the considered system.

This approach has an additional advantage. By prescribing the distribution of timescales, one must not solve a nonlinear ill-posed problem (by solving Eq. 7 for τ_i and g_i) but only a linear ill-posed problem (by solving Eq. 8 only for g(τ)), for which the mathematical theory is fairly well developed (Groetsch, 1984; Engl et al., 1996). Because the problem is linear, the solution is even given analytically (see Sect. 3.3.1), which makes the method very transparent.

3.2 Discretized problem

In view of applications to geophysical systems like the climate or the carbon cycle (Part 2 of this study) that are known to cover a wide range of timescales (Ghil and Lucarini, 2020; Ciais et al., 2013, Box 6.1), it is useful to switch to a logarithmic scale (Forney and Rothman, 2012) by rewriting Eq. (8) in terms of log ₁₀τ:

Hereafter, q(τ) and its discrete version q (see below) will be called the spectrum.

In order to apply the basic Eq. (6) together with Eq. (9) to experiment data, the whole problem needs to be discretized in time and also with respect to the spectrum of timescales. Here we assume the data to be given at equally spaced time steps $t_k=t_{\mathrm{0}}+k\mathrm{\Delta }t$, $k=\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{1}$, where N is the number of data, while the timescales are assumed to be equally spaced at a logarithmic scale between maximum and minimum values τ_max and τ_min, i.e.,

where M is the number of timescales. As shown in Appendix B, the resulting discretized equations corresponding to Eqs. (6) and (9) are

and

where η_k stands for the noise. Combining the response data ΔY_k, the spectral values q_j, and the noise values η_j into vectors ΔY∈ℝ^N, q∈ℝ^M, and η∈ℝ^N, these equations can be written in vector form as

with the components of matrix A given by

Matrix A is known from the prescribed spectrum of timescales τ_i and the forcings f_i. Considering η as a fitting error, in principle one can apply standard linear methods to solve Eq. (13) for the desired spectrum by minimizing

where $\left|\right|\cdot \left|\right|$ denotes the Euclidean norm, i.e., $\left|\right|\mathit{x}\left|\right|=\sqrt{{\sum }_i{x}_i^{\mathrm{2}}}$. Here we denoted the spectrum as q_η instead of q to emphasize that the spectrum found in this way can only be an approximation to the original q depending on the noise present in the data.

Unfortunately, it turns out that solving Eq. (15) is not a trivial task. The first difficulty is that the finite information provided by the data makes the problem underdetermined: ideally one wants to obtain a spectrum q(τ) defined for $\mathit{\tau }\in \left[\mathrm{0},+\mathrm{\infty }\right[$, but the data ΔY are discrete and cover only a limited time span. However, the most serious issue in identifying χ(t) arises because Eq. (1) is a special case of a Fredholm equation of the first kind (Groetsch, 1984, 2007, see also Appendix A), where the quest for the integral kernel is well known to be an ill-posed problem (see, e.g., Bertero, 1989; Hansen, 1992). This basically means that any solution q_η of Eq. (15) obtained via classical numerical methods such as lower-upper (LU) or Cholesky decomposition will be extremely sensitive to even small errors in the data (Hansen, 1992). Therefore, to solve Eq. (15) for the spectrum q_η, we invoke regularization.

3.3 Regularization

To treat the ill-posedness of Eq. (15), our RFI method combines techniques from regularization theory with a novel approach to estimate the noise level in the data. To facilitate the understanding of the method, in this section we briefly review these techniques along with some other aspects of the theory that are relevant for our method development.

3.4 Estimating the noise level δ

To introduce our approach, in the following we assume that data from an unforced experiment (control experiment) are available – as is typically the case in applications to Earth system models (see Part 2) – that allow for an independent estimate of the noise level δ.

A naive way to invoke these data to determine λ would be to take δ essentially as the standard deviation of the control experiment – more precisely: $\mathit{\delta }:=\mathit{\sigma }\sqrt{N}=\left|\right|\mathbf{\Delta }{\mathit{Y}}_{\text{ctrl}}\left|\right|$. Technically, to find λ, one way is to start with a large value for λ and decrease it until the left-hand side of Eq. (22) matches δ (as suggested by Hämarik et al., 2011). That this procedure works is explained by the fact that the function $\mathit{\lambda }\mapsto \left|\right|\mathbf{A}{\mathit{q}}_{\mathit{\lambda }}-\mathbf{\Delta }\mathit{Y}\left|\right|$ is continuous, is increasing, and contains δ in its range (Groetsch, 1984, Theorem 3.3.1). Having found λ in this way, the desired solution q_λ is then obtained from Eq. (20). However, this approach is not as straightforward as one may think: because of the forcing, the noise in the perturbed experiment may have different characteristics from that in the control experiment. Therefore in the following we devise a method for how to account for this problem.

Formally in Eq. (13) ΔY is split into a “clean” part and noise η. Entering this into Eq. (19) gives

Accordingly, the first term in the sum gives the “true” solution q, while the second term gives the noise contribution to q_η. As already pointed out when discussing regularization, the “true” solution of ill-posed problems can only be recovered if it is dominated by the projection onto the first singular vectors. This requirement is formally stated by the discrete Picard condition (Hansen, 1990), which demands that the size of the projection coefficients $|{\mathit{u}}_i•\mathbf{A}\mathit{q}|$ drops sufficiently quickly to zero, so that they become smaller than σ_i before σ_i levels off to a finite value because of numerical errors. To find a good estimate for the noise level δ, we use this in the following way. Let i_max be the value of the index i, where the singular values σ_i start to level off. Assuming that the Picard condition holds, one can infer that

Therefore,

This equation determines the high-frequency components of the noise η. It remains to determine also the low-frequency components to obtain an estimate for δ.

For this purpose, we take advantage of the data from the control experiment. The control experiment is an experiment performed for the same conditions as the perturbed experiment, with the only difference that the forcing f is zero, so that the resulting ΔY_ctrl can be understood as pure noise; therefore we write η_ctrl:=ΔY_ctrl. While in the forced experiment the low-frequency noise is obscured by the low-frequency response induced by the forcing, the low-frequency part of the control experiment data can to first order be expected to give an estimate of the low-frequency noise present in the forced experiment. Nevertheless, it is clear that due to the forcing the spectral characteristics of noise may be different in the forced and unforced experiments. More precisely, the spectrum of noise may differ in its overall level and spectral distribution (i.e., the “shape” of the spectrum). In the following, we account for a possible difference in the overall level. However, we will assume that the spectral distribution is approximately the same for η_ctrl and η; we call this the spectral similarity assumption.

After these considerations, λ can be determined as follows: take i_max as the last index i before the plateau σ_i≈0. This i_max distinguishes high-frequency ($i>i_{max}$) from low-frequency ($i\le i_{max}$) components. Then

are the levels of high-frequency noise in the perturbed (see Eq. 25) and control experiments, respectively. We now scale the spectral components of η_ctrl so that its high-frequency level matches the high-frequency level of ΔY:

In this way, the magnitude of the high-frequency components of η^′ matches that of ΔY and because of Eq. (25) also that of η. On the other hand, the spectral distribution of η^′ is the same as for η_ctrl and by the spectral similarity assumption approximately the same as for η. Because η^′ and η have similar spectral distributions, the fact that the magnitude of the high-frequency components of η^′ matches that of η implies that the magnitude of their low-frequency components also matches. Therefore, η^′ can be seen as an estimate of the noise η in the perturbed system not only at high, but also at low frequencies. Hence this corrected noise vector η^′ can be used to obtain an estimate of the noise level of the perturbed system by setting

Compared to taking for δ simply the noise level from the unperturbed experiment (as was insinuated above), taking it in this scaled way ensures that the high-frequency components are consistent with the Picard condition that must hold for q to be recoverable from the ill-posed problem tackled here. Having determined δ, λ can now be computed from Eq. (22) as described above, from which the q follows (Eq. 20) and hence χ(t) (Eq. 12).

3.5 Additional noise-level adjustment in the presence of a monotonicity constraint

In the application to the land carbon cycle in Part 2 of this study, we show that certain response functions χ(t) decrease monotonically to zero. In attempts to recover such response functions by employing the noise-level adjustment described in the previous section, it may turn out that the numerically found response function fails to be monotonic. There may be several reasons for this failure (strong nonlinearities, signal too obscured by noise, etc.). However, one additional reason may be that the low-frequency level of the noise was not properly estimated by assuming that the spectral distribution in the unperturbed experiment reflects the distribution in the perturbed experiment. For such cases one may try to improve the result by further adjustment of the low-frequency noise level to obtain a more reasonable result.

Figure 1Final RFI algorithm (see text for details).

Download

The idea is to adjust the low-frequency components of noise independently of the high-frequency components iteratively until the solution obeys the monotonicity constraint. To understand how to do so, several things must be explained.

1. A sufficient condition for χ(t) being monotonic is that all components q_i have the same sign (see Appendix C). Therefore, starting out from a numerical solution for χ(t), it would develop towards monotonicity if one could come up with a sequence of vectors q_λ with fewer and fewer sign changes.

2. From Eq. (20) it is seen that because of Hansen's observation explained in Sect. 3.3.1, that singular vectors v_i are less noisy for lower i, q_λ has fewer sign changes the fewer v_i contribute to the sum.

3. As seen from Eqs. (20) and (21), this is the case the more components the filter function is suppressing, i.e., the larger the value of λ.

4. To obtain larger values of λ, one sees from the discrepancy method (Eq. 22) that one has to increase δ. The proof for this can be found in Groetsch (1984) (Theorem 3.3.1), but it can also be made plausible as follows: starting from λ=0, q_λ=q_η, which is the solution of the minimization problem (Eq. 15). Hence, for λ=0 the discrepancy on the left-hand side of Eq. (22) is minimal. By increasing λ, one decreases all components of q_λ (Eq. 20), thereby increasing the discrepancy.

5. Following the reasoning of the previous section, in order to obtain a larger value for δ, one must increase the noise level $\left|\right|{\mathit{\eta }}^{\prime }\left|\right|$ (compare Eq. 29). In doing so, one must keep the high-frequency components of $\left|\right|{\mathit{\eta }}^{\prime }\left|\right|$ unchanged because they must keep matching the level of the high-frequency components of the noise in the perturbed experiment η (given by Eq. 25). Hence, to increase δ, one sees from Eq. (29) that this is achieved by scaling up the low-frequency components of $\left|\right|{\mathit{\eta }}^{\prime }\left|\right|$.

Summarizing these considerations, we have to increase the level of low-frequency contributions to η^′ to develop a given solution for χ(t) towards monotonicity.

This leads to the overall algorithm listed in Fig. 1. The first five steps have already been explained at the end of Sect. 3.3.2. To account for monotonicity, the additional step 6 combined with the loop back to step 4 has to be iteratively executed. To enhance the low-frequency noise level as explained above, we calculate in step 6 a new noise vector η_new by keeping the high-frequency part from η and enhancing its low-frequency components by a factor c>1. Then we recompute χ(t) from steps 4 and 5 and once more check for monotonicity.

For the RFI algorithm to be applicable, two conditions must be met: (1) a linear response exists for sufficiently weak perturbation and, (2) in addition to the response experiment, a control experiment is also available. The assumptions needed for the successful application of the algorithm are summarized in Table 1.

Table 1Summary of assumptions underlying the RFI algorithm.

Download Print Version | Download XLSX

4.1 Toy model and artificial experiments

As a toy model we take

Here the matrix $M\in {\mathbb{R}}^{D\times D}$ describes the relaxation of the unperturbed model. The second right-hand-side term represents the deterministic forcing constructed from the time-dependent forcing strength $f:\mathbb{R}\to \mathbb{R}$ and the coupling vector a∈ℝ^D. Additionally, the system is perturbed by the stochastic forcing $\mathit{n}:\mathbb{R}\to {\mathbb{R}}^D$, which for simplicity is assumed to be white noise. To make the relation to the carbon cycle considered in Part 2, the components of x may be understood as the carbon stored in plant tissues and soils at the different locations worldwide, so that the observable $Y\left(t\right):={\sum }_i{x}_i\left(t\right)$ is the analog of globally stored land carbon. The solution of the system is

We assume from the outset M to be diagonal with eigenvalues $-\mathrm{1}/{\mathit{\tau }}_i^{*}$, the ${\mathit{\tau }}_i^{*}$ being the relaxation timescales. Then

with the linear response function χ^*(t) and the noise term η^*(t) given by

To complete the description of the toy model, one has to specify its parameters. For the dimension D we take 70, and the timescales are assumed to be distributed logarithmically between ${\mathit{\tau }}_{min}^{*}=\mathrm{0.01}$ and ${\mathit{\tau }}_{max}^{*}=\mathrm{1000}$, i.e., ${\mathit{\tau }}_i^{*}=\mathrm{0.01}\times {\mathrm{10}}^{i\mathrm{\Delta }\log \mathit{\tau }}$ with $\mathrm{\Delta }\log \mathit{\tau }=({\log }_{\mathrm{10}}{\mathrm{10}}^{\mathrm{3}}-{\log }_{\mathrm{10}}{\mathrm{10}}^{-\mathrm{2}})/\mathrm{70}$. With carbon cycle applications in mind, the distribution of the components of the coupling vector is adapted from the log-normal rate distribution found by Forney and Rothman (2012) for the decomposition of soils:

with μ and σ chosen so that the peak timescale is around τ=5 and the limits of the log-normal distribution are approximately within τ=0.1 and τ=200 (see the “true” spectrum in Fig. 3). The components of n are taken as uncorrelated, i.e., 〈n_i(0)n_j(t)〉=ξδ_ijδ(t), with standard deviation ξ being chosen differently in different experiments.

Table 2Experiments considered in this study. Forcings are shown in Fig. 2. To standardize the types of experiments considered here and in Part 2, we select forcing functions that mimic those employed in climate change simulation experiments to whose data the RFI method is applied in Part 2. Note that in principle any type of forcing could be employed.

Download Print Version | Download XLSX

Figure 2Forcings for the experiments considered in this study. To standardize the type of experiments considered here and in Part 2, we select forcing functions that mimic those employed in climate change simulation experiments to whose data the RFI method is applied in Part 2. Note that in principle any type of forcing could be employed.

Download

In our experiments we explore how Y(t) behaves as a function of the forcing f(t). To this end, we choose a forcing function f(t) (see Table 2 and Fig. 2). The most obvious way to perform the toy model experiments would be to integrate Eq. (30). However, to have better control over the noise, it is for our purpose more appropriate to use the analytical solution (Eqs. 32–34). Hence, we numerically integrate Eq. (32), using the representation Eqs. (33) and (34). The data from these experiments are then used to investigate the performance of the RFI algorithm to recover χ^*(t). Since all a_i values are non-negative, the response function (Eq. 33) is monotonic, so that we apply the extended version of the algorithm (see Fig. 1, including step 6). In all experiments we generate N=140 data points to have a time series of similar length to the climate change simulations analyzed in Part 2 (140 years, one value for each year). To apply the RFI method, the noise from an associated control experiment is also needed. This is obtained from Eq. (34) by using another realization n_i(t) of white noise for each system dimension i.

4.2 Choice of parameters for the RFI method

To apply the RFI method, we choose M=30 timescales for the recovery of χ^*. Using ${\mathit{\tau }}_{min}=\mathrm{0.1}$ and ${\mathit{\tau }}_{max}={\mathrm{10}}^{\mathrm{5}}$, we distribute the spectrum of timescales according to Eq. (10). These parameters are also used for the application on the carbon cycle in Part 2 and for the comparison with previous methods in Sect. 5.

4.3 Ideal conditions

To gain trust in the numerics of our implementation of the RFI method, we present in this section a technical test considering conditions under which it is known that the linear response function should be quite perfectly recoverable. Such ideal conditions are characterized by perfect linearity and absence of noise. Hence we use the presented toy model (which is anyway linear) in the absence of noise (n=0) for this test. Actually, this will not be a full test of the algorithm but only of the implementation of its basic apparatus (Sects. 3.2 and 3.3.1), culminating in Eq. (20) since in the absence of noise, the method to determine the regularization parameter λ (Sects. 3.3.2 and 3.5) is not applicable. One might think that in the absence of noise one could use Eq. (19) to determine the linear response function, but even under such ideal conditions the ill-posedness of the problem calls for regularization to suppress the numerical noise that prevents one from obtaining a sensible solution from Eq. (19) (see the discussion in the paragraph after Eq. 19). However, choosing the small value of $\mathit{\lambda }={\mathrm{10}}^{-\mathrm{8}}$ for the regularization parameter when evaluating Eq. (20) is sufficient for this technical test.

Figure 3c shows the response of the noiseless toy model to the forcings shown in Fig. 2; i.e., we performed the experiments listed in Table 2, although for the present test the control experiment is not needed.

Applying Eq. (20) to the experiment data gives the spectrum q_λ shown in Fig. 3a. Here, we derived the spectrum q_λ for each experiment separately, although in the figure only single dots are seen, because all results coincide so closely and are almost indistinguishable from the “true” solution q^*, as was expected for this ideal case. The next Fig. 3b shows the response function obtained from the spectra q_λ using Eq. (12). Obviously from Fig. 3a the “true” response function is reconstructed perfectly from whatever experiment is used. As a final test we predict using in Eq. (1) the response function obtained from the 1 % experiment the responses of other experiments, and indeed, these predicted responses are indistinguishable from the responses obtained directly from the experiments (see Fig. 3c). This latter result demonstrates perfect robustness of the numerical approach to recover the responses in this ideal case.

Figure 3Demonstration of robust recovery for noise-free data from the toy model: (a) recovered q_λ; (b) recovered χ(t); and (c) original responses and predictions using χ(t) derived from the 1 % experiment. Reconstructed values are almost indistinguishable from original data. To plot the “true” spectrum of the toy model in subfigure (a), we used the relation ${\mathit{q}}^{*}=\mathit{a}/\mathrm{\Delta }{\log }_{\mathrm{10}}\mathit{\tau }$, which can be obtained by comparing Eq. (33) with Eq. (12). Since from the discrete spectrum the response function and the response may be obtained for any time t, the spectrum is plotted as dots, while the response function and response are plotted as continuous lines. The regularization parameter is chosen as $\mathit{\lambda }={\mathrm{10}}^{-\mathrm{8}}$.

Download

4.4 First complication: noise

The presence of noise may severely hinder the detailed recovery of χ_*(t) due to the ill-posed nature of the problem. In order to demonstrate the effect of the addition of noise on the quality of the derived χ(t), we define a relative error for the prediction of the responses from different experiments. Consider a particular experiment – which is in our case the 1 % experiment – from which we have obtained by the RFI method the response function, which we call here χ⁰(t). The relative error for the prediction of the response from an experiment “k” by the recovered χ⁰(t) via the convolution (Eq. 1) is

where ⋆ stands for the discrete form of the convolution operation (Eq. 1) used to predict the responses, i.e., $\mathrm{\Delta }t\sum _i{\mathit{\chi }}_{j-i}^{\mathrm{0}}{f}_i^k$. In the following we denote as ${\mathit{\epsilon }}_k^{\mathrm{0}}$ the prediction error for the experiment “k”. To measure the quality of the prediction across multiple experiments, we also define the mean prediction error

where K is the number of predicted responses. The reader may wonder why we quantify the quality of the recovery only indirectly from the responses found in different experiments and not directly from the recovery of χ(t). The reason is that in real applications the “true” χ(t) is not known but the responses are. The reliability of this indirect measure for the quality of the recovery is discussed in Sect. 5.

To study how the quality of the recovery depends on the noise level, we introduce the signal-to-noise ratio (SNR) of the response data from a perturbation experiment as

where δ is the final noise-level estimate obtained by the RFI method, as described in Sect. 3.3.2 (see Eq. 29).

Figure 4Mean prediction error (Eq. 37) of the recovery when deriving χ(t) for different values of the SNR. As the SNR increases, the recovery of χ(t) improves. To illustrate the most general case where χ(t) is not known to be monotonic, we do not apply the monotonicity check (step 6 of Fig. 1).

Download

To demonstrate the dependence of the mean prediction error (Eq. 37) on the SNR, we performed 1 % experiments using different noise levels. The resulting dependence is shown in Fig. 4. As expected, for a small error a sufficiently large SNR is needed; i.e., a good recovery may be hindered by a too high noise level.

Figure 5Demonstration of the operation of the RFI algorithm in the presence of noise using toy model data from a 1 % and control experiment. To demonstrate the relevance of the noise-level adjustment (step 3 from Fig. 1), the standard deviation of the noise in the control experiment was taken to be 10 times smaller than that for the noise in the perturbed experiment. (a) Picard plot showing the singular values σ_i and the projection coefficients of the data $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$, the “true” noise $|{\mathit{u}}_i•\mathit{\eta }|$, and the final noise estimate $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$; (b) coefficients of the regularized solution (Eq. 20); (c) “true” and recovered linear response functions. Since the RFI algorithm correctly adjusted the noise level to the “true” noise in the data, the resulting regularized solution has contributions only from the first few projection coefficients which are not completely obscured by noise. Overall, the recovery is almost perfect, because the SNR (chosen as about 520) is still sufficiently good and because the noise was chosen to conform with the spectral similarity assumption. The regularization parameter determined by the algorithm is λ≈ 30 364. Because the noise-level adjustment (step 3 from Fig. 1) already gave a good estimate of the “true” noise in the data, no monotonicity check was needed (step 6 from Fig. 1).

Download

In Fig. 5 we demonstrate how the overall noise-level adjustment in step 3 of the RFI algorithm (see Fig. 1) affects regularization to recover the correct response function. To guarantee that the overall level of the noise spectrum is indeed substantially different in the control and perturbed experiments (so that the adjustment is really needed), we take for the noise in the control experiment a standard deviation 10 times smaller than that for the noise in the perturbed experiment. To demonstrate how the adjustment works, it is helpful to consider the so-called “Picard plot”. This type of plot was originally introduced to analyze the spectral characteristics of an ill-posed problem (see, e.g., Hansen, 1992). In Fig. 5a we show the Picard plot for data obtained from a 1 % experiment with the toy model using a SNR of ≈ 520 to ensure a good recovery. The singular values σ_i decrease to extremely small values as the index i increases. This demonstrates that the problem of solving for the response function is indeed ill-posed and therefore regularization is needed for its solution (compare Eq. 19 with Eq. 20). The data labeled by $|{\mathit{u}}_i•\mathit{\eta }|$ are the “true” noise coefficients, obtained by subtracting the “clean” response Aq, known analytically from the toy model description, from the noisy toy model response ΔY. Comparing them to the projection coefficients of the response $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$, one sees that with the exception of the first few coefficients the response is dominated by its noise content. Accordingly, only the information contained in these first few coefficients is recoverable from this ill-posed problem whatever method is used. The data labeled by $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$ have been added to the Picard plot to demonstrate how the RFI algorithm operates: these data are the projection coefficients of the estimated noise content in the data, where η_est is the final value of η^′ obtained by the RFI method. Obviously, the RFI algorithm correctly estimates the “true” noise level not only at high frequencies – where it is correct by the noise-level adjustment in step 3 of the RFI algorithm (see Fig. 1) – but also at low frequencies, where it is predicted from the adjusted low-frequency components of the control experiment (also step 3). Accordingly, in this case the spectral similarity assumption holds, and there is no need to further adjust the noise level (step 6).

How the estimation of the noise in the data and the resulting regularization affects the projection coefficients of the spectrum q can be seen in Fig. 5b: only those few coefficients not dominated by noise contribute to the regularized solution. In this case these few coefficients selected by determining the regularization parameter λ from the noise level are sufficient for an almost perfect recovery of the response function, as seen in Fig. 5c.

It is important to note that in the situation of Fig. 5, where the overall noise level differs considerably in the control and perturbed experiments, a naive noise estimate taken from the control experiment without the adjustment in step 3 (as first suggested in Sect. 3.3.2) would severely underestimate the noise actually in the data. This would in turn lead to an underestimation of the regularization parameter (see Groetsch, 1984, Theorem 3.3.1). As a result, the wrong filtering by regularization would leave projection coefficients dominated by noise in the solution, likely leading to large errors in the recovered response function. This example therefore demonstrates the relevance of the noise adjustment in step 3.

Figure 6Demonstration of the additional noise-level adjustment in the presence of a monotonicity constraint using toy model data from a 1 % and control experiment: (a) Picard plot; (b) coefficients of the regularized solution (Eq. 20) and (c) recovered linear response function. All the figures are based on the same toy model experiments using a SNR = 1189. To demonstrate the effect of the noise-level adjustment, the spectral similarity assumption is broken by artificially increasing the low-frequency components of the noise in the experiments. The plots in the first row show the results from the RFI algorithm in the absence of additional noise-level adjustment (step 6 in Fig. 1). Although the “true” response function of the toy model is monotonic, the response function recovered by the RFI algorithm is non-monotonic (last figure in the first row). However, if the noise adjustment is switched on (second row), the response function is correctly recovered as monotonic (last figure in the second row). Arrows in subfigures (b) indicate the index i_critical that separates components of the solution that are only weakly suppressed (i<i_critical) from those that are almost completely suppressed (i≥i_critical). The regularization parameter determined by the algorithm is λ≈1 for the first row and λ≈ 11 450 for the second. For more details, see the text.

Download

Finally in this section, we demonstrate that by accounting for monotonicity of the linear response function, one may obtain a better estimate of the low-frequency components of the noise whereby the recovery of the response function is improved. In Fig. 6 we plot results from toy model experiments where the spectral similarity assumption does not hold. This was achieved by artificially enhancing the low-frequency components of the noise η^*(t) in Eq. (32). The top row plots show the results from the recovery when the additional noise-level adjustment was not used. Because the spectral similarity assumption does not hold, the estimated low-frequency components of the noise $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$ do not match those of the “true” noise $|{\mathit{u}}_i•\mathit{\eta }|$ (Fig. 6a1). Ideally, only those four projection coefficients of the data $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$ which are larger than the projection coefficients of the “true” noise $|{\mathit{u}}_i•\mathit{\eta }|$ should contribute to the recovered response function. Instead, as seen in Fig. 6b1, the coefficients with indexes between i=4 and i=7 give the dominant contributions because they are larger than the estimated noise coefficients $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$ (compare Fig. 6a1). Therefore, the recovery of the response function is poor (Fig. 6c1). However, since in this case the low-frequency components of noise are such that the recovered response function is non-monotonic although the “true” response function is known to be monotonic, one may further adjust the noise level to improve the results.

This further adjustment is the purpose of step 6 of the RFI algorithm (see Fig. 1). Its effect is demonstrated by the second-row plots of Fig. 6: the estimated noise components now match the “true” noise components better that had been underestimated in the first row (compare Fig. 6a2 and a1), so that only those four components that carry information (compare in Fig. 6a2 the projections $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$ for low index i with $|{\mathit{u}}_i•\mathit{\eta }|$) survive the regularization (Fig. 6b2). As a result, the quality of the recovery of the response function has considerably improved (Fig. 6c2).

4.5 Second complication: nonlinearity

The second difficulty in recovering the linear response function χ(t) from a perturbation experiment may arise from nonlinearities present in the considered system. Generally it must be suspected that nonlinearities are present, so that they should not hurt as long as they are small, and indeed, from the viewpoint of regularization, contributions from nonlinearities can be considered an additional noise, so that in principle they can also be filtered out. However, as with noise, when getting stronger they cause a deterioration of the recovery of the response function. In the following, we show this more formally and discuss in detail how the RFI algorithm behaves in the presence of nonlinearities.

To understand how contributions from nonlinearities affect the recovery of the response function, we write the nonlinear terms in Eq. (5) collectively as $\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(t\right)$. This formally gives

instead of Eq. (13). Plugging this into Eq. (20), the spectrum is obtained as

Accordingly, the nonlinear contributions can be understood as an additional noise in the spectrum q_λ, so that the theory of regularization fully applies when replacing η by the combined noise $\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}$. Hence, as in their absence, nonlinearities do not prevent the application of regularization as long as the signal is not buried under this combined noise.

However, for the RFI algorithm to give good results, a second condition is that the contributions from $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ must not be large compared to those from η. To understand this, one must realize that the response and with it the nonlinear contributions $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ are dominated by low-frequency components because of the low-frequency nature of the forcing for the problems of interest (for instance in % experiments). The RFI algorithm uses an estimate for the noise level in the perturbation experiment obtained from the control experiment assuming that the spectral distribution is approximately the same in the noise from the control experiment and the noise in the data from the perturbation experiment (spectral similarity assumption; step 3 of Fig. 1). However, the control experiment does not contain any contributions from nonlinearities because the forcing is zero. Therefore, if in the data from the perturbation experiment the contributions from nonlinearities $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ are not small compared to those from η, the spectral similarity assumption does not hold. Since this assumption is at the heart of the RFI algorithm, its breakdown leads to a poor recovery of the linear response function.

All this is demonstrated in the following by toy model experiments. For this purpose, we artificially consider the response of the toy model not in Y but in its nonlinear transform

where the parameter a determines the strength of the nonlinearity. The particular functional form chosen for Y_nonlin(t) mimics the nonlinear effect of saturation encountered for instance in the land carbon sink when atmospheric CO₂ rises to high values. In the following, to demonstrate the effect of nonlinearities, we set the noise level in the toy model experiments to a rather small value in order to have a good SNR in the experiments considered.

Figure 7Mean prediction error (Eq. 37) of the recovery when deriving χ(t) for different values of the nonlinearity factor a of the toy model. As a increases, the recovery of χ(t) deteriorates because the level of the contributions from nonlinearities $\left|\right|\stackrel{\mathrm{̃}}{\mathit{\eta }}\left|\right|$ becomes large compared to the noise level $\left|\right|\mathit{\eta }\left|\right|$; how these terms are computed for the toy model is explained in Appendix D. To demonstrate here the pure effect from the breakdown of the spectral similarity assumption, the RFI algorithm is used here without the additional noise-level adjustment enforcing monotonicity.

Download

In Fig. 7 we show by plotting the mean prediction error (see Eq. 37) how the recovery of the response function deteriorates as the nonlinearity parameter a increases. To demonstrate that this is indeed caused by a breakdown of the spectral similarity assumption, we plot in addition the ratio $\left|\right|\stackrel{\mathrm{̃}}{\mathit{\eta }}\left|\right|/\left|\right|\mathit{\eta }\left|\right|$. It is seen that, indeed, as claimed above, the recovery works well only when this ratio is not large, i.e., when the contributions from nonlinearities $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ are not large compared to those from the noise η.

Figure 8Demonstration of how nonlinearities affect the recovery of the response function: (a) Picard plot; (b) coefficients of regularized solution (Eq. 20) and (c) recovered linear response function. First row: nonlinearity factor $a={\mathrm{10}}^{-\mathrm{10}}$ (no monotonicity check); second row: nonlinearity factor $a=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ (no monotonicity check); third row: nonlinearity factor $a=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ (with monotonicity check). The noise is overestimated in the low-frequency spectrum in the third row because nonlinearities yield a derived χ(t) that does not obey the monotonicity constraint. As a consequence, the method increases the level of low-frequency components until the monotonicity constraint is obeyed. The failure to obey the monotonicity constraint and consequent large overestimation of noise in this case can be taken as an indication of the presence of nonlinearities in the response. Note that the “true” linear response function in this nonlinear case a≠0 is obtained analytically from the linear case a=0 via Eq. (41) (see Appendix D). The regularization parameter determined by the algorithm is λ≈3120 for the first row, λ≈74 for the second, and λ≈14 611 873 for the third. For more details, see the text.

Download

More insight into how nonlinearities affect the recovery is obtained from the more detailed SVD analysis shown in Fig. 8. The first row of subfigures was obtained from the toy model assuming a rather small nonlinearity ($a={\mathrm{10}}^{-\mathrm{10}}$). In the Picard plot (Fig. 8a1) it is seen that in this case both conditions necessary for a good recovery are met: first, the signal $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$ is clearly visible above the combined noise $|{\mathit{u}}_i•(\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left)\right|$ (see the first four components). Second, in this case $|{\mathit{u}}_i•\stackrel{\mathrm{̃}}{\mathit{\eta }}|/|{\mathit{u}}_i•\mathit{\eta }|$ is small over the whole spectrum; i.e., the contributions from $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ are small compared to those from η. As explained above, because this second condition is also met, the noise estimate from the RFI algorithm η_est is a good approximation to the combined noise across all frequencies (compare in the Picard plot $|{\mathit{u}}_i•(\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left)\right|$ to $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$). As a result, the four components selected by the regularization for the recovered solution (Fig. 8b1) are precisely those dominated by the signal (compare $f_i\left(\mathit{\lambda }\right)|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|/{\mathit{\sigma }}_i$ with $f_i\left(\mathit{\lambda }\right)|{\mathit{u}}_i•(\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left)\right|/{\mathit{\sigma }}_i$). This example demonstrates that as long as these two conditions are met, small contributions from nonlinearities do not prevent a good recovery of the response function (see Fig. 8c1).

In the second row of Fig. 8, we demonstrate how the violation of the second condition obstructs the recovery. In this case the nonlinearity parameter has been given a larger value ($a=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}$). As a consequence, one sees in the Picard plot that the low-frequency components of the combined noise are enhanced. The first condition is still met: the signal $|{\mathit{u}}_i•\mathbf{\Delta }\mathit{Y}|$ is visible above the combined noise $|{\mathit{u}}_i•(\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left)\right|$ (see the first two components). However, now the ratio $|{\mathit{u}}_i•\stackrel{\mathrm{̃}}{\mathit{\eta }}|/|{\mathit{u}}_i•\mathit{\eta }|$ becomes large at low frequencies, violating the second condition. As explained, the violation of the second condition leads to the breakdown of the spectral similarity assumption. As a result, the RFI algorithm underestimates the combined noise at low frequencies (compare in the Picard plot $|{\mathit{u}}_i•(\mathit{\eta }+\stackrel{\mathrm{̃}}{\mathit{\eta }}\left)\right|$ to $|{\mathit{u}}_i•{\mathit{\eta }}_{\text{est}}|$). Using this wrong noise estimate, regularization selects components for the recovered solution that are to a large extent dominated by the combined noise (see components i=2 to i=6 in Fig. 8b2). The result is that the strong low-frequency contributions from nonlinearities deteriorate the recovery of the response function at long timescales (Fig. 8c2).

In the third row, we demonstrate for this type of nonlinearity that by accounting for monotonicity one can remove from the recovered solution all components dominated by noise. For this purpose, we set the nonlinearity parameter to the same value as for the second row ($a=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}$) but employ the additional noise-level adjustment (step 6 of Fig. 1); i.e., the low-frequency range of the noise estimate is now automatically adjusted in order to recover a response function that decays monotonically to zero. As seen in the Picard plot, the additional noise-level adjustment results in an artificial enhancement of the low-frequency components of the noise estimate, with a large jump separating the low- from high-frequency ranges. In this case, such enhancement is able to better estimate the largest components of the combined noise (first few components in the Picard plot). As a consequence, regularization correctly selects for the recovered solution only the two first components which are not dominated by noise (Fig. 8b3). Unfortunately, as seen in Fig. 8c3, these two first components do not contain enough information for a perfect recovery, since the quality improves at long timescales but deteriorates at short timescales (compare Fig. 8c3 and c2). This is a consequence of how regularization works: it filters out components dominated by noise (or in this case nonlinearity) at the expense of also removing useful information contained in those components.