Mathematical Reconstruction of Land Carbon Models From Their Numerical Output: Computing Soil Radiocarbon From 
 12 C Dynamics

Radiocarbon (14C) is a powerful tracer of the global carbon cycle that is commonly used to assess carbon cycling rates in various Earth system reservoirs and as a benchmark to assess model performance. Therefore, it has been recommended that Earth System Models (ESMs) participating in the Coupled Model Intercomparison Project Phase 6 report predicted radiocarbon values for relevant carbon pools. However, a detailed representation of radiocarbon dynamics may be an impractical burden on model developers. Here, we present an alternative approach to compute radiocarbon values from the numerical output of an ESM that does not explicitly represent these dynamics. The approach requires computed 12C stocks and fluxes among all carbon pools for a particular simulation of the model. From this output, a time‐dependent linear compartmental system is computed with its respective state‐transition matrix. Using transient atmospheric 14C values as inputs, the state‐transition matrix is then applied to compute radiocarbon values for each pool, the average value for the entire system, and component fluxes. We demonstrate the approach with ELMv1‐ECA, the land component of an ESM model that explicitly represents 12C, and 14C in 7 soil pools and 10 vertical layers. Results from our proposed method are highly accurate (relative error <0.01%) compared with the ELMv1‐ECA 12C and 14C predictions, demonstrating the potential to use this approach in CMIP6 and other model simulations that do not explicitly represent 14C.


Introduction
Terrestrial ecosystems, and soils in particular, are important global carbon (C) reservoirs but with large uncertainties in the size of the stock, process rates, and their representation in Earth System Models (ESMs) (Luo et al., 2016;Sulman et al., 2018). The dynamics of radiocarbon in various components of terrestrial ecosystems can be used to interpret the rates at which CO 2 is removed from the atmosphere into plants through photosynthesis, moves from vegetation into soils through litter fall and mortality, and finally cycles back to the atmosphere via respiration. As a result, radiocarbon has been used in many empirical studies to evaluate rates of soil C cycling and, more recently, as a tool for model evaluation (Ahrens et al., 2015;Chen et al., 2019;Dwivedi et al., 2017) and intercomparison (He et al., 2016).
Measurements of radiocarbon offer great potential to constrain temporal dynamics of the land carbon cycle in ESMs as well as to assess their performance. However, to date, such comparisons have been limited to a small number of models that simulate C isotopes (e.g., Chen et al., 2019;Koven et al., 2013;Tifafi et al., 2018), as the majority of models do not explicitly represent isotope dynamics. The current coupled climate-carbon cycle model intercomparison project (C4MIP) of the Coupled Model Inter-comparison Project Phase 6 (CMIP6) requests model participants to explicitly model radiocarbon and include the corresponding output from the major land and ocean model compartments (Graven et al., 2017;Jones et al., 2016). However, this additional request can be time-consuming to implement and computationally intensive to run. Thus, the goal of this work is to develop a tool to make radiocarbon-based comparisons of C cycling practical and comparable across different ESMs even if they do not explicitly include these dynamics.
The availability of radiocarbon output from a range of models would allow an additional comparison and evaluation against global data sets, such as the International Soil Radiocarbon Database (ISRaD; Lawrence et al., 2019), and would provide a new way to differentiate patterns of C cycling across different model structures (Sierra et al., 2014). There are literally hundreds of models of soil carbon cycling described in the literature (Manzoni & Porporato, 2009), each relying on different hypotheses about internal system dynamics, with specific model structures and parameterizations. Generalizations of models have also been proposed based on basic assumptions of soil organic matter dynamics (Sierra & Müller, 2015). Computing radiocarbon output from these different classes of models, without the need to incorporate explicit isotopes dynamics in the model, would benefit data-model comparisons and improve the understanding of soil C dynamics in these models.
Linear response theory suggests that the impulse response function of a particular model can be used to compute isotope dynamics without explicitly simulating the isotope within the model but rather using the response function to convolve the isotope tracers (Thompson & Randerson, 1999). Therefore, it may not be necessary to explicitly represent isotope dynamics a priori in a model; rather, they can be obtained a posteriori through the convolution procedure. However, linear response theory relies on assumptions of steady-state and model linearity, which are not met for the type of transient simulations performed with ESMs. Metzler et al. (2018) recently developed a method to obtain a mathematical object, the state-transition matrix (see details below), which generalizes response function theory for nonlinear models out of equilibrium. This approach removes previous assumptions that prevented the use of the linear response approach for the computation of isotope and tracer dynamics a posteriori from simulations. Therefore, the state-transition matrix approach offers an opportunity to compute isotope dynamics from models that do not explicitly represent them and could be applied to models participating in CMIP6, for example. However, since output from these models is provided in discrete time steps often much larger than the internal model time step, it is a challenge to translate the state-transition-matrix approach from the continuous to the discrete case.
We introduce in this manuscript an approach to mathematically reconstruct land carbon models from simulation output and use the reconstruction to compute radiocarbon values. To this end, we extend the method presented in Metzler et al. (2018) for the reconstruction of a differential equation model from an existing solution and the computation of the state-transition matrix in continuous time. We show how to translate the approach to the discrete-time case and apply it to reconstruct tracer and isotope dynamics from model output. We first give a mathematical overview of the approach and then present examples of its use. We test our approach by comparing our computed values to output from the ELMv1-ECA model (Chen et al., 2019;Riley et al., 2018;Zhu et al., 2019), the land component of the Exascale Earth System Model (E3SM) model, which explicitly represents radiocarbon and is therefore an appropriate test case to evaluate our proposed method. We evaluate the modeled radiocarbon values along the entire vertically resolved soil profile and across all represented soil and litter pools over time. We believe this new framework has potential for application to the CMIP6 simulations and to models with more diverse structures.

Carbon Cycle Models as Compartmental Systems
Carbon cycle models are subject to the law of mass conservation and therefore can be represented as compartmental systems (Anderson, 1983;Jacquez & Simon, 1993). A compartmental system with d ≥ 1 compartments is of the general form of the initial value problem d dt Here, C(t) is the vector of compartment contents at time t, B(C, t) is a matrix-valued function that governs the cycling of carbon through the system, u(C, t) is the vector of external carbon inputs to the system, t 0 is some fixed initial time, and C 0 is a given vector of initial system contents. Note that the internal cycling as well as the external inputs can each depend on time (e.g., due to changing environmental conditions) and on the current carbon content of the system itself. Consequently, system (1) is nonautonomous (time dependent) and nonlinear.
To guarantee mass conservation, the matrix B(C, t) is required to be compartmental for all nonnegative C and t > t 0 ; that is, B(C, t) = (B i (C, t)) i, =1,2, … ,d satisfies the properties The diagonal values of B ii can be interpreted as the cycling rate of compartment i, the off-diagonal entries B i represent the rate of carbon flux from compartment to compartment i, and z = − ∑ d i=1 B i is the rate of carbon leaving the system through pool .
We now assume that we know a unique (numerical) solution C = C(t) of system (1). Plugging this solution into the system, we obtain a new compartmental system: whereB(t) = B(C(t), t) andũ(t) = u(C(t), t). This system is still time dependent, but linear. Under the assumption that it has a unique solutionC =C(t), the two solution trajectoriesC and C are identical. Then the linear compartmental system (2) describes exactly the same solution trajectory as the nonlinear compartmental system (1). For the sake of simplicity of notation, we omit the tilde in equation (2) in the remainder of the manuscript.
The unique semianalytical solution of the linear compartmental system (2) is given by (Brockett, 2015) We call the solution semianalytic because the system's state transition matrix Φ is only implicitly given as the unique solution of the matrix initial value problem: where Id denotes the identity matrix in R d . The state transition matrix contains information about the entire dynamics of the original system along the solution trajectory. It can be used to reconstruct the original nonlinear model with the specific forcings specified for the particular simulation.

Representation of Depth-Resolved Systems
This framework can also incorporate spatially explicit systems that move carbon along vertical or horizontal gradients by, for example, advection and diffusion. For instance, consider a depth-resolved compartmental model with l depth layers and c compartment types per layer. Let d ∶= c · l. Values of d from 1, 2, … , l belong to Layer 1, Layer 2, … , and Layer l of compartment c = 1. The numbers l + 1, l + 2, … , 2 · l belong to Layer 1 through Layer l of compartment c = 2 and so forth. In general, the numbers (c − 1) · l + 1, (c − 1) · l + 2, … , c · l belong to Layer 1 through Layer l of compartment type c. We introduce B ∶= (B i ) i, =1,2, … ,d and u ∶= (u i ) i=1,2, … ,d . Then, the depth-resolved system fits in the framework of equation (1), and the nonlinear structure of the system accommodates interactions among system variables that are commonly used to represent movement along spatial gradients. For a system with 10 depth layers and 7 compartment types per layer, B is a 70 × 70 matrix, and u is a vector with 70 components. For example, B 12,11 is the flux downward from layer l = 1 to layer l = 2 of compartment c = 1, and B 35,55 is the exchange rate from compartment c = 5 to compartment c = 3 within depth layer l = 5 (recall that the first index represents the target and the second index the source). Obviously, fluxes among different compartment types and vertical layers can be represented, for example, by B 23,52 .
The negative diagonal entries (B ii ) i=1,2, … ,70 are the total exit rates from each compartment, that is, the rates to other compartments in any layer plus the rate of material advecting or diffusing from the compartment.

The Radiocarbon Compartmental System
Once we know the representation (1) of a carbon cycle model, we can readily manipulate the equation to obtain a compartmental system for radiocarbon. We only need to substitute B by 14 B defined as where is the radioactive decay constant of 14 C scaled to the relevant time unit. For example, if the system is observed in time units of years, then = ln(2)/5,568 years. This manipulation leads to a state-transition matrix of the radiocarbon system given by from which we can directly calculate the effect of radioactive decay.

Reconstruction of a Continuous-Time System From Discrete-Time Model Output
The approach presented above works only when the original model is completely described by a system of ODEs, which by definition implies that the system is continuous in time. However, we are interested in using the approach for ESMs for which we only have access to model output reported in discrete time steps and for which the underlying process descriptions may be a combination of ODEs, PDEs, noncontinuous representations, and so forth. Therefore, it is necessary to produce a Continuous-Time Approximation (CTA) of the model from its discrete model output. First, we will illustrate the approach using a one compartment system as an example and then for multiple interconnected compartments.

One compartment Case
Suppose we are given a one compartment time-dependent linear compartmental system d dt where (t) and u(t) are unknown. Further, for discrete times t 0 < t 1 < t 2 < … < t n , we are given • the initial system state C 0 , Here, u k and r k denote the accumulated external system input and the accumulated external system output in the time interval We are interested in constructing a compartmental system, as simple as possible, that matches these given data as well as possible. To that end, we approximate the (unknown) system (3) on each interval I k by a linear time-independent system Suppose that at time step k the valueĈ k is already known from the previous time step, and for k = 0 we set C 0 = C 0 .
We write Δ k = t k+1 − t k and choose û k = u k ∕Δ k . Consequently, for the time interval I k the cumulative external system input u k of the unknown system (3) and the cumulative external system input Δ k û k of the approximating system (4) coincide. We are left with finding good approximationŝk for (t) on

10.1029/2019MS001776
I k , k = 0, 1, … , n − 1. We can compute the cumulative external system output of the approximating system (4) on I k , in dependence on the choice of̂k, bŷ is the unique solution of (4) on I k . Consequently, we obtain̂k by solving the one-dimensional nonlinear optimization problem We defineĈ k+1 =Ĉ(t k+1 ) and continue with the next time step k + 1.
In the special case of a time-independent unknown system (3), that is, (t) = < 0 and u(t) = u ≥ 0, the approximation by system (4) is perfect, independent of the number n of given data points.

Multiple Compartment Case
If the compartmental system consists of more than one compartment, internal cycling among the different compartments may occur. In particular, all modern soil carbon models are multicompartmental (Manzoni & Porporato, 2009;Sierra & Müller, 2015), and they often represent vertical structure with explicit representation of diffusive or advective fluxes that connect adjacent soil layers in the discretized model. An approximating system is then required to also match the internal fluxes. Suppose we are given a d-compartment time-dependent linear system d dt Here, u k and r k denote the accumulated external system input vector and the accumulated external system output vector in the time interval I k = [t k , t k+1 ], k = 0, 1, … , n − 1, respectively. Furthermore, F k i denotes the accumulated flux from compartment to compartment i during I k .
We are interested in constructing a compartmental system, as simple as possible, that matches these given data as well as possible. To that end, we approximate the (unknown) system (5) on each interval I k by a linear time-independent system Suppose that at time step k the vectorĈ k is already known from the previous time step; for k = 0 we setĈ 0 = C 0 .
Again, we write Δ k = t k+1 − t k and chooseû k = u k ∕Δ k . Consequently, for the time interval I k the cumulative external system input vector u k of the unknown system (5) and the cumulative external system input vector Δ kû k of the approximating system (6) coincide. We are left with finding good approximationsB k for B(t) on I k , k = 0, 1, … , n − 1 under the constraint thatB k is a compartmental matrix. To that end, we try to match the internal fluxes F k i and the external output flux vectors r k as well as possible. For i ≠ , we can compute the internal flux from to i of the approximating system (6), in dependence on the choice ofB k , during the is the unique solution of (6) on I k and e (t−t k )B k denotes the matrix exponential.
i ≥ 0 the rate of external outflow from compartment . Then we can compute the cumulative external system output through compartment of the approximating system (6) during I k , in dependence on the choice ofB k , bŷ Consequently, we obtainB k by solving the m-dimensional nonlinear optimization problem where the maximum difference between predicted and provided fluxes is minimal as well as the difference between predicted and provided release fluxes. The dimension m of the optimization problem is at maximum d 2 . If some internal fluxes or external output fluxes are inexistent, then m < d 2 . We defineĈ k+1 =Ĉ(t k+1 ) and continue with the next time step k + 1.

Reconstruction of a Discrete-Time Compartmental System
Model output is usually provided at discrete time points t 0 < t 1 < … < t n , either because there is no feasible method to represent continuous-time model output or the underlying process model works on the basis of discrete time steps. In either case, if we reconstruct the unknown original model using a CTA (6), we interpolate the given data for points in time that lie between t k and t k+1 . This interpolation can be done in many ways, some of which could lead to large errors.
To ameliorate this problem, we can restrict the model reconstruction only to the data provided by the model at those particular time steps. Consequently, we reconstruct the unknown model in a discrete time setting, on the time grid t 0 < t 1 < … < t n that is predetermined by the given data and make no interpolations for other time points. This approach leads to a reconstructed model that is represented by a d-compartment discrete-time system of the form The discrete-time compartmental matricesB(k) satisfy for all k = 0, 1, … , n − 1 two conditions:

10.1029/2019MS001776
Analogously to the CTA,ẑ (k) = 1 − ∑ d i=1Bi (k) is the rate at which mass at time step k leaves the system through compartment . The discrete-time state-transition matrix is, for 0 ≤ k 1 < k 2 ≤ n, given bŷ andΦ(k, k) = Id. The solution of the discrete-time initial value problem (7) is then given bŷ Suppose we are given a d-compartment time-dependent linear system with structure (5) and unknown B and u. For discrete times t 0 < t 1 < t 2 < … < t n we are given • the initial system content vector Here, u k and r k denote the accumulated external system input vector and the accumulated external system output vector in the time interval I k = [t k , t k+1 ], k = 0, 1, … , n − 1, respectively. Furthermore, F k i denotes the accumulated flux from compartment to compartment i during I k .
For cases where the full model output temporal resolution is too coarse, which can lead to depletion of a pool and negative pool sizes during DTA reconstruction, the Continuous-Time Approximation (CTA) is required.

Application of the Framework to Output From the ELMv1-ECA Model
As the land component of E3SM, ELMv1-ECA Zhu et al., 2019) simulates vertically resolved soil C, N, and P dynamics based on the Century soil model (Parton et al., 1993) with enhancements to include soil O 2 effects (Koven et al., 2013). Radiocarbon ( 14 C) is treated as a separate tracer in addition to 12 C. 14 C enters into a terrestrial ecosystem from the atmosphere (Graven et al., 2017), and it is stored in vegetation compartments through biomass growth. After vegetation turnover (e.g., natural phenology, mortality, and disturbances), 14 C flows into litter or coarse woody debris along with 12 C and is consequently decomposed or accumulated in soil organic matter. Besides microbial decomposition, 14 C also radioactively decays. ELMv1-ECA also added several new features that could affect C inputs to soil and SOC dynamics that potentially have great impacts on soil radiocarbon signals. Briefly, the changes relevant to the current study include (a) phosphorus limitation on soil organic carbon decomposition (Zhu et al., 2016); (b) application of the Equilibrium Chemistry Approximation of multiple-consumer-multiple-resource competition network to resolve microbial immobilization of soil available nutrients  10.1029/2019MS001776 Zhu et al., 2017); and (c) dynamic plant stoichiometry that affects litter nutrient content and decomposability. ELMv1-ECA has been evaluated against multiple global-scale observations of ecosystem C, water, and energy stocks and fluxes using ILAMB (Collier et al., 2018;Zhu et al., 2019) and short-term nutrient cycling observations and global-scale partitioning of N losses Zhu & Riley, 2015). Overall, model benchmarking shows improvement from the precursor CLM4.5 predictions, and in particular for this study, more accurate estimates of spatially distributed soil C stocks .
In this study, we present simulation results from ESMv1-ECA at three sites with contrasting background climate dynamics: boreal ( . These sites were selected due to their representative background climates and plant functional types representing high-, middle-, and low-latitude ecosystems. Characteristics for the three sites were extracted from the standard global model datasets, which are used to set plant functional types, soil properties, and so forth. . The simulation protocol for all three sites follows the standard approach described in Oleson et al. (2010). The model is first spun up (i.e., brought close to equilibrium) using an accelerated soil decomposition approach (Koven et al., 2013) for 1,000 years, followed by a 400-year regular spinup with regular soil decomposition. Soil phosphorus pools were initialized from an observationally derived data set (Yang et al., 2013) at the beginning of the regular spinup. The spinup simulations were forced with repeated meteorology (GSWP reanalysis; Dirmeyer et al., 2006;from years 1901to 1920 and constant atmospheric CO 2 mole fraction (285 ppm). We evaluated the model's spun-up soil carbon stock to ensure changes at the end of the spinup are less than 0.01%/year. After spinup, the model was run in a transient simulation from 1850 to 2010 with GSWP reanalysis forcing (Dirmeyer et al., 2006), transient CO 2 concentrations, N deposition (Lamarque et al., 2005), and P deposition (Mahowald et al., 2008). Atmospheric Δ 14 C is prescribed to be zero during spinup simulations and using observed values during transient simulations (Levin et al., 2010). The spinup procedure ultimately determines the distribution of radiocarbon in the vector of initial carbon contents. Other spinup procedures may lead to different initial radiocarbon distributions; however, this is an issue particular for each simulation setup without relevance for the reconstruction of the model in the transient simulation.
We performed analyses using different approaches to reconstruct ELMv1-ECA predictions from model output. To characterize the effect of temporal discretization, we applied the DTA approach using ELMv1-ECA output averaged over 1-and 10-day time steps. This analysis will inform required output frequency for other models to apply model reconstruction.

Results and Discussion
We first give two examples with simplified systems (i.e., one-and two-compartment systems) to demonstrate characteristics of the Continuous-Time Approximation (CTA) and then present comparisons using Discrete-Time Approximation (DTA) to the full ELMv1-ECA predictions.

One compartment Example
In some cases, it is useful to represent soil carbon as one single compartment with an aggregated decomposition rate (Olson, 1963). This assumption may have important limitations to represent complex dynamics in soils (Manzoni & Porporato, 2009;Sierra & Müller, 2015), but it can be very useful to introduce complex concepts using a simplified system. We consider first a one compartment system with constant inputs and rates that represents the accumulation of carbon over time starting with a known initial amount of carbon d dt C(t) = −0.04 C(t) + 3, t > 1909, The CTA approach exactly matches the system's analytical solution despite the relatively small number of data points available for the reconstruction, n = 3, 6 ( Figure 1). However, the CTA can give approximations with large errors in the time-dependent case. Consider the one compartment system with time-dependent inputs beginning in year 1909: The CTA improves with increasing number n of data points used in the approximation (Figure 2). For small number of data points, the CTA can miss important dynamics occurring inside the time step. For instance, for the reconstruction using time steps of 30 years, the convexity of the curve is completely different between the original model and the reconstruction. Therefore, care must be taken in using the CTA when important changes in dynamics occur inside the time step.

Multicompartment Test
As a preliminary test before demonstrating the approach for the full ELMv1-ECA model reconstruction, we evaluated the approach for several configurations of multiple C pools and nonlinear time dependencies. Figure 2. Transient C stocks of the one-compartment system with time-dependent inputs reconstructed by CTA for different data time steps Δt. Mean annual absolute bias is 3.92%, 0.96%, and 0.43% for Δt = 30, 20, and 10 years, respectively. The solid black line represents the carbon content of the original system (9); the dashed lines represent the carbon contents of the CTAs with different time steps.
We found that the accuracy of the CTA approach does not depend on the number of pools or the nonlinear time dependencies. The time step required for accurate reconstruction depends on the frequency of the forcing and nonlinear interactions among pools. For example, consider a two compartment system, which can represent the dynamics of a fast and a slow soil carbon pool (Andren & Katterer, 1997), initialized at year t = 1909 and propagated for 110 years: , where 1 (t) = 0.1 + 0.05 sin(0.1 t), 2 = 0.05.
In this case, data from the original model extracted at 5-year time steps leads to a very accurate CTA, whereas coarser temporal time steps lead to progressively worse approximations (Figure 3).

ELMv1-ECA Approximation
As long as the time-step is short enough, we can use either CTA or DTA to reconstruct the ELMv1-ECA model. CTA is generally more demanding in terms of computations because it needs to solve a system of ODEs for each time step. For practical reasons, we present here results using the DTA approach as described previously, compute radiocarbon values for all pools, and compare with the computations of radiocarbon performed in the original model. To evaluate errors related to the size of the time step, we performed simulations using 1-and 10-day ELMv1-ECA model outputs of 12 C inputs, fluxes between soil pools, and atmospheric Δ 14 C.
Our DTA accurately reproduced the ELMv1-ECA bulk soil 12 C and Δ 14 C values across the three sites and at all soil depths (Figure 4). The bulk soil Δ 14 C values shown in Figure 4 combined seven SOM component pools (coarse woody debris [CWD]; metabolic, cellulose, and lignin litter pools [Litters 1, 2, and 3, respectively]; and fast, intermediate, and passive soil pools [Soils 1, 2, and 3, respectively]), which were separately reproduced with the DTA. The mean absolute biases were 0.0003%, 0.0014%, and 0.0003% for soil, litter, and coarse woody debris, respectively.
The DTA also accurately reproduced temporal variability of ELMv1-ECA modeled Δ 14 C values for all these pools over the 110-year simulation from 1901 to 2010 ( Figure 5). For all sites, the mean absolute error was <0.001% for soil carbon stocks and Δ 14 C values at all depths (Table 1). These extremely low reconstruction errors are well within model uncertainty and far less than within-gridcell observed spatial heterogeneity of soil carbon stocks and radiocarbon values (Chen et al., 2019;Lawrence et al., 2019) therefore, tolerable within the context of soil C modeling.
During the 1960s, the radiocarbon concentration in the atmosphere increased sharply due to nuclear weapons testing (generally known as the bomb spike). The CWD, litter, and fast soil pools responded relatively rapidly to this change, as atmospheric 14 CO 2 rapidly cycled through these pools. This bomb carbon 14 C signal was also evident in the top meter of soil ( Figure 4). In contrast, the deeper and slower soil pools (i.e., intermediate and passive) had a more lagged and damped response to the atmospheric bomb-radiocarbon perturbation.
Although the relative errors were extremely low across all pools over time, the percent error was larger for the CWD and litter pools than for the soil pools, with the error in the litter pool reaching a maximum during the bomb-spike period at around 0.2% ( Figure 6). The mean absolute percent errors for DTA 14 C values were, in most cases, at least 2 orders of magnitude larger than those for 12 C across sites and pools (Table 1). The DTA 14 C error increased for larger time steps, being an order of magnitude larger for CWD and litter pools in the 10-day versus 1-day time step cases. However, the mean absolute percent error for the 10-day time step case remained very low: < 1 × 10 −4 % for 12 C and < 0.01 % for 14 C for all pools across different sites (Table 1), remaining well within model uncertainty bounds and soil heterogeneity.

Performance and Future Application
Using a 1-day time step, the DTA accurately matched ELMv1-ECA modeled Δ 14 C values (<0.001% relative error for most sites and C pools). Although the error rate increases for larger time steps, acceptable low error (<0.01% relative error) is still achieved with a 10-day time step (Table 1). These results suggest that the error would likely remain within reasonable bounds for a 30-day time step, which is common output from most ESMs. We demonstrated here that the CTA and DTA approaches can be applied to numerical output from land carbon models to reconstruct a linear version of the model that very closely matches the output of the original model. The approximations can then be used to simulate isotope (e.g., 13 C or 14 C) dynamics consistent with the original model formulation of 12 C dynamics. Furthermore, the approximations can be used for both single-depth and vertically resolved soil radiocarbon models run at different spatial and temporal resolutions. The approach is not only restricted to soils but can also be used to reconstruct whole ecosystem carbon and isotope dynamics. We recommend using the DTA for reconstructing ESM-scale models, since the CTA makes implicit assumptions about what the original model does between time steps, and these assumptions may be unjustified in some cases depending on time-step size.
This DTA approach has widespread applicability, as it can be used across a range of models independent of their underlying structure. The approach does not require the model to implement additional isotope dynamics or tracers, but simply to output 12 C stocks and fluxes. For maximum accuracy, all model stocks and fluxes are required (e.g., all individual transfers between all pools represented in the model). Aggregating C stocks and fluxes results in greatly increased uncertainties. However, given the uncertainties inherent in ESM land model soil carbon dynamics, reconstruction error may be lower than uncertainty in model structure and parameterization (Chen et al., 2019). Error arising from pool aggregation is of particular importance when models report carbon stocks only for aggregated pools (e.g., sum of all soil carbon pools in a particular soil layer). In such cases, the original and reconstructed model will differ in structure, leading to error in the reconstruction. This error would probably depend on the number and size of pools being aggregated, and on the temporal resolution of the output. Numerical errors and instabilities may also arise from empty or very small pools. For instance, we observed larger relative errors in the relatively small deep litter and CWD pools.
The potential of this approach to produce radiocarbon predictions across models that do not explicitly represent it will allow us valuable data-model comparisons and enable use of radiocarbon as a constraint for parameter estimation (He et al., 2016) Although radiocarbon is commonly associated with C age, it cannot be directly interpreted as such in soils, due to mixing of C inputs of different ages in the open system . Instead, radiocarbon measurements and model output are primarily of value as a tracer, and for direct model-data comparison. In addition, by using a state-transition matrix approach as we did here, there is the added possibility to calculate age and transit time distributions for the system (Metzler et al., 2018). In the future, this approach could be applied to compute radiocarbon values, ages, and transit times across a suite of models and evaluate Figure 6. Relative error in radiocarbon values for CWD, aggregated litter, and soil pools over time for the temperate forest site. 14 C relative error is approximately 2 orders of magnitude larger than 12 C error but is always small in comparison to model uncertainty and spatial heterogeneity ( the time carbon spends in terrestrial ecosystems. This approach would provide meaningful comparisons of C ages and transit times between model structures and transient scenarios. In addition, the reconstruction we propose of numerical model output as a dynamical system can be used to compute other system-level diagnostics of the terrestrial carbon cycle. For instance, Luo et al. (2017) proposed a set of diagnostics to evaluate carbon storage capacity and potential from transient simulations and developed a framework to assess model performance based on the ability to trace particular components of ecosystems (Luo et al., 2015;Xia et al., 2013). From the numerical model reconstructions we proposed here, it is possible to perform these analyses without the need to analytically rewrite existing models in matrix form. This means that it is still possible to obtain time-dependent versions of these matrices and vectors in numeric form, without the need to explicitly write by hand all equations of the model. This can be an advantage for complex models with multiple interactions between the carbon cycle and the water, nutrient, and energy cycles.

Conclusions
Numerical output from carbon cycle models can be used to reconstruct a linear version of the model that exactly matches the solution trajectory of the original model. This reconstruction can be used to compute other quantities not included in the original model such as isotope dynamics and other system-level metrics. The approach consists of computing a state transition matrix that contains all information of the dynamics of the system. Since available model output is always reported in discrete time steps, we proposed here two methods to approximate the original system, approximating a continuous system between time steps or approximating a discrete dynamical system. As long as the time step of the numerical output is sufficiently small to capture complex dynamics, both approaches provide very good approximations to the solution trajectory of the original system.
Radiocarbon is a valuable tracer to evaluate, understand, and improve soil carbon models. However, most land models do not currently include explicit representation of soil radiocarbon dynamics, precluding comparison with observations. Here we demonstrated how continuous and discrete time approaches can reconstruct soil radiocarbon profiles using only land model 12 C stocks and fluxes and atmospheric 14 C values. We applied the approach here to the ELMv1-ECA model Zhu et al., 2019), which explicitly represents soil organic matter 14 C dynamics with seven soil pools and ten vertical layers, and showed that our approach very accurately reproduces the full model simulations (bias < 0.01%). The approach can also be used for the new suite of CMIP6 simulations to evaluate whether the models are accurately representing soil C age and transit time, in addition to soil carbon stock (He et al., 2016;Lawrence et al., 2019). For applications to CMIP6 models, it is required that the models report all pools and all fluxes among pools as requested in the tier-2 output requirements of C4MIP (Jones et al., 2016).