Regional Changes in Carbon Dioxide Fluxes of Land and Oceans Since 1980 Philippe Bousquet, Philippe Peylin, Philippe Ciais, Corinne Le Quéré, Pierre Friedlingstein, and Pieter P. Tans

Supplementary Material

Calculation of the Deseasonalized Trend from the Raw Data
In Fig. 1 of the report, we present the variance of the deseasonalized trend at each measurement site. To compute this quantity, the raw data are fitted using the sum of a second-order polynomial function and four annual harmonics. Residuals are then defined as the difference between the raw data and the fitted curve. The deseasonalized trend is calculated as the sum of the same second-order polynomial function and the smoothed residuals using a low-pass filter [an effective width of 1 year; see (1)].

Method and Control Inversion
We have developed an inverse model to retrieve the net CO₂ fluxes every month from 1980 to 1998. The method extends the work of (2). We divide the land surface into 11 regions and the ocean surface into eight regions (Fig. 1 of the report), with regional fluxes being assigned a priori monthly values, monthly uncertainties, and spatial patterns (2). Land fluxes are assigned prior values as in previous inversions and prior geographic patterns within each region from the SIB-2 global biosphere model (3, 4). Ocean fluxes are assigned prior values and patterns from a global synthesis of air-sea flux measurements (5). Prior errors on net fluxes over each region (ocean or land) are set to ±1 GtC year^-1. Such loosely defined prior errors do not nudge the solution of the inversion to the prior estimates. We next calculate the atmospheric CO₂ distribution caused by atmospheric transport acting on a pulsed source of 1 GtC emitted at a constant rate by region i during month j. We do this for all regions in a 3D global transport model, and archive the resulting CO₂ concentration field for the following 2 years (6). This concentration pattern, initially with a maximum over the source region, becomes progressively homogeneous in the atmosphere, as CO₂ emitted over the source region gets diluted globally by the atmospheric transport. Let B_ij(x), O_ij(x), and F_ij(x) denote respectively the land, oceanic, and fossil CO₂ concentration patterns induced by region i at month j. The modeled spatial CO₂ pattern P caused by all sources (N regions and N months) is

	Nregions	Nmonths
P =	S	S	(bijBij + wijOij + fijFij)
	i = 1	j=1

The coefficients f_ij, corresponding to the regional magnitude of fossil emissions, are set to fixed monthly values, based on interannual fossil fuel emissions statistics. Note that the value of b_ijreflects the sum of land use and of other biospheric carbon sources and sinks. As changes in the land use-induced CO₂ emissions over 1980-1998 are poorly documented, they were not subtracted from the inversion results.

We solve for the coefficients b_ij and w_ij in order to minimize a cost function J. This cost function comprises the sum of the distances between model responses and observations, and the sum of the distances between a priori flux estimates and optimized flux estimates. All distances are weighted by the inverse of the a priori variances respectively on observations and fluxes. We assume all uncertainties to be Gaussian and to be uncorrelated spatially and temporally. Errors on the data are calculated monthly at each site either from the standard deviation of raw flask measurements (if more than four individual flasks are sampled during a month), or from the mean standard deviation GLOBALVIEW-CO₂ (if not). A minimum error value of ±0.5 ppm is also applied at all sites. Additional constraints to limit unrealistic monthly fluctuations of the source magnitudes are also used as in (2). We also add a constraint on the global oceanic uptake of 2.0 ± 1.2 GtC year^-1, based on O₂/N₂ data (7). The choice of 1.2 GtC year^-1 for the error on this global ocean uptake is the maximum error reported by (7).

Solving for 20 years of atmospheric data (1979-1998) implies the inversion of a 25,000 × 8000 matrix. To reduce the size of the matrix, we apply a Choleski decomposition to 6-year subsets of the CO₂ time series. We next drop the two first and the two last years of each estimated flux time series. Such truncation avoids jumps in the estimated fluxes between adjacent 6-year subsets. This approach relies on the hypothesis that the concentration of a tracer emitted at one point on Earth's surface requires 2 years to become homogenized within the troposphere. With such a hypothesis, sources of a given month do not influence atmospheric concentrations 2 years after their emission.

Net Fluxes

At a global scale, we infer a long-term mean uptake of 1.3 ± 0.8 GtC year^-1into the land biosphere and of 1.4 ± 0.7 GtC year^-1 into the oceans. These values are in the range of the estimates given by the IPCC (8). Note that for the oceans, the IPCC uptake of 2.0 ± 0.8 GtC year^-1 does not account for the CO₂ source due to transport by the river, which is estimated to be an additional source of 0.4 GtC by Aumont (9). The uncertainties of our flux estimates are calculated as the square root of "overall variances." For each flux, we calculate the overall variance as the sum of (a) the residual variance given by the a posteriori covariance matrix of the control inversion, and (b) the variance produced by the eight inversions performed. We thus assume that these two types of errors are independent.

Web table 1 below gives the long-term mean fluxes over the 1980s and 1990s for four land regions and four ocean regions (fossil fuel emissions are subtracted). Globally, we find an increase of around 30% for both ocean and land uptakes between the 1980s and the 1990s. Regionally, Eurasia is found to be a mean large sink over the two decades, whereas North America is rather neutral. However, the inferred "overall variances" are very large (~1.0 GtC year^-1), which indicates that the partition of the long-term mean Northern Hemisphere sink is not safely retrieved by the inversions, contrary to flux anomalies. Other inverse studies estimating regional CO₂ net fluxes show conflicting results. For instance, for the 1988-1992 period, Fan et al. (10) estimate a large terrestrial carbon sink in North America of 1.7 ± 0.5 GtC year^-1, whereas we infer a small source of 0.1 ± 0.9 GtC year^-1. An analysis of those differences is currently under investigation (11). Tropical lands are found to be a mean net sink of 0.6 GtC year^-1 in the 1980s and of 0.3 GtC year^-1 in the 1990s, with large uncertainties (above 1.0 GtC year^-1). For the oceans, we infer net fluxes with smaller uncertainties than for the land fluxes. This partly reflects the denser network over the oceans than over the continents. In our inversions, the increase of the global ocean sink between the 1980s and the 1990s is mainly due to a reduction of 50% of the tropical ocean source.

Supplemental Table 1. Long-term mean fluxes inferred by our eight inversions for four land regions and four ocean regions.
LANDS	1980-1989	1990-1998	OCEANS	1980-1989	1990-1998
North America	0.2 ± 0.9	-0.1 ± 0.8	North Atlantic	-0.7 ± 0.5	-0.8 ± 0.4
Eurasia	-1.1 ± 1.2	-1.3 ± 1.1	North Pacific	-0.4 ± 0.4	-0.5 ± 0.4
Tropical lands	-0.6 ± 1.2	-0.3 ± 1.1	Tropical oceans	0.9 ± 0.5	0.4 ± 0.4
Nontropical SH lands	0.5 ± 0.8	0.3 ± 0.7	Nontropical SH oceans	-0.9 ± 0.3	-0.9 ± 0.3
Total lands	-1.0 ± 0.9	-1.4 ± 0.8	Total oceans	-1.1 ± 0.7	-1.8 ± 0.6

---

Multivariate Regression Coefficient and Time Lag Calculation
In order to evaluate which fluxes dominate year-to-year variations, we have calculated a regression coefficient R_m. R_m is the coefficient of a multivariate regression of the total flux anomaly (global land + global ocean) against the time series of both land and ocean fluxes for three latitudinal bands (north of 20°N, 20°N to 20°S, and south of 20°S). Doing so, we find that tropical land surfaces strongly dominate the total flux anomaly in the 1980-1989 period (R_m > 0.99). Northern hemisphere land surfaces slightly dominate the total flux anomaly between 1991 and 1995 (R_m = 0.75). In the 1996-1998 period, tropical land surfaces and Northern hemisphere land surfaces contribute almost equally to the total flux anomaly (R_m = 0.6).

We have computed time lag correlations between the time series of the mean CO₂ growth rate for tropical stations (20°S to 20°N) and the time series of the mean CO₂ growth rate for stations poleward of 40°N. For each month of the period 1980-1998, we search for the time lag that produces the maximum correlation over a restricted period of 2 years centered at the current month. Such optimal time lags, if different from zero, indicate anomalies in the CO₂ growth rate that propagate from the tropics to the high latitudes or vice versa.

North American Carbon Balance
Changes in the difference of CO₂ concentrations between Atlantic and Pacific stations correlate with inferred changes in the North American carbon balance. This correlation is a necessary, albeit not sufficient, verification of the inversion results. One could argue that a negative difference between Atlantic and Pacific stations could reflect an anomalous source over Eurasia, or an anomalous sink over the Atlantic Ocean, rather than just an anomalous sink over North America. We have carried out an additional inverse calculation in which the North American flux anomaly is forced to be zero over the past 20 years. To match the observed changes in the difference between the Pacific and the Atlantic, that inversion generated an anomalous source (+1 GtC) over Eurasia in 1989-1990 combined with an anomalous sink (1.5 GtC) over the northern oceans. The latter flux would translate into an increase in DpCO₂ by a factor of 2 within 12 months poleward of 20°N. For relatively well-sampled areas such as the North Atlantic, the statistical uncertainty of the DpCO₂ climatology for the two northern basins is as small as ±0.1 PgC year^-1 (12), making this inversion result very unlikely.

References and Notes

1. K. W. Thoning, P. P. Tans, W. D. Komhyr, J. Geophys. Res. 94, 8549 (1989).

2. P. Peylin, P. Bousquet, P. Ciais, P. Monfray, in Inverse Methods in Global Biogeochemical Cycles, P. Kasibhatla et al., Eds. (American Geophysical Union, Washington, DC, 1999).

3. P. J. Sellers et al., J. Clim. 9, 676 (1996).

4. S. Denning et al., Tellus 48B, 521 (1996).

5. T. Takahashi et al., Proc. Natl. Acad. Sci. U.S.A. 94, 8292 (1997).

6. Pulses of regional monthly sources are prescribed into the TM2 global 3D atmospheric transport model (M. Heimann, Report 10, Max-Planck-Institut für Meteorologie, Hamburg, 1995). The horizontal resolution of the model is 7.5° × 7.5° and the air column is divided into nine levels in sigma coordinates. The model output every 3 hours is sampled at each station and averaged into monthly means. Analyzed ECMWF meteorological parameters for year 1990 are used as an input of TM2 model in the control inversion. Year 1993 is used in one case as a sensitivity study. We do not model explicitly year-to-year changes in the atmospheric transport over 1980-1998. A recent study [R. J. Dargaville, R. M. Law, F. Pribac, Global Biogeochem. Cycles 14, 931 (2000)] has shown that the impact of year-to-year variations in atmospheric transport on inversion results is small compared to year-to-year changes in sources and sinks.

7. M. Battle et al., Science 287, 2467 (2000).

8. International Panel on Climate Change, Climate Change 1995. The Science of Climate Change (Cambridge Univ. Press, Cambridge, 1995).

9. O. Aumont, thesis, Université Paris VI, France (1998).

10. S. Fan et al., Science 282, 442 (1998).

11. P. Peylin, in preparation.

12. T. Takahashi, personal communication.