Using generalized cross-validation to select parameters in inversions for regional carbon fluxes

[ 1 ] Estimating CO 2 fluxes from the pattern of atmospheric CO 2 concentrations with atmospheric transport models is an ill-posed inverse problem, whose solution is stabilized using prior information. Weights assigned to prior information and to CO 2 concentrations at different locations are quantified by parameters that are not well known, and differences in the choice of these parameters contribute to differences among published estimates of the regional partitioning of CO 2 fluxes. Following the TransCom 3 protocol to estimate CO 2 fluxes for 1992–1996, we find that the partitioning of the CO 2 sink between land and oceans and between North America and Eurasia depends on parameters that quantify the relative weight given to prior flux estimates and the extent to which CO 2 concentrations at different stations are differentially weighted. Parameter values that minimize an estimated prediction error can be chosen by generalized cross-validation (GCV). The GCV parameter values yield fluxes in northern regions similar to those obtained with the TransCom parameter values, but the GCV fluxes are smaller in the poorly constrained equatorial and southern regions.


Introduction
[2] To complement direct measurements of, for example, land or ocean carbon uptake, inverse modeling is widely used to estimate CO 2 fluxes from the observed spatial and temporal variations of atmospheric CO 2 concentrations. Different inverse modeling studies, however, have reached apparently contradictory conclusions on the longitudinal and land/ocean partitioning of CO 2 fluxes. For example, studies using similar data sets have variously placed a sink of order 1 Pg C yr À1 in temperate North America [Fan et al., 1998], in north Asia [Bousquet et al., 1999a], or in the north Atlantic and Pacific [Rayner et al., 1999].
[3] Estimating CO 2 fluxes involves the solution of a linear system where b is a vector of observed variations in CO 2 concentrations, E is a vector of random errors with zero mean and with covariance matrix cov(E) = C b , x is an unknown vector of CO 2 fluxes, and A is a transport operator that relates CO 2 fluxes to CO 2 concentrations.
[4] The solution of this inverse problem is generally not well constrained by the CO 2 concentrations (i.e., the problem is ill-posed) and must be stabilized through the use of prior information or regularity constraints [e.g., Hansen, 1998]. Estimates of CO 2 fluxes are usually constrained to be close to CO 2 fluxes specified a priori [e.g., Enting, 2002] by minimizing an object function of the form consisting of the sum of the least squares object function (first term) and a penalty term (second term) that penalizes deviations of the solution x from a given prior estimate x 0 . The covariance matrix C x represents uncertainty about the prior estimate x 0 , and the regularization parameter l sets the weight of the term involving the prior information relative to the least squares term. (See also the online supplement 1 .) [5] In CO 2 inversions, the covariance matrices C b and C x are usually taken to be diagonal, with diagonal entries c b and c x equal to assumed variances of the local CO 2 concentration errors and of the regional prior flux distributions. The prior fluxes x 0 and their assumed variances c x are typically chosen ad hoc from a range of reasonable values, as are the error variances c b of the CO 2 concentrations. It is known that inversion results can sensitively depend on the choice of such inversion parameters [Bousquet et al., 1999b;Rayner et al., 1999;Law et al., 2003], but a unified approach to quantifying this source of uncertainty and to choosing inversion parameters systematically has not been pursued.
[6] Several methods are available to choose inversion parameters systematically [e.g., Hansen, 1998, chapter 7]. One method applied in inverse problems in such fields as meteorological data assimilation [Wahba et al., 1995] and geodesy [Ditmar et al., 2003] is generalized cross-validation (GCV) [Golub et al., 1979]. A form of leave-one-out crossvalidation, GCV chooses parameters by minimizing an estimated mean-squared error of predictions with the model specified by the parameters. For the CO 2 problem, this means that GCV chooses inversion parameters by minimizing an estimated mean-squared error of predictions of CO 2 concentrations given estimated CO 2 fluxes and the model of atmospheric transport. Although such predictions are generally not the objective of inverse modeling of CO 2 fluxes, GCV provides a useful heuristic for choosing inversion parameters.
[7] We applied GCV to estimate two parameters in the TransCom 3 framework for inverse modeling of CO 2 fluxes [Gurney et al., 2002]. The two inversion parameters considered control the weighting of the prior flux estimates and the relative magnitudes of error variances of CO 2 concentrations at different locations. We examine how the CO 2 flux estimates depend on the values assigned to these parameters, choose parameter values by GCV, and discuss how the flux estimates thus obtained differ from those obtained with the original TransCom parameter values.

Inverse Modeling Framework
[8] The TransCom 3 annual-mean inversion [Gurney et al., 2002] uses measured CO 2 concentrations to estimate CO 2 fluxes from 11 land and 11 ocean regions after subtracting the estimated effects of fossil fuel emissions, of ocean carbon uptake, and of covariances between the seasonality of the biosphere and that of atmospheric transport (the seasonal ''rectifier effect''). The data used are the mean 1992 -1996 CO 2 concentrations at 76 stations and the mean growth rate of atmospheric CO 2 concentrations during that period, the latter constraining the sum of the regional fluxes. We followed the TransCom 3 protocol except for excluding one station, Darwin, whose mean concentration is anomalously high and apparently reflects local sources [Law et al., 2003], leaving 75 stations. (The TransCom 3 data and flux patterns are available at http://transcom. colostate.edu/TransCom_3/T3_Input_Data/. The prior fluxes and variances used in the TransCom 3 inversion are given in the online supplement to Gurney et al. [2002].) [9] The transport operator A for this inverse problem specifies the modeled impact of the magnitude of each of the TransCom 3 regional flux patterns on each data point. We modeled atmospheric CO 2 transport with the Model of Atmospheric Transport and Chemistry (MATCH) [Rasch et al., 1994;Mahowald et al., 1997], run at T21 resolution (about 5.5°* 5.5°) with 26 vertical levels. The model was driven with one year of winds from a simulation with the National Center for Atmospheric Research's Community Climate Model version 3 [Kiehl et al., 1998]. In this configuration, MATCH has been used to study column CO 2 variability [Olsen and Randerson, 2004]. Other configurations of MATCH participated in the TransCom model intercomparisons ].

Parameters Varied
[10] Like other recent inverse modeling studies, TransCom 3 uses prior estimates of regional CO 2 fluxes to stabilize the flux estimates. The weighting assigned to the prior estimates of regional CO 2 fluxes determines the relative influence of the data and of the prior estimates on the flux estimates. Assigning a very large weight to the prior estimates leads to a solution that is close to the prior estimates and fits the data poorly; a very small weight leads to an unstable solution that typically exhibits spurious largeamplitude fluxes.
[11] Here we scale the weighting given to the prior estimates by a parameter l that we vary from 0.5 to 3 (equation (2)). The parameter l can be interpreted as the regularization parameter in what is known as Tikhonov regularization or ridge regression [e.g., Hansen, 1998, chapter 5]. Alternatively, in a Bayesian inversion framework, the parameter l can be thought of as a scaling factor of the prior standard deviations, so that instead of the error variances c x,0 assigned by TransCom to the prior estimates, the error variances are of the form c x = l À2 c x,0 . (See also the online supplement.) [12] The relative magnitudes of the error variances c b assigned to the CO 2 concentration data control the relative weights of CO 2 concentrations at different locations in the inverse model. The data error variances should reflect uncertainty due to error in the transport model and error in the regional source patterns assumed in the inversion as well as measurement and calibration error [Kaminski et al., 2001]. (Estimated errors in the forward model can be taken into account explicitly in the solution of inverse problems [Van Huffel and Vandewalle, 1989;Golub et al., 1999]. Using such methods to estimate CO 2 fluxes appears to be worth pursuing, but we will not do so here.) [13] Investigators have typically either assumed the same error variance for all stations [Fan et al., 1998;Rayner et al., 1999] or have, as in TransCom, taken the error variance at each station to be proportional to the high-frequency variance of CO 2 concentrations at the station. In fact, the error variance may well be the sum of components that scale in each of the two ways [cf. Rödenbeck et al., 2003]. Therefore, we allow a gradation from uniform station error variances c b,e to the TransCom error variances c b,0 by varying a parameter t between 0 (uniform error variances) and 1 (TransCom error variances). The i-th station variance is given by c b i = (c b,0 i ) t (c b,e i ) 1Àt , and the magnitude of the uniform error variance is chosen such that the diagonal covariance matrices C b,e and C b,0 have the same determinant.

Effects of Parameter Choices on Flux Estimates
[14] The large northern land sink inferred in the TransCom 3 inversions is stable with respect to changing the parameter values, decreasing by $25% (0.7 Pg C yr À1 ) with increasing l (Figure 1a). The longitudinal distribution of the sink changes as t is varied (an effect noted by Gurney et al. [2002] as explaining part of the difference between the TransCom flux distribution and that inferred by Fan et al. [1998]). Most of this shift is due to two stations at the ITN tower in North Carolina [Bakwin et al., 1995], whose low CO 2 concentrations imply a large North American sink but which are assigned large errors, hence little weight, in TransCom. At low t, the North American sink is inferred to be higher by 0.3 Pg C yr À1 (Figure 1b) and the estimated sinks in Europe (Figure 1c) and Asia correspondingly decrease.
[15] By contrast, the estimated flux distribution in equatorial and southern regions depends strongly on the parameter values chosen. Changing l can shift fluxes of 2 Pg C yr À1 between equatorial and southern land regions (Figures 1d and 1e), with particularly large flux magnitudes at low l, where the prior constraints are weaker.
[16] The sensitivity of the equatorial and southern regions' flux estimates to parameter choices affects the global partitioning of fluxes between land and ocean. The total ocean flux estimate varies by up to 0.8 Pg C yr À1 as t is changed (Figure 1f), with compensating changes in the land flux estimate. (Because of the global growth rate constraint, the sum of land and ocean fluxes stays approximately constant at À2.8 Pg C yr À1 as inversion parameters are varied.)

Choosing Parameters by GCV
[17] The GCV function is minimized when t = 0.7 and l = 2.1 (Figure 2). TransCom's error model (t = 1) or an intermediate error model (0 < t < 1) thus seems more appropriate than an error model with equal error variances for all stations (t = 0). The GCV choice of l implies that the TransCom standard deviations of the prior fluxes may be too large by about a factor of 2, intimating that TransCom's estimated CO 2 fluxes overfit the data.
[18] For northern regions, using the GCV parameter choices in place of the TransCom parameter values changes the estimated fluxes by less than 20% (Figures 1a -1c). For many equatorial and southern regions (Figures 1d and 1e), however, the GCV flux estimates are much closer to the prior values than the TransCom estimates, suggesting that the station network provides little information about CO 2 fluxes in these regions.

Implications for Future Inverse Modeling Studies
[19] Parameter-choice methods such as GCV offer some scope for reducing uncertainty in and variations between inverse modeling studies. Although selecting what parameters to estimate will remain arbitrary, the fact that even our restricted variation of parameters resulted in differences in estimated fluxes that are comparable to the TransCom inter-model variation and to the a posteriori uncertainties estimated from the data and prior flux uncertainties [Gurney et al., 2002] suggests that choosing parameters systematically is necessary.
[20] GCV could be particularly useful for choosing relative weights for blocks of data from different sources [Gao, 1994]. As inversions for carbon fluxes add additional types of data, such as CO 2 isotopic composition [Ciais et al., 1995;Rayner et al., 1999], spectroscopically determined column CO 2 [Yang et al., 2002], and measurements of other gases such as CO and O 2 , GCV can be used to select a weighting for each type that reflects not only measurement accuracies but also model capability to simulate each type of data.

Conclusion
[21] The choice of parameters can significantly contribute to variations in CO 2 flux estimates obtained in inverse modeling studies. Methods such as GCV that choose parameters systematically by optimizing a given objective  function can improve inversion results. We recommend that uncertainty in inversion parameters be considered in future inverse modeling protocols and that formal parameter choice methods be used where appropriate.