2.1 RFI method and C⁴MIP experiments

The RFI method identifies the response function χ(t) taking data from the response Y(t) – in our example the global land carbon – and the perturbation f(t) – atmospheric CO₂ or temperature (see below) – assuming the following relation (see Part 1):

where η(t) is a noise term. Here ΔY(t) and Δf(t) mean that we are taking only the change in the variables with respect to their equilibrium values from a control simulation. Following Forney and Rothman (2012a), the spectrum of internal timescales is obtained by assuming that the response function χ(t) can be represented by an overlay of exponential modes:

To account for the large range of timescales in the carbon cycle (Ciais et al., 2013) it is useful to rewrite Eq. (4) in terms of log ₁₀τ, so that

The spectrum of timescales is then given by q(τ), which following the terminology from Part 1 we call simply spectrum. The problem is solved by discretizing Eqs. (3) and (5), prescribing a distribution of timescales τ, taking the data on ΔY(t) and Δf(t), and solving a minimization procedure for the spectrum q(τ). The parameters to prescribe the distribution of timescales are taken identically to those chosen for the application to the toy model in Part 1. To treat the ill-posedness we employ Tikhonov–Phillips regularization (Phillips, 1962; Tikhonov, 1963) in a singular value decomposition (SVD) framework that gives the solution by the expansion (Hansen, 2010; Bertero, 1989)

where • is the usual scalar product, M is the number of timescales, u_i and v_i are the singular vectors, σ_i are the singular values, λ is the regularization parameter, and f_i(λ) are the filter functions

The regularization parameter λ is determined by the discrepancy method (Morozov, 1966) with noise level estimated from a SVD analysis of the data combined with information from the control simulation. Once the spectrum q_λ is found, the response function χ(t) follows from Eq. (5).

All linear response functions are identified by the RFI method taking data from C⁴MIP-type experiments performed with the MPI-ESM. We focus on identifying the response functions from standard 1 % experiments that are widely available in international databases. In addition, to examine the quality of the results, we also identify some response functions taking data from additional experiments. To investigate the robustness of the identified response functions, we employ them to predict the response of the MPI-ESM in several experiments not used for the identification. A summary of the experiments considered in the study is given in Table 1, with forcing scenarios shown in Fig. 1.

The variables taken for the identification of the response functions are the ones relevant for the quantification of the land carbon sensitivities γ and β, respectively:

• (a)

  the change in global land carbon in response to changes in global land temperature;

• (b)

  the change in global land carbon in response to changes in atmospheric CO₂.

Global land carbon is computed as the sum of the total land carbon over all grid cells of the model. Global land temperature is calculated as the mean near-surface air temperature over land at 2 m height. The changes are computed as $\mathrm{\Delta }Y\left(t\right)=Y\left(t\right)-Y_{\text{PI}}$, with Y_PI being the mean value of observable Y from a control simulation at pre-industrial conditions. Since the main interest lies in long-term variations, annual mean data are used.

As demonstrated in Appendix C, the generalization of the β sensitivity can be shown to be monotonic. Therefore, in the following we will derive it employing the additional noise-level adjustment in the RFI algorithm (step 6 of Fig. 1 in Part 1). Since the generalization of the γ sensitivity is not known to be monotonic, for this sensitivity the RFI algorithm will be applied without the additional adjustment.

Table 1C⁴MIP-type experiments considered in this study. Forcings are shown in Fig. 1. Abbreviations “rad” and “bgc” refer to standard CMIP model setups used to calculate the climate-carbon cycle sensitivities. In the “rad” (radiatively coupled) setup only the radiation code of the model is affected by changes in atmospheric CO₂. This setup is used to calculate γ. In the “bgc” (biogeochemically coupled) setup only the biogeochemistry of the model is affected by changes in atmospheric CO₂. This setup is used to calculate β. In brackets are names of standard CMIP experiments.

Download Print Version | Download XLSX

Figure 1Forcings for the C⁴MIP-type experiments considered in this study.

Download

2.2 Estimating the linear regime from “percent” experiments

As described in Part 1, the recovery of response functions is complicated by the presence of noise and nonlinearities. The RFI method is designed to cope with the former. In the present section we present a technique to cope with the latter, but at the expense of performing additional response experiments. This technique will serve as a complement to the RFI method in the application to the land carbon cycle in the following sections. By these additional experiments we will determine the range of forcing strengths for which the response can be considered linear – an analysis in general not possible using only the control and perturbation experiments on which the RFI method is based.

Before we introduce this technique to cope with nonlinearities, it is useful to specify more clearly what we mean by “nonlinearity”. Dynamical systems like the Earth system are crowded with nonlinearities. Our notion of “nonlinearity” is much more specific: the linear response formula (Eq. 3) can be interpreted as the linear term in a Volterra expansion of the response ΔY(t) of the system into the perturbation Δf(t) (see Sect. 2 in Part 1). This expansion is given by

where 𝒪((Δf)²) represents terms of order 2 and higher in the perturbation Δf. For small enough perturbations the term 𝒪((Δf)²) is small compared to the linear-order term and can be ignored, so that Eq. (8) yields Eq. (3). Throughout this study, nonlinearities are defined as the nonlinear terms 𝒪((Δf)²) in this expansion. Indeed, these nonlinear terms arise from the intrinsic nonlinearities of the underlying system. In the language of statistical mechanics one would call those intrinsic nonlinearities as “microscopic”. This in mind, our notions of linearity and nonlinearity should then be termed “macroscopic”. Since the idea here is to study the linear response of the system, we will be interested only in the case where perturbations are so small that the nonlinear terms 𝒪((Δf)²) can be ignored. However, the question is whether they can be ignored. It is clear that nonlinearities may gain importance with increasing forcing strength. Hence, to study only the linear response, experiment data can be used for our purpose only up to such perturbation strengths where nonlinearity start to contribute to the response; i.e., we have to identify the regime where the response is still linear. This is the purpose of the technique that we describe in the following.

To introduce our technique we use the example of simulations with the linear toy model employed in Part 1. To demonstrate the effect of nonlinearities on the recovery of χ(t), following Part 1 we artificially add a nonlinear term −aY²(t) to the linear response Y(t) of the toy model:

In this way we obtain a nonlinear response Y_nonlin(t), with nonlinearity strength controlled by the parameter a, from which χ(t) is derived. In addition to including this nonlinear term in the toy model response, to introduce the technique we will need to quantify the quality of the recovery of the response function. Since in our application to the land carbon the “true” response function is not known a priori, following Part 1 we quantify the quality of the recovery indirectly by measuring the quality with which the recovered response function can predict the response of the model in experiments not used for the recovery itself. For this purpose we introduce the relative prediction error

where ⋆ stands for the convolution operation, ΔY^k and Δf^k are the response and the perturbation in experiment “k”, and χ is the response function recovered from the 1 % experiment.

We can now present our technique. Taking first a purely linear situation (a=0), we show in Fig. 2a the relative prediction error (Eq. 10) when using the response function obtained from a 1 % simulation to predict the response from two other % simulations with smaller growth rate. More precisely, performing a sequence of 1 % experiments for increasingly longer simulation periods, we calculated for each experiment the response function and used it to predict the response for a 0.5 % and 0.75 % experiment covering the same simulation period. Then we plotted in Fig. 2a the relative prediction error against the final forcing strength of the 1 % experiment. As a result, it is seen that with increasing final forcing strength the relative prediction error decreases. This happens because in this linear case the SNR is increasing with increasing simulation period, i.e., with increasing final forcing strength, so that the recovery of the response function continuously improves.

Figure 2Toy model example for the identification of the linear regime by using additional experiments. Shown is the relative prediction error (Eq. 10) for the response of 0.5 % and 0.75 % experiments as obtained from the response function calculated by the RFI method from 1 % experiments. The relative prediction errors are plotted against the final forcing strengths of a sequence of 1 % experiments with increasing time series length. The crosses at the minima indicate the final forcing strength for which the response function is optimally recovered (see text). (a) shows the behavior for the fully linear toy model (a=0) and (b) the behavior in the presence of nonlinear contributions to the response ($a=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$). For the purpose of demonstrating more clearly the increase in the relative prediction error for a decrease in the forcing strength, we include in the plot cases where the forcing strength is extremely small, corresponding to very small time series lengths. To deal with such cases, we set for a number of data points N<30 the number of timescales M=N (see parameters for the RFI method in Part 1). For such small number of timescales, usually no plateau in the singular values spectrum is found (step 2 of Fig. 1 in Part 1). Therefore, for these special cases we also modify the algorithm to interpret the two smallest singular values as a plateau, since their small magnitude makes them have a similar effect to those singular values belonging to the plateau itself. In addition, to illustrate the most general case where χ(t) is not known to be monotonic, we exclude here the monotonicity check (step 6 of Fig. 1 in Part 1).

Download

Calculating the relative prediction error only for experiments with smaller growth rate becomes important in the next case where nonlinearities are considered (Fig. 2b). This plot was obtained by the same procedure except that we took for the nonlinearity parameter a value a>0. As seen, in this nonlinear case the relative prediction error is first improving but deteriorating afterwards. For small forcing, nonlinearities are small and therefore the relative prediction error behaves as in the linear case; i.e., it decreases with final forcing strength. However, when forcing strengths become larger, nonlinearities start to contribute substantially to the response, thereby causing a deterioration of the recovery of the response function and consequently the relative prediction error once more increases. This increase of the prediction error can be unambiguously traced back to the presence of nonlinearities in the 1 % simulation because the prediction error was calculated only for experiments with smaller growth rate, i.e., smaller forcing strength throughout the whole simulation. Therefore, nonlinearities contribute already substantially to the response in the 1 % simulation before they become relevant in the other experiments. Accordingly, with this type of experiment setup we can be sure that the increase in the prediction error comes solely from the deterioration of the recovery of the response function and not from nonlinearities in the additional experiments used for calculating the prediction error.

Obviously, for forcing strengths smaller than at the minimum, nonlinearities do not hinder the recovery of the response function, so that one can consider this to be the regime of linear system behavior. In view of the trade-off between noise and nonlinearities, for the 1 % experiment the response function is thus optimally recovered when taking as final forcing strength the value at the minimum of the prediction error curve. Similarly, if the error curve has no minimum (as in the linear case shown in Fig. 2a), the optimal recovery is obtained from the experiment with the largest forcing strength.

With the presentation of this additional technique to identify the linear regime we are ready for the application to the MPI-ESM in the next sections.

3.1 Estimating the linear regime

As a first step in obtaining a proper approximation of χ_γ(t), we investigate what maximum forcing strength can be used to ensure that the recovered χ_γ(t) is not spoiled by the presence of nonlinearities. Using the technique introduced in Sect. 2.2, we show in Fig. 3 the relative prediction error (Eq. 10) for χ_γ(t) recovered from the 1 % rad experiment as a function of the final forcing strength in the 1 % rad experiment. There is no clear minimum, so that for the data available the recovery of χ_γ(t) seems not to be limited by nonlinearities. For optimal recovery we thus take the full time series, i.e., the maximum final forcing strength.

Figure 3Relative prediction error (Eq. 10) for the 0.5 % and 0.75 % rad experiments using χ_γ(t) recovered from the 1 % rad experiment. The error is shown as function of the final forcing strength used for the recovery of χ_γ(t). No clear minimum is found, so that the recovery seems not to be limited by nonlinearities.

Download

3.2 χ_γ and the quality of its recovery

The quality of the recovery can in principle be improved by taking an experiment with better SNR. To investigate whether the recovery from the 1 % rad experiment can be further improved, we also applied the RFI algorithm to recover χ_γ from the 2×CO₂ rad experiment. We chose the 2×CO₂ rad experiment because as shown in Fig. 4 it has smaller forcing strengths than the maximum forcing strength for the 1 % rad experiment – therefore nonlinearities should also not limit the recovery – but is expected to carry useful information over a larger range of the response spectrum. This expectation can be justified as follows (MacMartin and Kravitz, 2016). Taking the Laplace transform of Eq. (11) gives

where the tilde denotes Laplace transformed functions. From Fig. 4 it is seen that for the 1 % rad experiment the temperature ΔT behaves approximately as a linear function, which gives a Laplace transform $\mathrm{\Delta }\stackrel{\mathrm{̃}}{T}\left(p\right)$ proportional to $\mathrm{1}/p^{\mathrm{2}}$. For the 2×CO₂ rad experiment, the temperature behaves approximately as a step function (ignoring the transient in the first 20 years), which gives a Laplace transform $\mathrm{\Delta }\stackrel{\mathrm{̃}}{T}\left(p\right)$ proportional to $\mathrm{1}/p$. This means that for the same χ_γ, the first term on the right-hand side of Eq. (13) – the “clean” response – decays to zero faster for the 1 % rad experiment than for the 2×CO₂ rad experiment. Hence, assuming the same noise η for both cases, the response from the 2×CO₂ rad experiment is buried in the noise only at larger p, meaning that this experiment carries useful information until higher rate values p than the response from the 1 % rad experiment.

Figure 4Forcing temperature ΔT(t) for 1 % and 2×CO₂ rad experiments.

Download

The response function χ_γ(t) recovered from the two types of experiments is shown in Fig. 5a. As expected, the different experiments indeed result in different recoveries. Because we know from the analysis of Fig. 3 that nonlinearities do not limit the recovery, the difference between the two response functions probably results from the difference in the quality of the data from the two types of experiments. To compare the robustness of each recovery, we analyze how well they predict the response from other experiments (Fig. 5b and c). If the response function is correctly recovered, it should be able to predict not only experiments with smaller, but also experiments with higher forcing rates. Therefore, we also include in the analysis 1.5 % and 2 % rad experiments. To exclude errors that may be caused by nonlinearity, we take as a conservative estimate of the linear regime forcing strengths smaller than the final forcing strength at the end of the 1 % rad experiment (which is approximately the maximum forcing strength; see temperature value at t=140 years for the 1 % rad experiment in Fig. 4). We take these values as an estimate of the linear regime because the 1 % rad experiment has only 140 years, so that no estimate for higher forcing strengths is available. Hence, for the 1.5 % and 2 % rad experiments the responses are expected to be reasonably predictable at least until the values marked with circles, where their respective forcing strengths reach this maximum forcing strength. All other experiments have forcing strengths smaller or equal to this maximum forcing strength, so that they should be predictable for the whole time series.

Figure 5χ_γ recovered from 1 % and 2×CO₂ rad experiments (a) and prediction of model responses using these recoveries of χ_γ. Circles indicate the maximum value for which 1.5 % and 2 % responses are predictable according to the estimate of the linear regime (see text). At the right of (b, c) the relative prediction error (see Eq. 10) is indicated for the different experiments, calculated for the 1.5 % and 2 % rad experiments by considering only values preceding the circles.

Download

Figure 5b shows the quality of the prediction using χ_γ(t) recovered from the 1 % rad experiment. Visually, model response and prediction seem to have a comparable quality of agreement across the 1.1×CO₂, 0.5 %, 0.75 % and 1.5 % rad experiments, while for the 2 % and 2×CO₂ rad experiments there are larger discrepancies. For a quantitative analysis, we compute for the estimated linear regime the relative prediction error (Eq. 10) for each experiment (right side of the plot). It is seen that the error varies from less than 10 % for the 0.75 % and 1.5 % rad experiments to values between 10 % and 20 % for the 0.5 %, 2 % and 2×CO₂ rad experiments, and a significantly larger value of 57 % for the 1.1×CO₂ rad experiment. To better understand these differences, it is important to note that as long as nonlinearities are small, experiments with higher forcing strength are expected to have smaller prediction error because they have higher SNR. This can be made plausible by considering that a perfectly recovered response function predicts a “clean” linear response (infinite SNR) with zero error, whereas the same response function can predict a noisy response only with some finite error. Therefore, if χ_γ(t) is well recovered, we expect large prediction errors for experiments with small forcing strengths such as the 1.1×CO₂ rad – which is indeed the case – but small errors for experiments with large forcing strengths but still well within the linear regime such as the 2×CO₂ rad (compare the forcing strengths for the 2×CO₂ rad experiment and the maximum forcing strength for the 1 % rad experiment in Fig. 4). Since contrarily to the expectation the prediction error is not particularly small for the 2×CO₂ rad experiment, probably the recovery of χ_γ(t) derived from the 1 % rad experiment is not completely accurate and may still be further improved. As suggested above, such improvement may be achieved by taking data with better quality from the 2×CO₂ rad experiment.

Figure 5c shows the quality of the prediction using χ_γ(t) recovered from the 2×CO₂ rad experiment. As expected, results indicate an improvement in the recovery (compare to Fig. 5b). The prediction error decreases for all experiments present in both plots. In addition, it also decreases if we compare the prediction of the 1 % rad response in Fig. 5c with that of the 2×CO₂ rad response in Fig. 5b.

This section therefore suggests two main conclusions: first, for χ_γ(t) the response seems to be approximately linear for temperature perturbations up to at least 6 K. Second, the overall improvement of the prediction in Fig. 5c compared to Fig. 5b confirms the expectation from the Laplace transform analysis that the step-like 2×CO₂ rad experiment indeed carries more information on the response function than the 1 % rad experiment. This suggests that step-like experiments may be more appropriate than the standard 1 % experiment for the identification of linear response functions.

4.1 Approaches to identify χ_β

Similarly to the last section, we identify χ_β(t) by several approaches to find the one that gives results with the best quality. For this purpose, we also consider in addition to Eq. (14) alternative formulas to derive χ_β(t), each taking a different forcing. The identification is performed in three different ways:

1. Using CO₂ as forcing (see Eq. 14);

2. Using the logarithm of CO₂ as forcing:

   $$\begin{array}{}\text{(16)}& \begin{array}{rl}& \mathrm{\Delta }{C}^{\text{bgc}}\left(t\right)\\ & =\underset{\mathrm{0}}{\overset{t}{\int }}{\mathit{\chi }}_{\mathit{\beta }}^{\ln }(t-s)C_{\text{atm}}^{\mathrm{0}}\ln \left({\displaystyle \frac{C_{\text{atm}}\left(s\right)}{C_{\text{atm}}^{\mathrm{0}}}}\right)\mathrm{d}s+\mathit{\eta }\left(t\right),\end{array}\end{array}$$

   where now the forcing is

   $$\mathrm{\Delta }f\left(t\right):=C_{\text{atm}}^{\mathrm{0}}\ln \left({\displaystyle \frac{C_{\text{atm}}\left(t\right)}{C_{\text{atm}}^{\mathrm{0}}}}\right),$$

   with $C_{\text{atm}}^{\mathrm{0}}$ being the pre-industrial value for atmospheric CO₂;

3. Using net primary production (NPP) as forcing:

   $$\begin{array}{}\text{(17)}& \begin{array}{rl}& \mathrm{\Delta }{C}^{\text{bgc}}\left(t\right)\\ & =\underset{\mathrm{0}}{\overset{t}{\int }}{\mathit{\chi }}_{\text{NPP}}(t-s)\mathrm{\Delta }\text{NPP}\left(C_{\text{atm}}\right(s\left)\right)\mathrm{d}s+\mathit{\eta }\left(t\right).\end{array}\end{array}$$

The first approach is the same used throughout the paper: χ_β(t) is identified using the ΔC_atm(t) forcing from Eq. (14).

The second and third approaches are motivated as follows. In the “bgc” setup, changes of atmospheric CO₂ affect land carbon exclusively via the resulting changes in primary productivity¹. Therefore, the reaction of land carbon to CO₂ perturbations may be understood as a two-step process: first, CO₂ affects primary productivity (carbon assimilation); second, the assimilated carbon adds to land carbon storage. This second step is in models typically described by an almost linear system of storage compartment forced by NPP (see, e.g., Xia et al., 2013; Huang et al., 2018). Thus, presumably the response of land carbon to NPP is well described by a linear response function. These considerations underlie the second and third approaches but are more explicitly realized in the latter: there we first derive the linear response function χ_NPP(t) describing the response of land carbon to changes in NPP, and then use the data of the simulated dependence of NPP on CO₂ to derive the desired response function χ_β(t) (see Eq. 21). In this way we try to keep the most critical step in the identification of the response function, namely the inversion of the response integral, largely free from nonlinearities, and indeed, for this reduced problem, without much interference from nonlinearity we can recover the response function from experiment data obtained for much higher perturbation strengths than in the first approach, and because of the higher SNR the quality of recovery is thereby largely improved (see below). Instead of NPP one could also try to split off the nonlinear part by taking GPP (gross primary production), which differs from NPP only by the autotrophic respiration. However, autotrophic respiration is lost to the atmosphere, so that it is NPP and not GPP that adds to land carbon storage. Thus, by taking NPP instead of GPP the linear part is arguably splitted off more completely. The second approach is motivated in the same way, except that the nonlinear part is guessed to be a logarithmic function of CO₂ as suggested by the Keeling formula (see, e.g., Alexandrov et al., 2003), so that here no additional NPP data are needed from simulations.

Mathematically, our procedure in each approach is the following. In the second approach (see Eq. 16), we employ as forcing a logarithmic expression because it is known that a large part of the nonlinear relation between CO₂ and carbon storage comes from the approximately logarithmic dependence of NPP on CO₂ (e.g., Alexandrov et al., 2003). Such formula has the advantage that the expansion of the forcing in C_atm gives

Taking ΔC_atm sufficiently small and comparing the result to Eq. (14) thus yields

Accordingly, the response function ${\mathit{\chi }}_{\mathit{\beta }}^{\ln }\left(t\right)$ from Eq. (16) also gives the desired χ_β(t).

In the third approach, we take directly the response to NPP (see Eq. 17). Expanding the forcing ΔNPP(C_atm) in C_atm gives

Taking ΔC_atm sufficiently small and comparing the result to Eq. (14) yields

Accordingly, by this approach we compute χ_β(t) in three steps: first, we identify the response function χ_NPP(t) using Eq. (17); second, we take the first derivative of NPP with respect to CO₂ at $C_{\text{atm}}=C_{\text{atm}}^{\mathrm{0}}$; third, we apply Eq. (21) to obtain χ_β(t) from χ_NPP(t).

4.2 Checking nonlinearities with the three approaches

Before analyzing the recovery of χ_β(t) employing the three approaches, one must verify that these two additional approaches indeed account for some of the nonlinearities in the response. If this is true, response formulas (Eqs. 16 and 17) should be able to predict the response to larger perturbation strengths than Eq. (14). To verify this expectation, in the following we compare the relative prediction error (Eq. 10) by applying Eqs. (14), (16) and (17).

Figure 6a shows the relative prediction error for χ_β(t) recovered from the 1 % bgc experiment with the first approach (Eq. 14). The prediction is also computed via Eq. (14). The clear minima indicate the presence of strong nonlinearities for forcing strengths above around 100 ppm (94 ppm for the 0.5 % and 133 ppm for the 0.75 %). Therefore, in contrast to the case of χ_γ discussed in the last section, we see here that one indeed has to cope with the additional difficulty of nonlinearities.

Figure 6Relative prediction error (Eq. 10) for the 0.5 % and 0.75 % bgc experiments obtained when using χ_β(t), ${\mathit{\chi }}_{\mathit{\beta }}^{\ln }\left(t\right)$ and χ_NPP(t) obtained from the 1 % bgc experiment to predict the response. The error is shown as a function of the CO₂ final forcing strength.

Download

Figure 6b shows the prediction error when using ${\mathit{\chi }}_{\mathit{\beta }}^{\ln }\left(t\right)$ recovered from the 1 % bgc experiment with the second approach (Eq. 16). To check how well nonlinearities are accounted for by taking the logarithmic forcing, the prediction is also computed via Eq. (16). Compared to Fig. 6a, we see a slight improvement in the results: the minima have smaller prediction errors and the prediction errors increase at a slower rate for increasing final forcing strength. This indicates that using the logarithm of CO₂ as forcing indeed accounts for some of the nonlinearities in the response. Accordingly, one can make predictions with smaller error for larger forcing strengths using Eq. (16) instead of Eq. (14).

Figure 6c shows the prediction error for χ_NPP(t) recovered via Eq. (17) from the 1 % bgc experiment (first step of the third approach). To check how well nonlinearities are accounted for by taking the NPP forcing, we also employ Eq. (17) for the prediction. Here, we see a substantial improvement in the results. The response is almost completely linear, with very “flat” minima. This indicates that a large part of the nonlinearity encountered in the response of the land carbon to changes in CO₂ can indeed be explained by the nonlinear relationship between NPP and CO₂. Accordingly, by employing Eq. (17) instead of Eq. (14) one can predict the response of the land carbon until forcing strengths as high as 800 ppm with an error smaller than 10 %.

4.3 χ_β and the quality of its recovery

So far, we have only considered the ability of Eqs. (14), (16) and (17) to predict the land carbon response. Now, we analyze how well the generalized sensitivity χ_β(t) can be identified by the three approaches. For the identification we took data from the 1 % bgc experiment until the CO₂ forcing strength reaches the first minimum for each case in Fig. 6: ΔC_atm=94 ppm for the first approach (30 years); ΔC_atm=133 ppm for the second approach (40 years); and ΔC_atm=279 ppm for the third approach (70 years). Since ΔC_atm=279 ppm is approximately the forcing strength for the 2×CO₂ bgc experiment and results from last section suggested that this type of experiment may carry more information for the identification, we employ the third approach, also taking the 2×CO₂ bgc experiment. For the present application where the recovery is limited by nonlinearities, taking the 2×CO₂ bgc experiment has the additional advantage that because its forcing strength has a constant value throughout the whole experiment, we can use the full time series (140 years). To compute the first derivative of NPP with respect to CO₂ (second step of the third approach), we fitted the function NPP=NPP(C_atm) to polynomials of order 4, 5, and 6 and then took the first derivative from the fits.

The results from the three approaches are shown in Fig. 7a. At short timescales there is an overall agreement among all recoveries with only small discrepancies. To be able to compare the results also for longer timescales, we extend the response functions recovered from the 1 % bgc experiment – obtained from time series with 30, 40, and 70 years, respectively, for the first, second, and third approaches – until 140 years (extensions are indicated by dotted lines). This can be done because with the RFI method we derive the response function from the ansatz (5), which formally gives the values of the response function for all times. The result is that all response functions recovered from the 1 % bgc experiment present relatively small discrepancies even at long timescales. Response functions derived from the 2×CO₂ bgc experiment with the third approach show a similar behavior among themselves but differ from the recoveries using the 1 % bgc experiment. The reason for this difference will be investigated below.

Figure 7Response function χ_β(t) derived by the three approaches (a) and the respective prediction errors for the first 30 years of the response (b). $\mathrm{\Delta }{C}_{\text{atm}}^{\mathrm{1}\mathit{\%}}$: recovery with first approach from 1 % bgc experiment; log (C_atm)^1 %: recovery with second approach from the 1 % bgc experiment; ${\text{NPP}}_{\text{pol}x}^{\mathrm{1}\mathit{\%}}$: recovery with third approach from the 1 % bgc experiment using for the derivative a polynomial fit of order x; ${\text{NPP}}_{\text{pol}x}^{\mathrm{2}\times }$: recovery with third approach from the 2×CO₂ bgc experiment using for the derivative a polynomial fit of order x. Continuous lines denote values of χ(t) that are within the time series length used for the recovery (30 years for $\mathrm{\Delta }{C}_{\text{atm}}^{\mathrm{1}\mathit{\%}}$, 40 years for log (C_atm)^1 % and 70 years for ${\text{NPP}}_{\text{pol}x}^{\mathrm{1}\mathit{\%}}$). Dotted lines denote extended parts of the response function, i.e., values not covered by the time series used for the recovery but obtained from the recovered spectrum using Eq. (5). Circles denote the response functions derived taking the full time series (${\text{NPP}}_{\text{pol}x}^{\mathrm{2}\times }$). For more details, see text.

Download

To quantitatively compare the quality of the recoveries, we plotted in Fig. 7b the relative prediction error (Eq. 10). Since the response functions were derived using different time series lengths, for a fair comparison we compute the error only at the minimal time series length of 30 years (the time series length used for the first approach). Results show no large discrepancies among the different approaches. For the third approach, there seems to be a small advantage in using a polynomial of order 6 for the computation of the derivative.

However, results from Fig. 7b reflect only the quality of the recoveries at short timescales, for which anyway no large discrepancies were encountered in Fig. 7a. To evaluate the quality of the recovery also at long timescales, one must take the extended version of the response functions that were derived from shorter time series. Since the only substantial difference at long timescales in Fig. 7a is found between the response functions recovered from the 1 % and 2×CO₂ bgc experiments, we take for this analysis exemplarily only the response functions recovered with the first approach (1 % bgc experiment) and the third approach (2×CO₂ bgc experiment). By choosing the response function recovered with the first approach, we evaluate for the worst case scenario (where only 30 years are used for the recovery) how reliable predictions are for timescales longer than the time series used for recovery. In contrast, by choosing the response function recovered with the third approach from the 2×CO₂ bgc experiment, we check whether the different values at long timescales are actually an improvement in the recovery. As mentioned above, such improvement is expected because this response function was recovered taking the full time series (140 years).

Following the same procedure as in the last section, in Fig. 8 we show the quality of the prediction for different experiments using the aforementioned recoveries of χ_β(t). Figure 8a shows the results for the predictions calculated via Eq. (14) using χ_β(t) recovered with the first approach. We take as an estimate of the linear regime forcing strengths smaller than the forcing strength at the first minimum in Fig. 6a. The response values corresponding to this forcing strength are marked in Fig. 8a with circles. It is seen that although χ_β(t) was recovered using a time series of only 30 years, it predicts the 1.1×CO₂ bgc experiment with only 3 % error over 140 years. Other experiments are predicted within the linear regime with error smaller than 10 % with exception of the 2 % bgc experiment, for which the error is around 14 %. Since the 2×CO₂ has a constant forcing strength outside the linear regime already from the beginning, its prediction fails as expected for the whole time series.

Figure 8Prediction of model responses employing response functions derived with the first approach from the 1 % bgc experiment (derived from data with 30 years length, but extended to 140 years) and with the third approach from the 2×CO₂ bgc experiment (derived from data with 140 years length). (a) Prediction by Eq. (14) taking the response function χ_β(t) derived from the 1 % bgc experiment with the first approach. (b) Prediction by Eq. (17) taking the response function χ_β(t) derived from the 1 % bgc experiment with the first approach and converted to χ_NPP(t) by Eq. (21). (c) Prediction by Eq. (17) taking the response function χ_NPP(t) derived in the first step of the third approach (see explanation after Eq. 21) taking data from the 2×CO₂ bgc experiment. Continuous lines are predictions and dashed lines are responses from the MPI-ESM. Circles indicate the maximum value for which responses are predictable according to our estimate of the linear regime (see text). The values printed to the right of the plots are the relative prediction errors (see Eq. 10) calculated for each experiment, considering when applicable only values preceding the circles.

Download

In Fig. 8a we could only evaluate the quality of the long timescales of χ_β(t) by the prediction of the 1.1×CO₂ bgc experiment, because this is the only experiment which has forcing strengths within the linear regime over the whole time series. To check the ability of χ_β(t) to predict also other experiments at long timescales, we account for some of the nonlinearity in the response by taking NPP instead of CO₂ as forcing (see discussion of Fig. 6). Therefore, we perform the prediction by employing Eq. (17) instead of Eq. (14). Since for employing Eq. (17) one needs χ_NPP(t) and not χ_β(t), we take the χ_β(t) derived with the first approach and compute χ_NPP(t) from it via Eq. (21). Because the conversion from χ_β(t) to χ_NPP(t) is a simple scaling, the timescales structure is maintained. Hence, we can evaluate the quality of the recovered χ_β(t) from the results given by the obtained χ_NPP(t). The prediction results computed via Eq. (17) with the obtained χ_NPP(t) are shown in Fig. 8b. Because errors at the minima in Fig. 6c are not substantially smaller than those at the maximum final forcing strength, we take as an estimate of the linear regime all values smaller than the last value of NPP for the 1 % bgc experiment (marked with circles in the responses). Once again it is seen that although χ_β(t) was recovered using a time series of 30 years, after conversion to χ_NPP(t) almost all experiments can be predicted with less than 10 % error for the whole time series. The 1.5 % and 2 % bgc experiments are predicted with errors of 17 % and 4 % within the linear regime. Results from Fig. 8a and b therefore suggest that although nonlinearities do restrict the recovery from Eq. (14), taking experiments with forcing strengths within the linear regime for the recovery leads to reliable prediction results even for times reasonably longer than the length of the time series from which the response function was recovered.

Finally, we investigate whether the different values seen in Fig. 7a for the χ_β(t) recovered from the 2×CO₂ bgc experiment indeed reflect a better quality of recovery. Following the same reasoning that led to Fig. 8b, since in the third approach χ_β(t) is obtained from a scaling of χ_NPP(t), we evaluate the quality of χ_β(t) from the results given by χ_NPP(t). Accordingly, in Fig. 8c we plot the prediction via Eq. (17) using the χ_NPP(t) recovered from the 2×CO₂ bgc experiment. As expected, the overall prediction indeed improves compared to Fig. 8b. Individual relative prediction errors decrease for all experiments, with the exception of the 1.1×CO₂. Since the response functions used for Fig. 8b and c differ only at long timescales, this improvement suggests that obtaining χ_β(t) from the 2×CO₂ bgc experiment indeed gives a better recovery over these timescales. Further, because all recoveries are similar at short timescales (see Fig. 7a), overall this recovery shows the best quality.

5.1 Obtaining the detailed spectrum

However, before we investigate the recovered spectrum, we demonstrate how such a detailed structure may arise from a better recovery of the response function. In Fig. 9a we plot the spectrum q(τ) of the response function used for Fig. 8b, i.e., χ_β(t) recovered with the first approach (see beginning of Sect. 4) and converted to χ_NPP(t) via Eq. (21). Because the data used for the recovery have a relatively low quality – the response function was recovered from the 1 % bgc experiment taken for only 30 years – regularization filters out most of the SVD components of the spectrum in Eq. (6). Since the low-index SVD components that are not filtered out tend to be smooth (Hansen, 1989, 1990, see also Part 1), the final result of this filtering is the smooth spectrum seen in Fig. 9a. Obviously, although this smooth spectrum is a sufficiently good approximation to make the predictions shown in Fig. 8b, it is not very informative about the internal timescale structure. Instead, the spectrum of the response function χ_NPP(t) used for Fig. 8c has a more detailed structure (Fig. 9b). In this case, the higher quality of the data used for the identification (2×CO₂ bgc experiment taken for 140 years) results in less filtering by regularization, thereby revealing more details of the underlying spectrum. The result is a spectrum with two peak timescales, at around 4 and 100 years².

5.2 Checking the robustness of the spectrum via a Gregory-type approach

Some trust in this result may be gained via a linear-regression procedure following the same logic underlying a “Gregory plot” (Gregory et al., 2004) but applied here not to radiative forcing but to land carbon. As for a Gregory plot, we look here at how the nonequilibrium fluxes relax to equilibrium after a step-type perturbation. In our case this is the 2×CO₂ bgc experiment. The idea here is to try to identify from an independent method important timescales in the response, so that they can be compared against the spectrum in Fig. 9b. For this analysis, we plot the global land carbon against its first time derivative (this is the net land–atmosphere carbon flux). Using the 2×CO₂ bgc experiment where the forcing is constant, the first time derivative vanishes as the land carbon approaches a new equilibrium value. The rate at which the derivative changes can be associated with a timescale τ_i, which should show up with a large weight value q_i in the spectrum. Interestingly, the plot shows that the transient behavior towards equilibrium develops approximately at two different rates: a higher rate from the starting point until around 3520 Pg C, and a lower rate from 3520 Pg C onwards. To determine these rate values, we fitted a linear function $\mathrm{d}C/\mathrm{d}t=a+b{C}^{\text{bgc}}$ for each of these two ranges of the land carbon. Then, from each rate value b we computed a timescale by $\mathit{\tau }=-\mathrm{1}/b$. The computed timescales taking 1 standard deviation into account are shown by the two ranges highlighted in Fig. 9b. As seen, this linear-regression procedure also reveals a timescale structure dominated by two timescales. While the shortest timescale peak of the spectrum partially overlaps with the corresponding timescale range from the linear-regression plot, the longest timescale peak is almost perfectly matched³. This result suggests that the recovered spectrum indeed reflects internal characteristics of the global response of the land carbon cycle.

Figure 9Spectra associated with χ_NPP derived with different resolutions and a Gregory-type plot to identify the timescale structure of land carbon. (a) Spectrum derived with the first approach (taking CO₂ as forcing) from the 1 % bgc experiment; (b) spectrum derived with the third approach (taking NPP as forcing) from the 2×CO₂ bgc experiment; (c) Gregory-type plot of global land carbon against land–atmosphere carbon flux. Dots are the data, with the color scale from white to black indicating the evolution from 1 to 140 years. Values of b indicate the rate at which the time derivative of land carbon changes with respect to the land carbon itself. Ranges of timescales corresponding to each rate accounting for 1 standard deviation are shown in (b).

Download

5.3 Checking the robustness of the spectrum via regional response functions

The robustness of this detailed spectrum can be further checked by a different method. In the following, we test this robustness by checking the consistency between the timescales of global and regional carbon responses.

To study regional carbon responses, we considered two regions: tropics and extra-tropics. Tropics were defined as the region between latitudes 30^∘ S and 30^∘ N, and extra-tropics as the remaining part of the globe. We then determined separately the linear response functions that characterize the land carbon response to NPP in the tropics and extra-tropics, defined, respectively, by

The data were taken once more from the 2×CO₂ bgc experiment.

Figure 10Investigation of the land carbon response in the tropics and extra-tropics and how the regional response functions combine to the global response functions. The analysis is based on the 2×CO₂ bgc experiment. (a) Laplace transform ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{NPP}}\left(p\right)$ of global χ_NPP(t) obtained directly from the global carbon response and from combining the tropical and extra-tropical response functions; (b) χ_β(t) obtained directly from the global carbon response and from combining the tropical and extra-tropical response functions; (c) as (b) but for q_β(τ); (d) decomposition of q_β(τ) into tropical and extra-tropical spectra (Eq. 31). In (c, d) the dots and crosses indicate the computed values, while the connecting lines are only inserted to guide the eye. For more details, see text.

Download

Before assessing the robustness of the spectrum, we check the consistency between the regional and global response functions. In this test, we show that the global response function χ_NPP can be reconstructed by combining ${\mathit{\chi }}_{\text{NPP}}^{\text{tr}}$ and ${\mathit{\chi }}_{\text{NPP}}^{\text{et}}$. For this purpose, we write the global land carbon as the sum of the land carbon in the tropics and extra-tropics:

Plugging Eqs. (17), (22) and (23) in Eq. (24) and recognizing that $\mathit{\eta }\left(t\right)={\mathit{\eta }}^{\text{tr}}\left(t\right)+{\mathit{\eta }}^{\text{et}}\left(t\right)$ gives

Applying a Laplace transform to both sides of Eq. (25) and reorganizing the resulting equation gives

Therefore, the Laplace transformed ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{NPP}}\left(p\right)$ can be obtained from combining the NPP forcings and the response functions for the tropics and extra-tropics. Hence, if the response functions are correctly recovered by the RFI algorithm, Eq. (26) should hold at least approximately. In order to check this, we computed the Laplace transforms analytically by approximating the NPP forcings by step functions (since we take the 2×CO₂ bgc experiment), and using Eq. (5) for the response functions. Figure 10a shows the quality of the result. The small discrepancy between ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{NPP}}\left(p\right)$ obtained from the global carbon response and from combining the regional responses can be at least partially explained by the approximation error made to represent the forcings by a step function.

We now check the robustness of the land carbon spectrum by combining ${\mathit{\chi }}_{\text{NPP}}^{\text{tr}}\left(t\right)$ and ${\mathit{\chi }}_{\text{NPP}}^{\text{et}}\left(t\right)$ to obtain χ_β(t) and thereby q_β(τ). To this end, we first obtain the tropical ${\mathit{\chi }}_{\mathit{\beta }}^{\text{tr}}\left(t\right)$ and extra-tropical ${\mathit{\chi }}_{\mathit{\beta }}^{\text{et}}\left(t\right)$ by applying Eq. (21):

Note that because CO₂ is well mixed, $\mathrm{\Delta }{C}_{\text{atm}}^{\text{tr}}\left(t\right)\approx \mathrm{\Delta }{C}_{\text{atm}}^{\text{et}}\left(t\right)\approx \mathrm{\Delta }{C}_{\text{atm}}\left(t\right)$, so that, by characterizing the tropical and the extra-tropical response to CO₂, ${\mathit{\chi }}_{\mathit{\beta }}^{\text{tr}}\left(t\right)$ and ${\mathit{\chi }}_{\mathit{\beta }}^{\text{et}}\left(t\right)$ are describing the response to a regionally correct perturbation. Using the response functions obtained from Eqs. (27) and (28), one can now write Eq. (14) for global, tropical and extra-tropical carbon. Plugging the result into Eq. (24) gives

Hence, one can infer that

Therefore, the global response function χ_β(t) can be obtained from ${\mathit{\chi }}_{\mathit{\beta }}^{\text{tr}}\left(t\right)$ and ${\mathit{\chi }}_{\mathit{\beta }}^{\text{et}}\left(t\right)$. However, in addition, since χ(t) is given by Eq. (5), Eq. (30) implies that one can also obtain the global spectrum q_β(τ) by combining the regional spectra:

Therefore, using the recovered response functions for tropical and extra-tropical carbon one can obtain the global response function χ_β and its associated spectrum q_β. Accordingly, in this test we check Eqs. (30) and (31). For the calculation of the derivatives in Eqs. (27) and (28) we fitted NPP^tr=NPP^tr(C_atm) and NPP^et=NPP^et(C_atm) once again by a polynomial of order 6 (which obtained the best results for global NPP in Fig. 7b) and took the derivatives from the fits. For χ_β(t) we used the recovery with best quality from Fig. 7 (“NPP${}_{\text{pol}\mathrm{6}}^{\mathrm{2}\times }$”). The spectra q_β(τ), $q_{\mathit{\beta }}^{\text{tr}}\left(\mathit{\tau }\right)$ and $q_{\mathit{\beta }}^{\text{et}}\left(\mathit{\tau }\right)$ are obtained by scaling the spectra from χ_NPP(t), ${\mathit{\chi }}_{\text{NPP}}^{\text{tr}}\left(t\right)$ and ${\mathit{\chi }}_{\text{NPP}}^{\text{et}}\left(t\right)$ by the respective derivatives $\frac{\partial \text{NPP}}{\partial C_{\text{atm}}}{|}_{C_{\text{atm}}=C_{\text{atm}}^{\mathrm{0}}}$, $\frac{\partial {\text{NPP}}^{\text{tr}}}{\partial C_{\text{atm}}}{|}_{C_{\text{atm}}=C_{\text{atm}}^{\mathrm{0}}}$ and $\frac{\partial {\text{NPP}}^{\text{et}}}{\partial C_{\text{atm}}}{|}_{C_{\text{atm}}=C_{\text{atm}}^{\mathrm{0}}}.$ Results are shown for χ_β in Fig. 10b and for q_β in Fig. 10c. As seen in Fig. 10b, the reconstruction matches almost perfectly the values of χ_β for times beyond about 20 years. Likewise, the spectrum q_β is almost perfectly reconstructed for timescales above 6 years, with a slight discrepancy between 15 and 25 years (see Fig. 10c). A larger disagreement is seen below 20 years for χ_β and below 6 years for q_β. One of the reasons is that only little information is available for timescales smaller than 1 year because this is the time step taken for the data. However, Appendix D shows that timescales much shorter than the time step affect only χ(0). The main reason for the disagreement is that high frequencies (and thus small timescales) are the most problematic to recover due to the ill-posedness of the problem that obscures the signal particularly at high frequencies.

Despite the disagreement at small timescales, Fig. 10b and c add confidence that (1) the response functions for global, tropical and extra-tropical carbon can be trusted over mid to long timescales (Fig. 10b) and that (2) the two-peak spectrum obtained for global land carbon indeed characterizes the response, since its computation from two independent approaches (combining regional spectra in Fig. 10c and the Gregory-type plot in Fig. 9c) yield similar results with characteristic timescales matching the peaks of the spectrum.

5.4 Investigating the two-peak structure of the spectrum

The reasons for the two-peak structure of the land carbon spectrum are conceivable. In principle, one possibility could be that the short timescales reflect the carbon dynamics in the tropics, a region with higher temperatures and thus larger heterotrophic respiration rates (see, e.g., Raich and Potter, 1995), while the long timescales may originate from the carbon dynamics in the extra-tropics, where respiration rates are smaller due to lower temperatures. This hypothesis may be checked by examining Fig. 10d, which shows the contribution from the tropics and the extra-tropics to the land carbon spectrum (Eq. 31). However, as seen, the two peaks arise both in the tropical and in the extra-tropical spectrum, so that one cannot attribute each peak to a particular region. Therefore, this cannot be the explanation.

An alternative hypothesis is that the different peak timescales originate from the very different characteristic timescales of functionally different elements in the land carbon cycle such as leaves, wood, litter and soil. In the following we investigate whether this hypothesis can explain the two-peak structure.

For this purpose, we split the land carbon pools of the MPI-ESM into two groups (see Fig. 11). In the first group are the pools whose dynamics is governed by fast processes such as shedding and decomposition of leaves, thus the pools that respond at short timescales. These are the pools representing non-woody tissues in living vegetation (leaves, fine roots, sugars, and starches) and the associated litter. In the second group are the pools with dynamics determined by slow processes such as the decomposition of woody plant parts, hence the pools that respond at long timescales. These are the pools representing the wood in stems and coarse roots, the associated litter pools, and the soil carbon (humus).

Now, we separate the land carbon spectrum into contributions arising from each of these two groups. The land carbon response can be described as the sum of the collective responses from the pools with fast and those with slow dynamics:

In the linear response framework, the response of each term to NPP is given by

which implies

Finally, assuming that each response function in Eq. (34) is described by Eq. (5), we obtain the separation of the land carbon spectrum in terms of the contribution from each pool group:

Hence, if our hypothesis is correct, then the peak at short timescales should be explained by $q_{\text{NPP}}^{\mathrm{f}}$ and the peak at long timescales should be explained by $q_{\text{NPP}}^{\mathrm{s}}$.

Figure 11Simplified pool scheme of the land carbon cycle in the MPI-ESM (see Appendix A) with pools split into two groups according to the characteristic timescales underlying their dynamics.

Download

To proceed with the analysis we now need to obtain the spectra $q_{\text{NPP}}^{\mathrm{f}}$ and $q_{\text{NPP}}^{\mathrm{s}}$. When investigating the tropical and extra-tropical land carbon, we obtained each regional spectrum individually by applying the RFI algorithm to the data from the tropical and extra-tropical responses. In principle one could proceed in the same way to separate contributions from the two pool groups, but, as will become clear below, this approach introduces slight inconsistencies between the separate recoveries that makes a quantitative comparison of the two contributions less reliable than the alternative method that we use in the following to separate the fast and slow components of q_NPP.

The idea of this alternative approach is the following. Numerically the land carbon spectrum is given by the regularized solution (Eq. 6), i.e.,

where the regularization parameter λ is determined by the RFI algorithm. The land carbon response ΔC^bgc is given by the sum (32) of the responses from the pools with fast and those with slow dynamics. Entering Eq. (32) into Eq. (36) yields

Therefore, by deriving the spectra ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{f}}$ and ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{s}}$ using the same regularization parameter λ employed for the land carbon spectrum q_NPP,λ, we can in principle accurately reconstruct q_NPP,λ from ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{f}}$ and ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{s}}$.

By obtaining ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{f}}$ and ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{s}}$ in this way and combining them via Eq. (35) we show that q_NPP,λ can indeed be very accurately reconstructed (Fig. 12a). This approach also leads to an almost perfect reconstruction of the response function χ_NPP via Eq. (34), as shown in Fig. 12b. Compared to our previous approach employed to reconstruct the global spectrum from the regional spectra (Eq. 31), this alternative approach gives a more accurate reconstruction because it takes the same regularization parameter λ for all q_NPP,λ, ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{f}}$, and ${\mathit{q}}_{\text{NPP},\mathit{\lambda }}^{\mathrm{s}}$, in contrast to the previous approach where λ was separately calculated by the RFI algorithm for each spectrum in Eq. (31).

Figure 12Investigation of the contribution from the pools with fast dynamics (non-woody pools) and those with slow dynamics (woody and humus pools) to the land carbon spectrum. The analysis is based on the 2×CO₂ bgc experiment. (a) Land carbon spectrum q_NPP and its reconstruction via Eq. (35) from the fast and slow components; (b) response function χ_NPP and its reconstruction via Eq. (34) from the response functions for the fast and slow components; (c) separation of spectrum q_NPP into $q_{\text{NPP}}^{\mathrm{f}}$ and $q_{\text{NPP}}^{\mathrm{s}}$; (d) separation of response function χ_NPP into ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{f}}$ and ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{s}}$; (e) separation of χ_NPP into contributions from the individual pools. For more details, see the text.

Download

Now, the question is whether the spectra $q_{\text{NPP}}^{\mathrm{f}}$ and $q_{\text{NPP}}^{\mathrm{s}}$ can indeed explain each peak in the land carbon spectrum q_NPP. To check this, in Fig. 12c we plot the three spectra. As seen, clearly the peak at short timescales of the land carbon spectrum arises mostly from the large peak in the fast-dynamics spectrum, while the peak at long timescales follows closely the large peak in the slow-dynamics spectrum. This result thus indicates that our hypothesis is correct, so that the two-peak spectrum indeed originates from the different characteristic timescales of functionally different elements in the land carbon cycle.

More insight into the dynamics of the land carbon can be gained by analyzing the different contributions to the response function χ_NPP shown in Fig. 12d and e. Figure 12d shows the contribution of the pools with fast dynamics ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{f}}$ and that of the pools with slow dynamics ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{s}}$. We see that ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{f}}$ only dominates the response at times smaller than about 2 years. Further, for times larger than 25 years only ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{s}}$ contributes to the response. In Fig. 12e we further separate the response functions into the contributions from the individual pools, i.e.,

where ${\mathit{\chi }}_{\text{NPP}}^{\text{nwt}}$ is the response function for the non-woody tissues in living vegetation, ${\mathit{\chi }}_{\text{NPP}}^{\text{nwl}}$ for non-woody litter, ${\mathit{\chi }}_{\text{NPP}}^{\text{wt}}$ for woody tissues in living vegetation, ${\mathit{\chi }}_{\text{NPP}}^{\text{wl}}$ for woody litter, and ${\mathit{\chi }}_{\text{NPP}}^{\mathrm{h}}$ for the humus pool. This more detailed separation shows that for long times the humus pool starts to dominate the response, so that at times larger than 100 years it gives the only contribution.