Detection and attribution of global change and disturbance impacts on a tower-observed ecosystem carbon budget: a critical appraisal

Observations worldwide are providing an increasing amount of atmosphere–ecosystem flux data. Thus, the establishment of a data mining methodology to detect significant trends and attribute changes to specific factors is important. This study examined the possibility of detecting significant trends in observed data at a test site with one of the longest records of flux measurements (Takayama, Japan). Statistical tests using non-parametric methods showed a ‘likely’ trend (i.e., detected at 66–90% confidence level) of increasing carbon sequestration. To investigate the change in carbon sequestration in relation to biological and environmental factors (ambient CO2, temperature, radiation, precipitation and disturbance), mechanistic and numerical methods were applied. A process-based model was used for the mechanistic attribution of change, and an optimal fingerprinting method in combination with model-based sensitivity simulations was used for numerical attribution. At the study site, local disturbances appeared to exert an impact on the observed carbon sequestration, whereas climatic factors made moderate contributions. These results indicate the feasibility of detection and attribution using current flux measurement data, although more evidence is needed to confirm global coherence.


Introduction
Elevated atmospheric CO 2 and anthropogenic climatic change exert various degrees of impact on natural and agricultural ecosystems (IPCC 2007) and propagate to human society through the degradation of ecosystem services, such as the provision of material and genetic resources, and the regulation of biogeochemical and biogeophysical processes. However, the detection and attribution of these impacts through observation alone are not easy as (1) these impacts may occur heterogeneously in space and time; (2) natural variability, such as year-to-year weather changes, conceals trends; and (3) local factors, such as stochastic natural disturbances and human intervention, affect ecosystems at the same time.
Continuous measurements of the atmosphere-ecosystem exchange of CO 2 are conducted at more than 500 sites, mainly using the eddy covariance method (FLUXNET;Baldocchi et al 2001). Flux data provide indispensable evidence of the current state of the terrestrial carbon budget, including that of croplands and urban areas, and provide implications about the variability of environmental change and mechanisms responding to it (e.g., Jung et al 2010, Yi et al 2010. However, the main difficulty in using the data to detect and attribute the impacts of global warming on the carbon budget is the brevity of the observation period, as the eddy covariance method of flux measurement was standardized in the 1990s, and most sites have an observation length of less than 10 yr. The longest measurement period is that of the Harvard Forest (Petersham, MA, USA), which began in October 1990 (Wofsy et al 1993, Urbanski et al 2007. In many cases, observation records contain gaps due to accidents or difficult conditions, and continuous (i.e., gap-filled) data reflect the gap-filling method. Nevertheless, examining the possibilities and difficulties of detecting and attributing the impact of global warming on ecosystem carbon budgets is important to obtain early warnings from observations and to derive evidence of climatic feedback from the terrestrial carbon cycle. The detection of global warming phenomena (e.g., temperature rise) and their impacts on social and natural systems is a prerequisite to implementing effective measures to mitigate and adapt to global change. However, due to the chaotic nature of the climate system and scarcity of observational data, chronological records still present uncertainties in trend detection (e.g., Brohan et al 2006). Few studies have attempted to detect ecosystem-scale impacts of global warming on the terrestrial carbon budget although several model studies (e.g., Piao et al 2009) have estimated the impacts of various factors at a regional scale. The attribution of the detected changes to individual factors is also an important issue in understanding global warming impacts. However, quantitative attribution is generally difficult, as multiple factors exert influences at the same time, and these factors are not always independent of one another.
The primary objective of the present study was to examine the current status (possibilities and difficulties) of detecting and attributing global warming impacts on a net ecosystem CO 2 budget, which is inevitably affected by local factors such as disturbance, using tower-observed flux data.

Overview
A case study was conducted at a flux measurement site (Takayama) with one of the longest observation records in the FLUXNET. Using the 16 yr (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) time series of CO 2 fluxes, several statistical methods commonly used in timeseries analyses were applied to trend detection. This study examined annual total values, and discussed the observation length necessary to detect a trend at a certain confidence level. Using process-based model simulations and a numerical algorism (optimal fingerprinting), I examined the attribution of the following observed changes to environmental factors: elevated CO 2 , temperature, precipitation, radiation and disturbance.

Description of site and flux data
This study was conducted at a cool-temperate deciduous broad-leaved forest in Takayama, central Japan: 36 • 08 N, 137 • 25 E, 1420 m above sea level. Annual mean temperature and annual precipitation during the period 1980-2002 were 7.2 • C and 2275 mm, respectively (Ohtsuka et al 2009; see figure S1 (available at stacks.iop.org/ERL/7/014013/mmedia) for interannual variability). The forest is mainly composed of the hardwood species: Quercus crispula Blume, Betula ermanii Cham and Betula platyphylla Sukaczev var. japonica Hara. The forest floor is densely covered with dwarf bamboo grass, as well as herbs and shrubs. Tree density was 1051 trees ha −1 in 2007, with a canopy height of 13-20 m. Leaf area indices of the tree canopy and understory were 5 and 1.6 at peak season, respectively (Nasahara et al 2008, Ohtsuka et al 2009. This secondary forest is approximately 40-50 yr old, and is recovering from a logging disturbance ca. 1960 (typical disturbance history in the area). The soil is forest brown soil with rich amounts of organic material and volcanic ash.
A 27 m tall flux tower that uses aerodynamic and eddy-covariance methods for continuous measurement was constructed in 1993 (Saigusa et al 2005). The fetch length of the forest is 500-1000 m, and the tower footprint is uniformly covered with hardwood species. Wind speed and direction are measured using a three-dimensional ultrasonic anemometer (DT-600; Kaijo, Japan) and atmospheric water vapour and CO 2 concentrations are measured using a closed-path infrared gas analyser (LI-6262; LI-COR, Lincoln, NE, US). Raw data at 5 Hz are corrected for spikes, wind rotation and air density, and used to calculate heat, vapour and CO 2 fluxes at 30 min intervals. Hirata et al (2008) filled observation gaps using the non-linear regression method, in which data at insufficient atmospheric turbulence conditions (friction velocity < 0.2 m s −1 ) were screened out. Gap-filled data from 1994 to 2009 were used.

Model description
This study used the process-based terrestrial model: Vegetation Integrated Simulator for Trace gases (VISIT; Inatomi et al 2010, Ito 2010. The model consists of submodules simulating radiation budget, hydrology, phenology, and carbon and nitrogen cycles. The ecosystem structure is simplified into four sectors of carbon stock: canopy trees, floor plants, dead biomass (litter) and mineral soil (humus). Each carbon stock sector is then divided into several carbon pools; for example, the tree sector is composed of leaf, stem and root pools (see Ito et al 2005 for a schematic diagram). Photosynthetic CO 2 assimilation, which is gross primary production (GPP) at the ecosystem scale, is estimated as a function of leaf area index (LAI; estimated from the leaf carbon pool and specific leaf area), canopy light absorption coefficient and leaf-level photosynthetic parameters (Ito and Oikawa 2002). Plant and soil respiration rates are estimated from the pool-specific respiration rate, pool size and response to temperature. Total respiration is called ecosystem respiration (RE). Net ecosystem production (NEP; equivalent to net CO 2 exchange with the atmosphere) is calculated as the difference between GPP and RE.
At the Takayama site, the VISIT model has been examined using various observational data and has been employed for biogeochemical studies. Ito et al (2005) compared the simulated gross and net CO 2 fluxes with flux measurement data; Ito et al (2007) compared the simulated soil CO 2 efflux with chamber-based data; Inatomi et al (2010) compared the simulated CH 4 /N 2 O fluxes with chamber-based data; and Ito (2010) compared the simulated carbon stock and LAI with biometric data (Ohtsuka et al 2009) and satellite data, including the impacts of cyclone-induced defoliation. These validations have confirmed that the model captures atmosphere-ecosystem exchange with high precision. Inatomi et al (2010) used time-series data of nitrogen deposition to simulate N 2 O fluxes at the site, and confirmed that both nitrogen input and gaseous efflux did not increase significantly during the study period. Therefore, nitrogen factor was not addressed in this study. For model-simulated daily CO 2 fluxes, the coefficient of determination (r 2 ) was 0.851, and the root mean square error (RMSE) of NEP was 0.674 g C m −2 d −1 (see figure S2 (available at stacks.iop.org/ERL/7/014013/ mmedia) for a comparison between observed and simulated monthly fluxes in 1999-2009).

Detection
Non-parametric statistical methods were used to determine significant trends in the observed time-series flux data at a certain (typically, 95%) confidence level without assuming any particular distribution: Sen's trend analysis and the Mann-Kendall test. These methods were applied to GPP, RE and NEP time-series data observed at the study site and also to environmental variable data. Sen's robust slope estimator (Sen 1968) was defined by the following: where y i and y j are data at time i and j, respectively. The Mann-Kendall rank statistic test (Kendall 1938, Mann 1945) is a non-parametric trend detection method that is flexibly applicable to data with non-Gaussian error. The null hypothesis for the Mann-Kendall test is that the data are independent and randomly ordered (i.e., no trend), and rank statistics (S and τ ) were calculated. The discussion of the necessary length of a time series to detect a trend at a certain confidence level is based on a method often used in climate data analyses (e.g., Kiesewetter et al 2010). Tiao et al (1990) established that the necessary data length is affected by the slope (ω), noise variance (σ N ), and autocorrelation (φ, time-shift by −1 step) of the data. Weatherhead et al (1998) derived the following equation to calculate the number of years (n * ) of data required to detect a real trend with a given confidence level: where α is a coefficient of the prescribed confidence level (e.g., 3.3 for 90% probability). It is assumed that |φ| < 1 and that σ N and ω are on the same scale. Intuitively, this equation indicates that longer records are required to detect a significant trend, when data have a weaker slope, higher variability and stronger autocorrelation.

Attribution
A mechanistic and a numerical method were used to separately evaluate the effects of radiation, temperature, precipitation, atmospheric CO 2 concentration and disturbance. Both methods used 'constrained' model simulations, in which the temporal variability of the specific factors was individually considered. Thus, five simulations of the five environmental factors were separately conducted. In each case, the other four factors were assumed to be temporally static (i.e., long-term average climate or no disturbance condition).
In the mechanistic attribution, the annual carbon budgets (GPP, RE and NEP) of the 'constrained' simulations of each factor were compared with that of a control that included all factors. Note that in reality and in the model, environmental factors are not independent of one another. Thus, the sum of the contributions of the five factors was not necessarily identical to the control result. For example, severe disturbance would affect ecosystem responsiveness to other environmental factors, and meteorological variables tend to co-vary. Here, the residue is called a non-linearity term.
The numerical attribution was conducted using the optimal fingerprinting method (Allen and Stott 2003), which has been used to separate the contribution of multiple factors on observed climate change. In the method, the observation (y) is represented using a linear model, where x k is the model simulation considering the kth factor (K = 5 for temperature, precipitation, radiation, CO 2 and disturbance), β k are scaling coefficients, and v represents noise in both the observation and model. The optimal values and uncertainty range of β k were obtained using the Monte Carlo approach, so that variance of the noise term v is minimized (Allen and Stott 2003). Using the simplified Metropolis-Hasting algorithm (Hastings 1970), we estimated the mean and variance of β k through 10 000 repetitions, each of which started from a different initial value. In each trial, optimal values of parameters, i.e., mean and variance of β k , were searched by adding random perturbation terms (200 000 trials) and comparing the variance of v to determine adoption or rejection of the updated parameters. The average of the 10 000 estimates should represent the strength and uncertainty of each environmental factor.

Detection
The most significant trend shown by the Mann-Kendall test was annual precipitation at the site and secondarily mean temperature at a nearby weather station. The increasing carbon budget trends were detected by the Mann-Kendall test at confidence levels of 80-90% (table 1). The IPCC Assessment Report (IPCC 2007) termed confidence levels of 66-90% as 'likely', 90-95% as 'very likely', and >95% as 'extremely likely'. In this regard, the precipitation trend is 'very likely', and the increasing trend of carbon sequestration is 'likely'. It should be noted that the three methods did not always agree in the significance of trend detection. I calculated the data length necessary to detect real trends (equation (2)) in the observed carbon fluxes. Detecting a trend in the observed GPP record at 90%, 95% and 99% confidence levels may require 11.3, 14.2 and 15.2 yr, respectively. The data contain small autocorrelation (φ = 0.029) and low variability (σ N = 84.9 g C m −2 yr −1 ); thus, current observations may be adequate to detect a trend in GPP. On the other hand, the linear trend in RE was less evident (φ = −0.11 and σ N = 30.2 g C m −2 yr −1 ), necessitating longer observation periods of 14.9, 18.7 and 20.0 yr at 90%, 95% and 99% confidence levels, respectively. The trend in NEP was intermediate between those in GPP and RE (φ = 0.032 and σ N = 69.1 g C m −2 yr −1 ), thus implying that data lengths of 11.4, 14.6 and 15.5 yr are required for 90%, 95% and 99% confidence levels, respectively.

Attribution
The mechanistic (i.e., process-model-based) attribution separated the contributions of environmental factors to the carbon budget change at the Takayama site. Figure 2 shows that carbon sequestration was primarily attributable to disturbance (62% of total NEP in 1994-2009, 277 g C m −2 yr −1 ) and secondarily to environmental changes (5.8% to temperature change and 11% to elevated CO 2 ). In contrast, the change in precipitation did not contribute to carbon sequestration at the study site because this site has a sufficiently humid (i.e., drought free) monsoon climate, with annual precipitation of about 2200 mm. The decreasing trend in solar radiation variation (table 1) contributed to carbon release, offsetting part of the carbon sequestration caused by other factors. As a result of the interaction among factors (non-linearity, figure 2), an additional 86 g C m −2 yr −1 (equivalent to 31% of total NEP) of net carbon sink occurred.
In the numerical attribution, which used the optimal fingerprinting method, a scaling coefficient (β k ) close to 1 indicates that the factor made a large contribution to the observed change. This analysis revealed that disturbance was the largest contributor to temporal GPP change and, secondarily, that changes in temperature and precipitation contributed equally (figure 3). Unexpectedly, elevated CO 2 showed a negative scaling coefficient for GPP, whereas the CO 2 fertilization effect (included in the model) should have increased GPP in response to elevated CO 2 . The estimated scaling coefficient for the effect of elevated CO 2 on GPP had a wide range of uncertainty (99% confidence range: −1.9-3.1). The change in RE was primarily attributed to elevated CO 2 and temperature change, whereas disturbance made only a small contribution. A previous study by Ito et al (2007) indicated that soil respiration, which should be more stable against disturbance compared to aboveground respiration, accounted for 67% of RE at the site. In contrast to GPP, elevated CO 2 showed the largest scaling coefficient (1.7) for NEP, with a wide range of uncertainty and a relatively small range (0.4) for the effect of disturbance on NEP.

Discussion
This study examined the feasibility of detecting a significant trend in long-term observed CO 2 fluxes and attributing this trend to changes in specific environmental and disturbance factors. Statistical methods enabled us to detect a 'likely' trend in the 16 yr time series of observational data at the Takayama site, although finding a clear trend in the data was difficult due to natural stochastic variability (e.g., year-to-year weather) and local disturbances. If the disturbance impact is subtracted from the carbon budget change (62% of NEP, figure 2), longer observational periods should be required: e.g., calculated by equation (3) as 23.3 yr at 90% confidence level. Compared with meteorological records (e.g.,   (3)) estimated by the optimal fingerprinting method for the carbon budget change at the Takayama site. air temperature for >100 yr), ecosystem flux data have shorter coverage (∼20 yr) but less temporal variability: the GPP and RE coefficients of variance were just 9.7% and 4.2%, respectively. Hence, standard statistical methods could be applied to detect a significant linear trend in the data. The mechanistic attribution using appropriate process-based models enabled us to separate the contributions of different factors to observed change, although the results depended on the environmental sensitivity assumed in the model. On the other hand, the numerical attribution results appeared to contain discrepancies (figure 3), which suggested that the method may have limited applicability in separating multiple factors that have the same tendencies (i.e., parallel increasing trends caused by disturbance and elevated CO 2 ). Indeed, the numerical method does not consider ecophysiological and biogeochemical consistency, but does imply the magnitude of estimation uncertainty.
In addition to this case study of the Takayama site, I additionally conducted a similar detection analysis of Harvard Forest flux data (Urbanski et al 2007) from 1992 to 2006 to examine how well the statistical tests capture site-specific biases. Each of the three tests detected an 'extremely likely' (>95% probability) trend in GPP and NEP (table S1 available at stacks.iop.org/ERL/7/014013/mmedia); note that this forest is ca. 80 yr old and has been influenced by several hurricanes. These results also support the feasibility of deriving long-term signals from existing flux measurement data, although perturbation induced by local factors, such as disturbances and nitrogen deposition (e.g., Magnani et al 2007), must be addressed. A further accumulation of evidence would allow us to investigate the coherence of the impacts of elevated CO 2 and climate change on the terrestrial carbon budget.
Among the insights gained from these detection and attribution efforts were: (1) simulation studies using coupled climate-carbon cycle models imply that terrestrial ecosystems may exert biogeochemical feedbacks on anthropogenic climate change (e.g., Cox et al 2000). The detection and attribution results of this study may provide information about the direction of terrestrial carbon budget changes for simulation studies.
(2) Current countermeasures to global warming rely considerably on carbon sequestration by terrestrial ecosystems, such as forest carbon sinks and reduction of emissions from deforestation and forest degradation. However, they do not pay enough attention to temporal changes in forest carbon budgets as a result of global change. Detection and attribution studies may allow us to correct the amount of forest carbon sequestration potential.
(3) Detections studies such as Parmesan and Yohe (2003) require extensive data and analyses to determine coherent signals of climate change impacts on natural systems. In contrast, the present study exemplifies how meaningful trends can be derived from single-site data by using appropriate data mining methods. (4) Finally, this study indicates a new opportunity in data-model fusion of long-term ecosystem dynamics, in addition to short-term biophysical and physiological processes.