资源与生态环境（文献资料）Model–data synthesis in terrestrial carbon observation - methods, data requirements and data uncertainty specifications

Global Change Biology(2005)11,378-397doi: 10.1111/j.1365-2486.2005.00917.xModel-data synthesis in terrestrial carbon observation:methods, data requirements and data uncertaintyspecificationsM.R.RAUPACH,P.J.RAYNERt,D.J.BARRETT,R.S.DEFRIESS,M.HEIMANN,D.S.OJIMAll,S.QUEGAN**andC.C.SCHMULLIUStt*CSIRO Earth Obseruation Centre,GPO Box 3023, Canberra,ACT 2601,Australia,LSCE-CEA de Saclay Orme des Merisiers,91191Gif/Yoette,France,CSIROLand andWater,Canberra,ACT2601,Australia,SDepartmentofGeography,UniversityofMaryland,College Park,MD20742,USA,Department of"Biogeochemical Systems",Max-Planck-Institut fuir Biogeochemie,D-07701Jena,Germany,NaturalResourceEcologyLaboratory,ColoradoStateUniversity,FortCollins,CO80523-1499,UuSA,**CentreforTerrestrialCarbonDynamics,UniversityofSheffield,SheffieldS37RH,UK,+InstitutefurGeographie,Friedrich-Schiller-Universitat,D-07743 Jena, Germany.AbstractSystematic, operational, long-term observations of the terrestrial carbon cycle (includingits interactions with water, energy and nutrient cycles and ecosystem dynamics) areimportant forthe prediction and management of climate, waterresources,foodresources, biodiversity and desertification. To contribute to these goals, a terrestrialcarbon observing system requires the synthesis of several kinds of observation intoterrestrial biosphere models encompassing the coupled cycles of carbon, water, energyand nutrients.Relevant observations include atmospheric composition (concentrationsof CO2and othergases); remotesensing;fluxandprocessmeasurementsfromintensivestudy sites; in situ vegetation and soil monitoring; weather, climate and hydrologicaldata; and contemporary and historical data on land use,land use change and disturbance(grazing,harvest, clearing,fire).A review of model-data synthesis tools for terrestrial carbon observation identifies'nonsequential' and 'sequential'approaches as major categories, differing according towhether data are treated all at once or sequentially. The structure underlying bothapproaches is reviewed, highlighting several basic commonalities in formalism and datarequirements.An essential commonality is that for all model-data synthesis problems,bothnonsequential and sequential, data uncertainties are as important as data valuesthemselves and have a comparable role in determining the outcome.Given the importance of data uncertainties, there is an urgent need for soundly baseduncertainty characterizations for the main kinds of data used in terrestrial carbonobservation.The first requirement is a specification of the main properties of the errorcovariance matrix.As a step towards this goal, semi-quantitative estimates are made of the mainproperties of the error covariance matrix for four kinds of data essential for terrestrialcarbonobservation:remote sensing ofland surfaceproperties,atmosphericcompositionmeasurements, direct flux measurements, and measurements of carbon stores.Received 8 June 2004; accepted 25August 2004Correspondence: Dr Michael Raupach, tel. + 61 2 6246 5573,e-mail:Michael.Raupach@csiro.au3782005Blackwell Publishing Ltd

Model–data synthesis in terrestrial carbon observation: methods, data requirements and data uncertainty specifications M. R. RAUPACH *, P. J . R AY NE R w , D. J. BARRETT z, R. S. DEFRIES§, M. HEIMANN } , D. S. OJIMA k, S. QUEGAN ** and C . C . S C H M U L L I U S w w *CSIRO Earth Observation Centre, GPO Box 3023, Canberra, ACT 2601, Australia, wLSCE-CEA de Saclay Orme des Merisiers, 91191 Gif/Yvette, France, zCSIRO Land and Water, Canberra, ACT 2601, Australia, §Department of Geography, University of Maryland, College Park, MD 20742, USA, }Department of ‘‘Biogeochemical Systems’’, Max-Planck-Institut fu¨r Biogeochemie, D-07701, Jena, Germany, kNatural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO 80523-1499, USA, **Centre for Terrestrial Carbon Dynamics, University of Sheffield, Sheffield S37RH, UK, wwInstitute fu¨r Geographie, Friedrich-Schiller-Universita¨t, D-07743 Jena, Germany. Abstract Systematic, operational, long-term observations of the terrestrial carbon cycle (including its interactions with water, energy and nutrient cycles and ecosystem dynamics) are important for the prediction and management of climate, water resources, food resources, biodiversity and desertification. To contribute to these goals, a terrestrial carbon observing system requires the synthesis of several kinds of observation into terrestrial biosphere models encompassing the coupled cycles of carbon, water, energy and nutrients. Relevant observations include atmospheric composition (concentrations of CO2 and other gases); remote sensing; flux and process measurements from intensive study sites; in situ vegetation and soil monitoring; weather, climate and hydrological data; and contemporary and historical data on land use, land use change and disturbance (grazing, harvest, clearing, fire). A review of model–data synthesis tools for terrestrial carbon observation identifies ‘nonsequential’ and ‘sequential’ approaches as major categories, differing according to whether data are treated all at once or sequentially. The structure underlying both approaches is reviewed, highlighting several basic commonalities in formalism and data requirements. An essential commonality is that for all model–data synthesis problems, both nonsequential and sequential, data uncertainties are as important as data values themselves and have a comparable role in determining the outcome. Given the importance of data uncertainties, there is an urgent need for soundly based uncertainty characterizations for the main kinds of data used in terrestrial carbon observation. The first requirement is a specification of the main properties of the error covariance matrix. As a step towards this goal, semi-quantitative estimates are made of the main properties of the error covariance matrix for four kinds of data essential for terrestrial carbon observation: remote sensing of land surface properties, atmospheric composition measurements, direct flux measurements, and measurements of carbon stores. Received 8 June 2004; accepted 25 August 2004 Correspondence: Dr Michael Raupach, tel. 1 61 2 6246 5573, e-mail: Michael.Raupach@csiro.au Global Change Biology (2005) 11, 378–397 doi: 10.1111/j.1365-2486.2005.00917.x 378 r 2005 Blackwell Publishing Ltd

MODEL-DATASYNTHESISINTERRESTRIALCARBONOBSERVATION379Introductionservations of quantities which are not directly obser-vable (suchas carbon stores and fluxes overlarge areas)Systematicearth observation impliesthecollection andand (4) forecasting (prediction forward in time on theinterpretation of multiple kinds of data about thebasis of past and current observations).evolvingstate of the earth system acrosswidespatialThe present paper arose from a workshop held inSheffield, UK, 3-6 June 2003, to further the develop-domains and over extended time periods.Three factorshavecausedamassiveacceleration inearthobservationment of a Terrestrial Carbon Observation Systemactivities over recent years. The first is need: global(TCOS) with a particular emphasis on model-datasynthesis.Antecedents for this effort were (1)pre-change is raising issues - such as greenhouse-inducedclimate change, water shortages and imbalances, landliminary steps toward aTCOS (Cihlar etal.,2002a,b,c);(2)a wider conceptforan Integrated Global Carbondegradation, soil erosion, loss of biodiversity-whichrequire informed human responses at both global andObserving Strategy including atmosphere, oceans, landregional levels. Second, technological advances inand human activities (Ciais et al., 2004) and (3) thesensors, satellite systems and data storageand proces-research program of the Global Carbon Project (Globalsing capabilities aremakingpossible observationsandCarbon Project, 2003).interpretations which were out of reach only a fewThe paper is founded on three themes arising fromyears ago and unimaginablea few decades ago.Third,the Sheffield workshop. First, model-data synthesis,based on terrestrial biospheremodels constrained withthesynthesisofformerlydiscretedisciplinesintoaunified Earth System Science is driving new hypothesesmultiplekinds of observation,is an essential compo-nent ofaTCOS.Second,from the standpoint ofmodel-about the dynamics of the earth system and theinterconnectednessof itscomponents,includinghu-data synthesis, data uncertainties are as important asmans.Systematic earthobservation motivates and testsdata values themselves and havea comparablerole inthesehypotheses.determining the outcome.Third, and consequently,The focus of this paper is observation of the carbonthere is an urgent need for soundly based uncertaintycycle, and in particular its land-atmosphere compo-specifications for the main kinds of data used innents,as onepart of an integrated earth observationterrestrial carbon observation.These themes are devel-system. It is a significant part because of the couplingoped as follows:thenext section summarizes majorbetween the carbon cycle and theterrestrial cycles ofpurposes and attributes of a TCOS. Model-datawater,energy and nutrients,and the connections ofallsynthesis:methods'providesan overviewof model-these biospheric processes with global climate anddata synthesis in the context of terrestrial carbonhuman activities (Field & Raupach, 2004; Raupach et al.,observation,by brieflydescribing some of the main2004).The carbon cycle is integral to the growth andmethods, indicating their common characteristics, anddecay of vegetation, maintains the water cycle throughhighlighting the key role of data uncertainty.Model-transpiration and provides habitat for maintainingdata synthesis: examples' provides some examples.biodiversity. Thus, terrestrial carbon observation isData characteristics: uncertainty in measurement andimportantfor climate observation and prediction,forrepresentation' undertakes a survey of the uncertaintythe management of water resources, nutrients andcharacteristics of the main kinds of relevant data.biodiversity, and for monitoring and managing theenhanced greenhouse effect.Purposes and attributes of a TCOSIt is increasingly recognized that strategies for earthA succinct statement of the overall purpose of a TCOSobservation (including terrestrial carbon observation)require methods for combiningdata andprocessmight be:to operationallymonitor the cycles of carbonmodels in systematic ways.This is leading to researchand related entities (water, energy, nutrients) in thetowards the application in terrestrial carbon observa-terrestrial biosphere, in support of comprehensive,tion (and inearth observationmoregenerally)ofsustained earth observation and prediction,and hence'model-data synthesis', the combination of the infor-sustainable environmental management and socio-mation contained in both observations and modelseconomic development.These words are congruentthrough bothparameter-estimation and data-assimila-with theFrameworkDocument emerging fromtheSecondEarthObservation Summit,Tokyo,April 2004tion techniques.Motivations for model-data synthesisapproaches include (1) model testing and data quality(http://earthobservations.org/docs/Framework%20control (through systematic checks for agreementDoc%20Final.pdf),whichcallsforaGlobalEarthwithin specified uncertaintybands for bothdata andObservationSystem of Systems'toservenineareasofsocio-economic benefit.ATCOSis a contributor to suchmodel); (2)interpolation of spatially and temporallysparse observations;(3) inference from available ob-asystemwithrelevancetoatleastsixoftheseareas:2005Blackwell PublishingLtd,Global ChangeBiology,11,378-397

Introduction Systematic earth observation implies the collection and interpretation of multiple kinds of data about the evolving state of the earth system across wide spatial domains and over extended time periods. Three factors have caused a massive acceleration in earth observation activities over recent years. The first is need: global change is raising issues – such as greenhouse-induced climate change, water shortages and imbalances, land degradation, soil erosion, loss of biodiversity – which require informed human responses at both global and regional levels. Second, technological advances in sensors, satellite systems and data storage and proces￾sing capabilities are making possible observations and interpretations which were out of reach only a few years ago and unimaginable a few decades ago. Third, the synthesis of formerly discrete disciplines into a unified Earth System Science is driving new hypotheses about the dynamics of the earth system and the interconnectedness of its components, including hu￾mans. Systematic earth observation motivates and tests these hypotheses. The focus of this paper is observation of the carbon cycle, and in particular its land-atmosphere compo￾nents, as one part of an integrated earth observation system. It is a significant part because of the coupling between the carbon cycle and the terrestrial cycles of water, energy and nutrients, and the connections of all these biospheric processes with global climate and human activities (Field & Raupach, 2004; Raupach et al., 2004). The carbon cycle is integral to the growth and decay of vegetation, maintains the water cycle through transpiration and provides habitat for maintaining biodiversity. Thus, terrestrial carbon observation is important for climate observation and prediction, for the management of water resources, nutrients and biodiversity, and for monitoring and managing the enhanced greenhouse effect. It is increasingly recognized that strategies for earth observation (including terrestrial carbon observation) require methods for combining data and process models in systematic ways. This is leading to research towards the application in terrestrial carbon observa￾tion (and in earth observation more generally) of ‘model–data synthesis’, the combination of the infor￾mation contained in both observations and models through both parameter-estimation and data-assimila￾tion techniques. Motivations for model–data synthesis approaches include (1) model testing and data quality control (through systematic checks for agreement within specified uncertainty bands for both data and model); (2) interpolation of spatially and temporally sparse observations; (3) inference from available ob￾servations of quantities which are not directly obser￾vable (such as carbon stores and fluxes over large areas) and (4) forecasting (prediction forward in time on the basis of past and current observations). The present paper arose from a workshop held in Sheffield, UK, 3–6 June 2003, to further the develop￾ment of a Terrestrial Carbon Observation System (TCOS) with a particular emphasis on model–data synthesis. Antecedents for this effort were (1) pre￾liminary steps toward a TCOS (Cihlar et al., 2002a, b, c); (2) a wider concept for an Integrated Global Carbon Observing Strategy including atmosphere, oceans, land and human activities (Ciais et al., 2004) and (3) the research program of the Global Carbon Project (Global Carbon Project, 2003). The paper is founded on three themes arising from the Sheffield workshop. First, model–data synthesis, based on terrestrial biosphere models constrained with multiple kinds of observation, is an essential compo￾nent of a TCOS. Second, from the standpoint of model– data synthesis, data uncertainties are as important as data values themselves and have a comparable role in determining the outcome. Third, and consequently, there is an urgent need for soundly based uncertainty specifications for the main kinds of data used in terrestrial carbon observation. These themes are devel￾oped as follows: the next section summarizes major purposes and attributes of a TCOS. ‘Model–data synthesis: methods’ provides an overview of model– data synthesis in the context of terrestrial carbon observation, by briefly describing some of the main methods, indicating their common characteristics, and highlighting the key role of data uncertainty. ‘Model– data synthesis: examples’ provides some examples. ‘Data characteristics: uncertainty in measurement and representation’ undertakes a survey of the uncertainty characteristics of the main kinds of relevant data. Purposes and attributes of a TCOS A succinct statement of the overall purpose of a TCOS might be: to operationally monitor the cycles of carbon and related entities (water, energy, nutrients) in the terrestrial biosphere, in support of comprehensive, sustained earth observation and prediction, and hence sustainable environmental management and socio￾economic development. These words are congruent with the Framework Document emerging from the Second Earth Observation Summit, Tokyo, April 2004 (http://earthobservations.org/docs/Framework%20- Doc%20Final.pdf), which calls for a ‘Global Earth Observation System of Systems’ to serve nine areas of socio-economic benefit. A TCOS is a contributor to such a system with relevance to at least six of these areas: MODEL –DATA SYNTHESIS IN TERRESTRIAL CARBON OBSERVATION 379 r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

380 M.R.RAUPACH etal.Model-data synthesis: methods. Understanding climate, and assessing and mitigat-ing climate change impacts;In this section, we survey a range of model-data.Improvingglobal waterresourcemanagementandsynthesis methods potentially applicablein a TCOS.understanding of the water cycle;More detail and further references can be found in a.Improving weatherinformation and prediction;growing number of excellent sources, for instanceTarantola (1987)and Evans &Stark (2002)for high-. Monitoring and managing inland ecosystems, in-level treatmentsof thegeneral statistical problemofcluding forests, and land use change;inverseestimation,Grewal&Andrews(1993)and. Supporting sustainable agriculture and combatingDrecourt (2003)for introductions to theKalmanFilter,desertification;Reichle et al.(2002)for hydrological applications with.Understanding,monitoring and preventing loss ofanemphasisontheKalmanFilterandEnting(2002)biodiversity.andKasibhatla etal. (2000)for applications of a range ofmethodstobiogeochemical cycles.Tomakethesecontributionseffectively,a TCOSmustOverviewhave a number of attributes (see also Running et al.,1999; Cihlar et al.,2002a; Ciais et al.,2004).First, scientificThe central problem is: using appropriate observationscredibility is needed to maintain methodological andand models,we must determine the spatial distribu-observational rigour,and to include procedures fortions andtemporalevolutionsof theterrestrialstoresestimating uncertainties or confidence limits.Second,and fluxes of carbon and related entities (waterconsistency withglobal budgets is necessaryto respectnutrients, energy)across the earth.Important fluxesconstraints from global-scale carbon and related bud-includeland-air exchanges (atmospheric sources andgets incorporating terrestrial, atmospheric and oceanicsinks), exchanges with rivers and groundwater, andpools and anthropogenic sources such as fossil fuelexchanges betweenterrestrial pools such asbiomassburning.Third, sufficient spatial resolution is necessaryand soil. We also need to determine the main processesto resolve spatial variations in patterns of land useinfluencing the fluxes,including those under human(typicallytens of metres, consistent with high-resolu-management.No singlemodel or set of observationstion remote sensing).Fourth,enough temporal resolutioncan supplythis amountofinformation-hencethe needis needed to resolve the influence of weather, inter-for a synthesis approach.The task of combiningannual climate fluctuations and long-term climateobservations and modelscanbecarriedoutinmanychange on carbon and related cycles. Fifth, the systemways,encompassed by the umbrellaterms'model-dataneedstoencompassa broad range ofentities,eventuallysynthesis'or'model-data fusion'.Thegeneral principleincluding CO2, CH4, CO,volatile organic carbonsis to find an 'optimal match'between observations and(VOCs)and aerosol black carbon.Ofthese,thehighestmodel by varying one or more 'properties' of thepriority is CO.Water is also a high priority because ofmodel.(Words in quotes have specific meaningsits importance in modulating other terrestrial GHGdefined below). The optimal match is a choice of modelfluxes. Sixth, a sufficient range of processes must beproperties, which minimizes the'distance' between theencompassed.Ahighpriorityis resolution of netland-model representations of a system and what weknowair fluxes of greenhousegases in which all terrestrialabout thereal biophysical system from observationalsources and sinks are lumped together. However, thereand prior'data'.At this high level of generality,model-is an equally high demand for identification of thedata synthesis encompasses both'parameter estima-terms contributing to thenetfluxes,for exampletotion'and 'data assimilation'.All applications rest onpartition a net flux between vegetation and soil storagethree foundations: a model of the system, data aboutchanges.Finally,quantification of uncertainty is required.the system, and a synthesis approach.The'demand side' of the uncertainty issue is:whatlevel of uncertainty is acceptablefora TCOS to offerModel.For a TCOS, the model is a terrestrial biosphereuseful information?The answer is not simpleandmodel describing the evolving stores and fluxes ofdepends on theapplication,forexample,fromtheareascarbon, water, energy and related entities. This dynamicmentioned above.This paper does not attempt tomodel hastheformanswer the demand-side question, but rather concen-dxtrates on the'supplvside'of uncertaintv:that is,how=f(x, u,p) + noiseordtuncertainty can be determined in a TCOS based onx"+1 =p(x",u", p) + noise =x"+At f(x",u",p)+ noise,model-datasynthesisandmultipleobservation(1)sources, each with its own specified uncertainty.2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

Understanding climate, and assessing and mitigat￾ing climate change impacts;
Improving global water resource management and understanding of the water cycle;
Improving weather information and prediction;
Monitoring and managing inland ecosystems, in￾cluding forests, and land use change;
Supporting sustainable agriculture and combating desertification;
Understanding, monitoring and preventing loss of biodiversity. To make these contributions effectively, a TCOS must have a number of attributes (see also Running et al., 1999; Cihlar et al., 2002a; Ciais et al., 2004). First, scientific credibility is needed to maintain methodological and observational rigour, and to include procedures for estimating uncertainties or confidence limits. Second, consistency with global budgets is necessary to respect constraints from global-scale carbon and related bud￾gets incorporating terrestrial, atmospheric and oceanic pools and anthropogenic sources such as fossil fuel burning. Third, sufficient spatial resolution is necessary to resolve spatial variations in patterns of land use (typically tens of metres, consistent with high-resolu￾tion remote sensing). Fourth, enough temporal resolution is needed to resolve the influence of weather, inter￾annual climate fluctuations and long-term climate change on carbon and related cycles. Fifth, the system needs to encompass a broad range of entities, eventually including CO2, CH4, CO, volatile organic carbons (VOCs) and aerosol black carbon. Of these, the highest priority is CO2. Water is also a high priority because of its importance in modulating other terrestrial GHG fluxes. Sixth, a sufficient range of processes must be encompassed. A high priority is resolution of net land￾air fluxes of greenhouse gases in which all terrestrial sources and sinks are lumped together. However, there is an equally high demand for identification of the terms contributing to the net fluxes, for example to partition a net flux between vegetation and soil storage changes. Finally, quantification of uncertainty is required. The ‘demand side’ of the uncertainty issue is: what level of uncertainty is acceptable for a TCOS to offer useful information? The answer is not simple and depends on the application, for example, from the areas mentioned above. This paper does not attempt to answer the demand-side question, but rather concen￾trates on the ‘supply side’ of uncertainty: that is, how uncertainty can be determined in a TCOS based on model–data synthesis and multiple observation sources, each with its own specified uncertainty. Model–data synthesis: methods In this section, we survey a range of model–data synthesis methods potentially applicable in a TCOS. More detail and further references can be found in a growing number of excellent sources, for instance Tarantola (1987) and Evans & Stark (2002) for high￾level treatments of the general statistical problem of inverse estimation, Grewal & Andrews (1993) and Dre´court (2003) for introductions to the Kalman Filter, Reichle et al. (2002) for hydrological applications with an emphasis on the Kalman Filter and Enting (2002) and Kasibhatla et al. (2000) for applications of a range of methods to biogeochemical cycles. Overview The central problem is: using appropriate observations and models, we must determine the spatial distribu￾tions and temporal evolutions of the terrestrial stores and fluxes of carbon and related entities (water, nutrients, energy) across the earth. Important fluxes include land–air exchanges (atmospheric sources and sinks), exchanges with rivers and groundwater, and exchanges between terrestrial pools such as biomass and soil. We also need to determine the main processes influencing the fluxes, including those under human management. No single model or set of observations can supply this amount of information – hence the need for a synthesis approach. The task of combining observations and models can be carried out in many ways, encompassed by the umbrella terms ‘model–data synthesis’ or ‘model–data fusion’. The general principle is to find an ‘optimal match’ between observations and model by varying one or more ‘properties’ of the model. (Words in quotes have specific meanings defined below). The optimal match is a choice of model properties, which minimizes the ‘distance’ between the model representations of a system and what we know about the real biophysical system from observational and prior ‘data’. At this high level of generality, model– data synthesis encompasses both ‘parameter estima￾tion’ and ‘data assimilation’. All applications rest on three foundations: a model of the system, data about the system, and a synthesis approach. Model. For a TCOS, the model is a terrestrial biosphere model describing the evolving stores and fluxes of carbon, water, energy and related entities. This dynamic model has the form dx dt ¼fðx; u; pÞ þ noise or xnþ1 ¼uðxn; un; pÞ þ noise ¼ xn þ Dtfðxn; un; pÞ þ noise; ð1Þ 380 M. R. RAUPACH et al. r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

MODEL-DATA SYNTHESISINTERRESTRIALCARBONOBSERVATION381where x is a vector of state wariables (such as stores ofhydrological data on river flows, groundwater, andcarbon, water and related entities, or store attributesconcentrations of C,N and other entities;(6) soilsuch as age class distributions); f is a vector of rates ofpropertiesand topography;(6)disturbancerecordschange (net fluxes where components of x are stores); (both contemporary and historical) including landis the discrete analogue for f; u is a set of externallymanagement, land use, land use change and fire andspecified time-dependent forcing variables (such as(8)climate and weather data (precipitation, solarmeteorological variables and soil properties) and p is aradiation,temperature and humidity).Of these, someset of time-independent model parameters (such as rate(especiallythefirstfive)typicallyprovideobservationalconstantsandpartitionratios).Inthediscreteconstraints (z), while others provide model drivers (u).formulation, time steps are denoted by superscripts.Examples of observation models(Eqn (2))includeThenoiseterms account for both imperfectionsinradiative transfer models to map modelled surfacemodel formulation and stochastic variability in forcingsstatesintotheradiancesobserved bysatellites;(u)or parameters (p).Once themodelfunction f(x,u,p)atmospherictransportmodelstotransformmodelledor (x",u",p) is specified, then the system evolution x(t)surfacefluxes to measured atmospheric concentrations;can be determined by integrating Eqn (1) in time (withandallometricrelationstotransformmodelledbiomasszero noise), from initial conditions x(0), with specifiedto observed treediameters.external forcing u(t) and parameters p.Symthesis.The final requirement is a synthesis process,ora systematicmethod forfindingtheoptimal matchData. These are generally of two broad kinds: (1)between the data (including observations and priorobservations or measurements of a set of quantitieszestimates)and themodel.Thisprocess needstoprovideand (2)priorestimates formodel quantities (x,u and p).three kinds of output: optimal estimates for the modelBoth include uncertainty,through errors and noise. Inproperties to be adjusted, uncertainty statements aboutthis paper, the term'data'includes both observationsthese estimates, and an assessment of how well theand prior estimates, and incorporates the uncertaintymodel fits the data, given the data uncertainties. In anyinherent in each.synthesis process, there are three basic choices: (1)theThe measured quantities (z) are related to the systemmodel properties to be adjusted or 'target variables', (2)state and external forcing variables by an obseroationthe measure of distance between data and model ormodel of the form'cost function'and (3) the search strategyfor finding the(2)z= h(x,u) +noise,optimum values.Search strategies can be classifiedwhere the operator h specifies the deterministicbroadly into (3a)‘nonsequential'or“batch'strategies inrelationship between the measured quantities and thewhich the data are treated all at once, and (3b)system state.The noise term accounts for both'sequential' strategies in which the data arrive in a'measurement error(instrumental and processingtime sequence and are incorporated into the model-errors in the measurements z), and 'representationdata synthesis step by step. The rest of this sectionerror'(errors in the model representation of z,explores the choices (1),(2),(3a)and (3b).introduced by shortcomings in the observation modelh). In the rare case where we can observe all stateTarget variablesvariables directly,h reduces to the identity operator, soThe target variables are the properties of the model toz=x+ (measurement) noise. In time-discrete form, Eqn(2)becomes z"=h(x"u")+noise.Note the inter-be adjusted in the optimization process.They includepretation of the time-step superscripts:x"and u" areany model property considered to be sufficientlysimply the model state and forcings at time step n,uncertain as to benefit from constraint by the data.whereas z" is the set of new observations introduced atModel properties which can be target variables include:time step n, whatever the actual time of its measu-(1) model parameters (p); (2) forcing variables (u"), ifthere is substantial uncertainty about them; (3) initialrement.However, no observations may be used moreconditions on the state variables (x) and (4)time-than once.Examples of potential observations in a TCOSdependent components of the state vector x".Theinclude(1)atmospheric composition (concentrationsinclusion of the state vector x" as a possible targetof CO2 and other gases); (2) remote sensing of terrestrialvariable is for the following reason: in a purelyandatmosphericproperties;(3)fluxesof carbon anddeterministic model the trajectory x" is determined byrelated entities,with supportingprocessobservations,the dynamical model (f or ), the values of p and u",and the initial value x It might seem sufficient,at intensive study sites; (4) vegetation and soil stores ofcarbon from forest and ecological inventories; (5)therefore, to estimate these and allow integration of2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

where x is a vector of state variables (such as stores of carbon, water and related entities, or store attributes such as age class distributions); f is a vector of rates of change (net fluxes where components of x are stores); u is the discrete analogue for f; u is a set of externally specified time-dependent forcing variables (such as meteorological variables and soil properties) and p is a set of time-independent model parameters (such as rate constants and partition ratios). In the discrete formulation, time steps are denoted by superscripts. The noise terms account for both imperfections in model formulation and stochastic variability in forcings (u) or parameters (p). Once the model function f(x, u, p) or u(xn , un , p) is specified, then the system evolution x(t) can be determined by integrating Eqn (1) in time (with zero noise), from initial conditions x(0), with specified external forcing u(t) and parameters p. Data. These are generally of two broad kinds: (1) observations or measurements of a set of quantities z and (2) prior estimates for model quantities (x, u and p). Both include uncertainty, through errors and noise. In this paper, the term ‘data’ includes both observations and prior estimates, and incorporates the uncertainty inherent in each. The measured quantities (z) are related to the system state and external forcing variables by an observation model of the form z ¼ hðx; uÞ þ noise; ð2Þ where the operator h specifies the deterministic relationship between the measured quantities and the system state. The noise term accounts for both ‘measurement error’ (instrumental and processing errors in the measurements z), and ‘representation error’ (errors in the model representation of z, introduced by shortcomings in the observation model h). In the rare case where we can observe all state variables directly, h reduces to the identity operator, so z 5 x 1 (measurement) noise. In time-discrete form, Eqn (2) becomes zn 5 h(xn , un ) 1 noise. Note the inter￾pretation of the time-step superscripts: xn and un are simply the model state and forcings at time step n, whereas zn is the set of new observations introduced at time step n, whatever the actual time of its measu￾rement. However, no observations may be used more than once. Examples of potential observations in a TCOS include (1) atmospheric composition (concentrations of CO2 and other gases); (2) remote sensing of terrestrial and atmospheric properties; (3) fluxes of carbon and related entities, with supporting process observations, at intensive study sites; (4) vegetation and soil stores of carbon from forest and ecological inventories; (5) hydrological data on river flows, groundwater, and concentrations of C, N and other entities; (6) soil properties and topography; (6) disturbance records (both contemporary and historical) including land management, land use, land use change and fire and (8) climate and weather data (precipitation, solar radiation, temperature and humidity). Of these, some (especially the first five) typically provide observational constraints (z), while others provide model drivers (u). Examples of observation models (Eqn (2)) include radiative transfer models to map modelled surface states into the radiances observed by satellites; atmospheric transport models to transform modelled surface fluxes to measured atmospheric concentrations; and allometric relations to transform modelled biomass to observed tree diameters. Synthesis. The final requirement is a synthesis process, or a systematic method for finding the optimal match between the data (including observations and prior estimates) and the model. This process needs to provide three kinds of output: optimal estimates for the model properties to be adjusted, uncertainty statements about these estimates, and an assessment of how well the model fits the data, given the data uncertainties. In any synthesis process, there are three basic choices: (1) the model properties to be adjusted or ‘target variables’, (2) the measure of distance between data and model or ‘cost function’ and (3) the search strategy for finding the optimum values. Search strategies can be classified broadly into (3a) ‘nonsequential’ or ‘batch’ strategies in which the data are treated all at once, and (3b) ‘sequential’ strategies in which the data arrive in a time sequence and are incorporated into the model– data synthesis step by step. The rest of this section explores the choices (1), (2), (3a) and (3b). Target variables The target variables are the properties of the model to be adjusted in the optimization process. They include any model property considered to be sufficiently uncertain as to benefit from constraint by the data. Model properties which can be target variables include: (1) model parameters (p); (2) forcing variables (un ), if there is substantial uncertainty about them; (3) initial conditions on the state variables (x0 ) and (4) time￾dependent components of the state vector xn . The inclusion of the state vector xn as a possible target variable is for the following reason: in a purely deterministic model the trajectory xn is determined by the dynamical model (f or u), the values of p and un , and the initial value x0 . It might seem sufficient, therefore, to estimate these and allow integration of MODEL –DATA SYNTHESIS IN TERRESTRIAL CARBON OBSERVATION 381 r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

382 M. R. RAUPACH et al.the model to take care of x".However, the model itselfEquation (3) defines the generalized least squares costfunction minimized by theminimum-variance estimatemay not be perfect, as indicated by the noise term inEqn (1), so there may be advantage in adjusting values(y)fory.Forany distribution of the errors in the dataof x" through the model integration.(observations z and priors y),this estimate is unbiased,To maintain generality,wedenote thevector of targetand hasthe minimum error covariance among all linear(in z), unbiased estimates (Tarantola 1987). Use of Eqnvariables by y.This vectormay or may not be a functionof time, and will usuallybea subset of all model(3)has another,additional foundation:provided thatthe probability distributions for data errors are Gaus-variables (x",u",p). Broadly speaking, parameter esti-sian, it yields a maximum-likelihood estimate for y,mation problems are those where the target variablesare restricted to model parameters (p), while dataconditional on the data and the model dynamics (Presset al.,1992,p.652;Todling2000).Outside the restrictionassimilationproblemsmayincludeanymodelpropertyas a target variable,usually with an emphasis on stateof Gaussian distributions,y as defined by minimizingavariables (x").quadratic Jis not exactly the maximum-likelihoodestimate,but it is often not far from it.AquadraticJis widely used even when the data errors are notCostfunctionGaussian; see Press et al. (1992, p. 690) for discussion.There are alternative cost functions Jin whichmodel-The cost or objective function J (a function of the targetmeasurementdifferences(z-h(y))areraisedtopowersvariablesy)definesthemismatchor distancebetweenother than 2,the choicein Eqn(3)(Tarantola,1987;the model and the data. It can take a wide range ofGershenfeld, 1999).For example, in flood event model-forms,but musthave certain properties (for example, itling,theabsolutemaximumerrorisneededtocapturemust bemonotonic in the absolutedifferencebetweenpeak flowrates,whilefor modelling baseflow rates, thedata and model-predicted values).Acommon choice ismean absolute deviation (Iz-h(y)I to the power 1)hasthe quadratic cost function:the desirable property of being less sensitive to outliersJ(y) =(z - h(y)[Covz)-'(z - h(y)than a power 2. Different powers for Iz-h(y)I produce(3)maximum-likelihood estimates for y with different+ (y -y)"[Covy]-'(y-),distributionsfordata errors; for example,a power 1Jwhere y is the vector of 'priors' (a priori estimates) foryields a maximum-likelihood estimate when the datathe target variables, and [Covz] and [Covy] areerrors are distributed exponentially,and a high-powerJcovariancematricesforzand y,respectivelypreferentially weights outliers such as peak flows.(Covz)mm=(2mzh),with Zm= zm-(zm),angle bracketsHere, we use a power 2 Jexclusively.denoting the expectation operator).The first term inEqn (3) is a sum of the squared distances betweenSearch strategies for nonsequential problemsmeasured components oftheobservationvector(z)andtheir model predictions (h(y), while the second is aIn nonsequential or batch problems,all data are treatedcorresponding sum of distances between target vari-simultaneouslyandtheminimizationproblem is solvedables and their prior estimates. The matrices [Covz]-1only once. A familiar case is least-squares parameterand [Covy]-'represent the weights accorded to theestimation.observations and the priors, and thus scale theExample.Someoftheattributesoftheseproblemsareconfidencesaccordedtoeach.Theirrolecanbeclarifiedby considering the simple casein which components zmdemonstrated by considering a simple linear exampleof the observation vector z are independent, withwhichextendstheparameter-estimation:problem.variances om; then [Covz]-1 is the diagonal matrixAlthoughmathematically straightforward,this casediag [1/o2] and the squared departures of the measure-finds important application in the atmosphericments(zm)fromthepredictions (hm(y)areseentobeinversion methods used to estimate trace gas sourcesweighted by the confidence measure 1/for eachfromatmosphericcompositionobservations(seeModel-data synthesis: Examples'). Here the targetcomponent.Themodel-data synthesis problem now becomes:variables (y) are a set of surface-air fluxes, averagedvaryy to minimize J(y), subject to the constraint thatover suitable areas; there is no dynamic model relatingx(t) must satisfy the dynamic model, Eqn (1). The valuefluxes at different times and places to each other; andof y at the minimum is the a posteriori estimate of y,theobservationoperator(h)isamodelofatmosphericincluding information from the observations as well astransport. From the linearity of the conservationthe priors. We denote it by y (so frowns and smilesequation for an inert trace gas, it follows that h isrespectivelydesignate priorand posterior estimates)linear and can hence be represented by a matrix H2005Blackwell Publishing Ltd, Global Change Biology,11,378-397

the model to take care of xn . However, the model itself may not be perfect, as indicated by the noise term in Eqn (1), so there may be advantage in adjusting values of xn through the model integration. To maintain generality, we denote the vector of target variables by y. This vector may or may not be a function of time, and will usually be a subset of all model variables (xn , un , p). Broadly speaking, parameter esti￾mation problems are those where the target variables are restricted to model parameters (p), while data assimilation problems may include any model property as a target variable, usually with an emphasis on state variables (xn ). Cost function The cost or objective function J (a function of the target variables y) defines the mismatch or distance between the model and the data. It can take a wide range of forms, but must have certain properties (for example, it must be monotonic in the absolute difference between data and model-predicted values). A common choice is the quadratic cost function: JðyÞ ¼ðz hðyÞÞT½Cov z 1 ðz hðyÞÞ þ ðy y _Þ T½Cov y _ 1 ðyy _Þ; ð3Þ where y _ is the vector of ‘priors’ (a priori estimates) for the target variables, and [Cov z] and ½Cov y _ are covariance matrices for z and y _, respectively (½Cov z mn ¼ z0 mz0 n
, with z0 m ¼ zm h i zm , angle brackets denoting the expectation operator). The first term in Eqn (3) is a sum of the squared distances between measured components of the observation vector (z) and their model predictions (h(y)), while the second is a corresponding sum of distances between target vari￾ables and their prior estimates. The matrices [Cov z] 1 and ½Cov y _ 1 represent the weights accorded to the observations and the priors, and thus scale the confidences accorded to each. Their role can be clarified by considering the simple case in which components zm of the observation vector z are independent, with variances s2 m; then [Cov z] 1 is the diagonal matrix diag ½1=s2 m and the squared departures of the measure￾ments (zm) from the predictions (hm(y)) are seen to be weighted by the confidence measure 1=s2 m for each component. The model–data synthesis problem now becomes: vary y to minimize J(y), subject to the constraint that x(t) must satisfy the dynamic model, Eqn (1). The value of y at the minimum is the a posteriori estimate of y, including information from the observations as well as the priors. We denote it by y ^ (so frowns and smiles respectively designate prior and posterior estimates). Equation (3) defines the generalized least squares cost function minimized by the minimum-variance estimate (y ^) for y. For any distribution of the errors in the data (observations z and priors y _), this estimate is unbiased, and has the minimum error covariance among all linear (in z), unbiased estimates (Tarantola 1987). Use of Eqn (3) has another, additional foundation: provided that the probability distributions for data errors are Gaus￾sian, it yields a maximum-likelihood estimate for y, conditional on the data and the model dynamics (Press et al., 1992, p. 652; Todling 2000). Outside the restriction of Gaussian distributions, y ^ as defined by minimizing a quadratic J is not exactly the maximum-likelihood estimate, but it is often not far from it. A quadratic J is widely used even when the data errors are not Gaussian; see Press et al. (1992, p. 690) for discussion. There are alternative cost functions J in which model￾measurement differences (zh(y)) are raised to powers other than 2, the choice in Eqn (3) (Tarantola, 1987; Gershenfeld, 1999). For example, in flood event model￾ling, the absolute maximum error is needed to capture peak flow rates, while for modelling base flow rates, the mean absolute deviation (|z–h(y)| to the power 1) has the desirable property of being less sensitive to outliers than a power 2. Different powers for |z–h(y)| produce maximum-likelihood estimates for y ^ with different distributions for data errors; for example, a power 1 J yields a maximum-likelihood estimate when the data errors are distributed exponentially, and a high-power J preferentially weights outliers such as peak flows. Here, we use a power 2 J exclusively. Search strategies for nonsequential problems In nonsequential or batch problems, all data are treated simultaneously and the minimization problem is solved only once. A familiar case is least-squares parameter estimation. Example. Some of the attributes of these problems are demonstrated by considering a simple linear example, which extends the parameter-estimation problem. Although mathematically straightforward, this case finds important application in the atmospheric inversion methods used to estimate trace gas sources from atmospheric composition observations (see ‘Model–data synthesis: Examples’). Here the target variables (y) are a set of surface-air fluxes, averaged over suitable areas; there is no dynamic model relating fluxes at different times and places to each other; and the observation operator (h) is a model of atmospheric transport. From the linearity of the conservation equation for an inert trace gas, it follows that h is linear and can hence be represented by a matrix H 382 M. R. RAUPACH et al. r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

MODEL-DATASYNTHESISINTERRESTRIALCARBONOBSERVATION383(Raupach, 2001), thus, z=Hy+noise.For now, thedata and the uncertainties are completely inseparable intheformalism.Toput thepoint provocatively,noiseisassumedtobeGaussianwithzeromeanandnotemporal correlation,and thuscompletely characterisedproviding data and allowing another researcher toby an observation error covariancematrix [Covz].Byprovide the uncertainty is indistinguishable fromminimizing Janalytically,one obtains the expressionallowing the second researcher to make up the data in(Tarantola,1987,p.196; Enting,2002):the first place.This realization informs the emphasis onuncertaintythroughout this paper.=y+[CovH'[Covz]-'(z-Hy),(4)where [Cov y], the estimated error covariance of the aposterioriestimatey,isgiven byAlgorithmsfornonsequential problems[Cov y]-1 = [Cov y]-1 + H'[Covz]-"H.(5)The task in general is to find the target variables ywhich minimize J(y). Clearly, the shape of J(y) is allThese expressions already tell us some importantimportant:itmay havea singleminimum or multiplethings. The posterior estimates are given by the priorseparated local minima, only one of which is the trueestimates plus a term depending on themismatchglobal minimum.Neartheminimum,Jmaybe shapedbetween the experimental observations and thelike a long,narrowellipsoidal valley.If this valley hasaobservations as predicted bythe prior estimates.Thisflatfloortracing out some line iny space,thenall pointsmismatch is weighted by our confidence in thealong that line are equally acceptable and these yobservations,[Covz]-1.Thus,observations with littlecoordinates cannot be distinguished in terms ofweight hardly shiftthe posteriorestimate from theoptimality,so such combinations of target variablesprior, and vice versa.Furthermore, the transpose of thecannot be resolved by model-data synthesis with theobservation operator (H) multiplies the weightedavailable data and model.Diagnostic indicators aboutmismatch. If this operator is very weak, that is if thethese issues are provided by the Hessian or curvatureavailableobservations are only weakly related tothematrix D='l/Qyjoyx,a measure of the local curvaturetargetvariables,then the update to the initial estimate isofJ(y).Thedegreeof orthogonalityamongcolumnsofalso small.Finallythe posterior covariance [Cov y] (EqnD indicates the extent to which it is possible to find a(5)) is bounded above,in some sense,by the priorunique local minimumtoJ(y)inthevicinityofthepointcovariance [Covy]. If the prior covariance is smallat which D is evaluated. A high ‘condition number'(suggestingsubstantial confidence in theinitial(ratioof largestto smallest eigenvalue)forDindicatesestimate)then the increment y-y (the differencethat some linear combination(s)of the columns of Darebetween the posteriorand prior estimates, a measurenearly zero, that is, that the curvature is nearly zero inof the information added by the observations z, andsomedirection(s),sothattheminimizationproblemisequal to the second term in Eqn (4) in the present case)ill-conditioned, as in the case of a valley witha flat floor.is also small.Giventheseconsiderations,classes ofmethod forAll the above is reasonable. More surprising is thefinding theminimum in J(y)include thefollowingrelationship between the data, its uncertainty and the1.Analytic solution is possiblewhenthe observationcost function.We candecompose the (positivedefinite)operator h(y) is linear (z=Hy +noise). In this case I(y)matrix [Covz]-1 into a matrix product A'A, using theis a quadratic form shaped like a parabolic bowl, andCholeskydecompositionforapositivedefinitematrix.theminimization can be carried out analytically as inFor a diagonal covariance matrix, diag [], thethe example of Eqns (4) and (5). This 'direct or one-decomposition is trivial:A=diag [1/2].Likewise, westep' solution is highly efficient when applicable;can write [Cov y]-1 =B'B. The cost function, Eqn (3),however,mostproblemsarenonlinearand requireacanthenberewrittenasnonlinearmethod.J(y) = (a - Ah(y)*(a - Ah(y)2. Gradient descent algorithms arethe most familiar+(b - b)(b - b)search algorithms for nonlinear optimization. They(6)include (for example) steepest-descent, conjugate-gra-where a=Az, b=By and b=By. Thus the costdient, quasi-Newton and Levenberg-Marquardt algo-function, and thence the entire minimization,takes arithms (Press et al.,1992).Gradient-descent methods areformin which neither theobservations nor the prioreasily implemented, provided that the gradient vectorestimatesappear; they are replaced byquantities a andVJ=OJ/Oykcanbe calculated.Themain advantages ofb scaled by the square roots of the inverse covariancegradient-descent algorithms are relative simplicity andmatrices, which are measures of confidence.This is nolowcost;themaindisadvantageisthatifthesurfacemathematical nicety;rather it demonstrates that theJ(y)hasmultipleminima,theytendtofind local2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

(Raupach, 2001), thus, z 5 Hy 1 noise. For now, the noise is assumed to be Gaussian with zero mean and no temporal correlation, and thus completely characterised by an observation error covariance matrix [Cov z]. By minimizing J analytically, one obtains the expression (Tarantola, 1987, p. 196; Enting, 2002): y ^ ¼ y _ þ Cov y ^ h iHT½ Cov z 1 z Hy_ ; ð4Þ where ½Cov y ^, the estimated error covariance of the a posteriori estimate y ^, is given by ½Cov y ^ 1 ¼ ½Cov y _ 1 þ HT½Cov z 1 H: ð5Þ These expressions already tell us some important things. The posterior estimates are given by the prior estimates plus a term depending on the mismatch between the experimental observations and the observations as predicted by the prior estimates. This mismatch is weighted by our confidence in the observations, [Cov z] 1 . Thus, observations with little weight hardly shift the posterior estimate from the prior, and vice versa. Furthermore, the transpose of the observation operator (HT ) multiplies the weighted mismatch. If this operator is very weak, that is if the available observations are only weakly related to the target variables, then the update to the initial estimate is also small. Finally the posterior covariance ½Cov y ^ (Eqn (5)) is bounded above, in some sense, by the prior covariance ½Cov y _. If the prior covariance is small (suggesting substantial confidence in the initial estimate) then the increment y ^ y _ (the difference between the posterior and prior estimates, a measure of the information added by the observations z, and equal to the second term in Eqn (4) in the present case) is also small. All the above is reasonable. More surprising is the relationship between the data, its uncertainty and the cost function. We can decompose the (positive definite) matrix [Cov z] 1 into a matrix product AT A, using the Cholesky decomposition for a positive definite matrix. For a diagonal covariance matrix, diag ½s2 m, the decomposition is trivial: A 5 diag ½1=s2 m. Likewise, we can write ½Cov y _ 1 5 BT B. The cost function, Eqn (3), can then be rewritten as JðyÞ¼ða AhðyÞÞTða AhðyÞÞ þðb b _ Þ Tðb b _ Þ; ð6Þ where a 5 Az, b 5 By and b _ ¼ By_. Thus the cost function, and thence the entire minimization, takes a form in which neither the observations nor the prior estimates appear; they are replaced by quantities a and b scaled by the square roots of the inverse covariance matrices, which are measures of confidence. This is no mathematical nicety; rather it demonstrates that the data and the uncertainties are completely inseparable in the formalism. To put the point provocatively, providing data and allowing another researcher to provide the uncertainty is indistinguishable from allowing the second researcher to make up the data in the first place. This realization informs the emphasis on uncertainty throughout this paper. Algorithms for nonsequential problems The task in general is to find the target variables y which minimize J(y). Clearly, the shape of J(y) is all important: it may have a single minimum or multiple separated local minima, only one of which is the true global minimum. Near the minimum, J may be shaped like a long, narrow ellipsoidal valley. If this valley has a flat floor tracing out some line in y space, then all points along that line are equally acceptable and these y coordinates cannot be distinguished in terms of optimality, so such combinations of target variables cannot be resolved by model–data synthesis with the available data and model. Diagnostic indicators about these issues are provided by the Hessian or curvature matrix D 5 @2 J/@yj@yk, a measure of the local curvature of J(y). The degree of orthogonality among columns of D indicates the extent to which it is possible to find a unique local minimum to J(y) in the vicinity of the point at which D is evaluated. A high ‘condition number’ (ratio of largest to smallest eigenvalue) for D indicates that some linear combination(s) of the columns of D are nearly zero, that is, that the curvature is nearly zero in some direction(s), so that the minimization problem is ill-conditioned, as in the case of a valley with a flat floor. Given these considerations, classes of method for finding the minimum in J(y) include the following. 1. Analytic solution is possible when the observation operator h(y) is linear (z 5 Hy 1 noise). In this case J(y) is a quadratic form shaped like a parabolic bowl, and the minimization can be carried out analytically as in the example of Eqns (4) and (5). This ‘direct’ or ‘one￾step’ solution is highly efficient when applicable; however, most problems are nonlinear and require a nonlinear method. 2. Gradient descent algorithms are the most familiar search algorithms for nonlinear optimization. They include (for example) steepest-descent, conjugate-gra￾dient, quasi-Newton and Levenberg–Marquardt algo￾rithms (Press et al., 1992). Gradient-descent methods are easily implemented, provided that the gradient vector HyJ 5 @J/@yk can be calculated. The main advantages of gradient-descent algorithms are relative simplicity and low cost; the main disadvantage is that if the surface J(y) has multiple minima, they tend to find local MODEL –DATA SYNTHESIS IN TERRESTRIAL CARBON OBSERVATION 383 r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

384M.R.RAUPACHetal.StateStateminima near thestartingvalue of yrather than theObservablevariableuncertaintyglobal minimum.3.Global searchmethods find theglobal minimum inaTime step n-ja-1[Cov yn-]2~1(posterior)function J(y) by searching (effectively) the whole of y1olo.oPredictionspace.They overcome the local-minimum pitfall (so tostepTimestepnspeak)of gradient-descent methods, but have the2Hyn[Covy"(prior)disadvantage of higher computational costs. SimulatedH,Rannealing and genetic algorithms are two examples.AnalysisdifferenThesemethodsareefficient atfindingthevicinityof asteps个global minimumwherethere may be multiplelocal-Time step nzynCovynminima, but do not locate an exact local minimum.(posterior)They may be combined withgradient-descent methodsFig.1 Information flow in the linear Kalman filter. Linearfor finding an exact global minimum once in the rightoperators are indicated nexttoarrows.Operationsinpredictionvicinity.and analysis steps are shown as dashed and solid lines,respectively.Search strategies for sequential problemsIn sequential problems, the task is to solve fora set ofvariablesy from one time step tothe next, by atarget variables y"associated with a particular timelinearized version of the model.The second term (Q) isstep, usually including the state variables of thethecovarianceofthenoiseterminthedvnamicmodeldynamic model (x").Theprocess is then repeatedEqn (1),which includes both model imperfections andsequentially to give a time history fory.Informationstochastic variability in forcings and parameters. Thisabout y" can come from two sources:evolution of theterm plays a crucial role in the Kalman filter: itdynamic model from the previous time step, andquantifies our lack of confidence in the ability of thecomparison between the observations at the currentdynamic model to propagate the model state, and istime step (z")and the model predictions (h(y").usually referred to as model error. In most imple-mentations of theKalman filter themodel errorisKalmanfilter.IntroducedbyKalman (1960),theKalmanassumed to be Gaussian with zero mean and notem-filter is by nowa group of algorithmsfor the sequentialporal correlation, and thus completely characterized bycombinationofdynamicandobservationalthecovariancematrixQ.information, using a ‘prediction' step and anIn the analysis step,the prior estimates are refined'analysis'step.In the prediction step,the dynamicby the inclusion of data. This is done using the priormodel is used to calculate prior estimates y for theestimateforthepredictedobservationvector,target variables at time step n,from the best (posterior)2"=h(y"), and its covarianceestimatesyatthepreviousstep.Intheanalysisstep[Cov 2"] = H[Cov y"JHT + R,(8)posterior estimates y at step n are obtained by"improvingtheprior estimateswith data.Themodelwhere H=Oh/oyistheJacobian matrix of thestate is then ready for evolution to the next (n+1) timeobservation model h(y), and R is the data covariancestep.Akeypoint isthatthe confidence in the currentmatrix [Covz], indicating lack of confidence in the datastate, embodied in the error covariance for the targetand often called the data error. Again it is usuallyvariables y,isalso evolved withthe dynamicmodel andassumed that thedata erroris Gaussianwith zeromeanimproved with observations.A schematic diagram ofand notemporal correlation,and thus completelythe information flow in the Kalman filter is given in Fig.characterized by R=[Covz].1.Theexpressionsforthefinal (posterior)estimatesforIntheprediction step,thetaskof evolving yisyand itscovariancearenowexactlyasforthecommon acrossall implementations of the Kalmannonsequential mode, except that the operation isfilter since it involves only a normal forward step of thecarried out for onetimestep only:dynamic model:y"=(y").Theprior estimatefor" -y" +[Covy"jH[Cov2"-(" - h(")the covariance at time step n evolves according to(9)=y" +K(z" - h(y"),n-ljoT +Q,[Cov y"] = @[Cov y"(7)where@=p/oy,theJacobianmatrix of thedynamic[Covy"]-} =[Cov y"]-" + H'RH or(10)model p(y).The first term on the right represents the[Cov y"] =(I - KH)[Cov y' ],propagation of the error covariance in the target2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

minima near the starting value of y rather than the global minimum. 3. Global search methods find the global minimum in a function J(y) by searching (effectively) the whole of y space. They overcome the local-minimum pitfall (so to speak) of gradient-descent methods, but have the disadvantage of higher computational costs. Simulated annealing and genetic algorithms are two examples. These methods are efficient at finding the vicinity of a global minimum where there may be multiple local minima, but do not locate an exact local minimum. They may be combined with gradient-descent methods for finding an exact global minimum once in the right vicinity. Search strategies for sequential problems In sequential problems, the task is to solve for a set of target variables yn associated with a particular time step, usually including the state variables of the dynamic model (xn ). The process is then repeated sequentially to give a time history for yn . Information about yn can come from two sources: evolution of the dynamic model from the previous time step, and comparison between the observations at the current time step (zn ) and the model predictions (h(yn )). Kalman filter. Introduced by Kalman (1960), the Kalman filter is by now a group of algorithms for the sequential combination of dynamic and observational information, using a ‘prediction’ step and an ‘analysis’ step. In the prediction step, the dynamic model is used to calculate prior estimates y _n for the target variables at time step n, from the best (posterior) estimates y ^n1 at the previous step. In the analysis step, posterior estimates y ^n at step n are obtained by ‘improving’ the prior estimates with data. The model state is then ready for evolution to the next (n 1 1) time step. A key point is that the confidence in the current state, embodied in the error covariance for the target variables y, is also evolved with the dynamic model and improved with observations. A schematic diagram of the information flow in the Kalman filter is given in Fig. 1. In the prediction step, the task of evolving y is common across all implementations of the Kalman filter since it involves only a normal forward step of the dynamic model: y _n ¼ u y ^n1 . The prior estimate for the covariance at time step n evolves according to ½Cov y _n ¼ U½Cov y ^n1 UT þ Q; ð7Þ where U 5 @u/@y, the Jacobian matrix of the dynamic model u(y). The first term on the right represents the propagation of the error covariance in the target variables y from one time step to the next, by a linearized version of the model. The second term (Q) is the covariance of the noise term in the dynamic model, Eqn (1), which includes both model imperfections and stochastic variability in forcings and parameters. This term plays a crucial role in the Kalman filter: it quantifies our lack of confidence in the ability of the dynamic model to propagate the model state, and is usually referred to as model error. In most imple￾mentations of the Kalman filter the model error is assumed to be Gaussian with zero mean and no tem￾poral correlation, and thus completely characterized by the covariance matrix Q. In the analysis step, the prior estimates are refined by the inclusion of data. This is done using the prior estimate for the predicted observation vector, z _n ¼ hðy _n Þ, and its covariance ½Cov z _n ¼ H½Cov y _n HT þ R; ð8Þ where H 5 @h/@y is the Jacobian matrix of the observation model h(y), and R is the data covariance matrix [Cov z], indicating lack of confidence in the data and often called the data error. Again it is usually assumed that the data error is Gaussian with zero mean and no temporal correlation, and thus completely characterized by R 5 [Cov z]. The expressions for the final (posterior) estimates for y and its covariance are now exactly as for the nonsequential mode, except that the operation is carried out for one time step only: y ^n ¼y _n þ ½Cov y _n HT½Cov z _n 1 ðzn hðy _n ÞÞ ¼y _n þ Kðzn hðy _n ÞÞ; ð9Þ ½Cov y ^n 1 ¼½Cov y _n 1 þ HTRH or ½Cov y ^n ¼ðI KHÞ½Cov y _n ; ð10Þ K H H Analysis steps zn−1 yn−1 zn zn yn yn Observable State variable State uncertainty Φ Φ, Q H, R difference Prediction step Time step n−1 (posterior) Time step n (prior) Time step n (posterior) Cov yn−1 Cov yn Cov yn Fig. 1 Information flow in the linear Kalman filter. Linear operators are indicated next to arrows. Operations in prediction and analysis steps are shown as dashed and solid lines, respectively. 384 M. R. RAUPACH et al. r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

MODEL-DATASYNTHESISINTERRESTRIALCARBONOBSERVATION385where K = [Covy"jH'[Cov2"j-'is the Kalman gainvariables are effectively the initial state variables forintegration of the model from step n to n+1. Thismatrix. The two equalities in Eqn (10)are equivalent.Time step n is now complete,and we are ready for theapproachfour-dimensionaldataunderpinsnext time step.We note that the ratio of the magnitudesassimilation (4DVAR)methods for assimilating dataof Q and R (model and data error covariances) isinto atmospheric and oceanic circulation models oncritical, since it largely determines how closely theweather and climatetime scales(Chen&Lamb,2000;evolution of y follows that suggested by the dynamicPark&Zupanski,2003).model (Q≤R)or the data.The tuning of Qand R is acrucial part of Kalman filter implementation; seeDiscussion of model-data synthesis methodsGrewal& Andrews (1993)for an excellent extendedDifferences between nonsequential and sequential strategies.discussion.The concepts underlying theKalman filter are nowParameterestimation and dataassimilation problemsimplemented in several different ways (seefor exampletend to beamenableto solution bynonsequential andGrewal&Andrews,1993;Evensen,1994,2003; sequential search strategies, respectively. However, thisKasibhatla et al., 2000; Reichle et al., 2002; Drecourt,is notan absolutecorrespondence:manyproblems can2003),including thefollowing:be solved using either nonsequential or sequential1. The linear Kalman filter (LKF), in which both (y)strategies.The most important advantage for sequentialand h(y)are linear in y,can be shown to be an optimalsolution for appropriate linear problems.methods is the ability of the optimal state to differ2. The extended Kalman filter (EKF) applies forfrom that embodied in the model equations.Thisrequires that the evolving model state x" be includednonlinear p(y)and h(y),bylinearizingthecovariancepropagation part of theanalysis step (Eqn (7),but notamongthetarget variablesyIn principle,y canalsothe prediction step, at each point.This is the algorithminclude x" in nonsequential methods but, since all timesketched above.steps are considered simultaneously,the size of the3.TheensembleKalmanfilter(EnKF)(Evensenproblem is usually intractable.Sequential methods also1994,2003)is appropriatefor high-dimensionalhave the computational advantages that their size doesproblems suchas dataassimilation into atmosphericnot grow with the length of the model integration, andand ocean models, where the error covariance matrixthat they can easily handle incremental extensions tofor y is too large to store, let alone integrate forward.time series observations.The advantages of nonsequential methods come,TheEnKFusesstochasticmethodsbasedonmultiplemodel runs to propagatethe covariancematrix withoutnaturally,fromtheirabilitytotreatalldataat once.Thisstoring it. Also, the EnKF does not explicitly require theisa directadvantagein itself.It is,forexample,difficultJacobian matrices (y) and h(y),which can be difficultfor a sequential method to treat the impact of a datumto derive analytically and expensiveto calculateon a state variable some time in the past, as can occurnumerically. Reichle et al.,(2002) sum marize thewhen, for example, signals are transported through thedifferencesbetweentheEKF and theEnKF.atmosphere so that the model state at sometime is only4.TheKalman smoother assimilatesmultitemporalobserved later. This problem is often handled with theinformation to constrain y"at each time point,byKalmansmoother.running both forward and backward in time (Todling,2000).It produces an estimateoftarget variables attimeModel and data error structures.The noiseterms in thestep n based on the entire record, rather than only thedynamic and observation models, Eqns (1) and (2), canrecord up to timestepn.This gives the Kalmanin principle be quite general in form, including biases,smoother the attributes of a nonsequential method, asdrifts, temporal correlations, extreme outliers and sodata at all times are used together.on. Many extant methods take these noise terms to beGaussian withzeromeanand no temporal correlation,as assumed in Eqns (4)-(10). However, more generalAdjoint methods.These form an additional group oferror.structuresareverycommon,andthemethodsapplicable tosequentialproblems.Thedevelopment of methods for dealing with such errorsprinciple (le Dimet & Talagrand, 1986; Giering 2000)isan activearea of current research.Inthe case ofis to update the target variables (including the modelbiases,drifts and temporal correlations,a promisingstate)by using measurementsat nearby times such asapproach is to introduce extra target variables tothe interval between steps n and n +1, and an estimaterepresent these features of the model or data error.Evensen (2003)showedhowthis approachcanbeusedof thegradientVyJobtained bybackward integration ofan 'adjoint model over that interval. The targetto treat both temporal correlation and bias in themodel2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

where K ¼ ½Cov y _n HT½Cov z _n 1 is the Kalman gain matrix. The two equalities in Eqn (10) are equivalent. Time step n is now complete, and we are ready for the next time step. We note that the ratio of the magnitudes of Q and R (model and data error covariances) is critical, since it largely determines how closely the evolution of y follows that suggested by the dynamic model (Q  R) or the data. The tuning of Q and R is a crucial part of Kalman filter implementation; see Grewal & Andrews (1993) for an excellent extended discussion. The concepts underlying the Kalman filter are now implemented in several different ways (see for example Grewal & Andrews, 1993; Evensen, 1994, 2003; Kasibhatla et al., 2000; Reichle et al., 2002; Dre´court, 2003), including the following: 1. The linear Kalman filter (LKF), in which both u(y) and h(y) are linear in y, can be shown to be an optimal solution for appropriate linear problems. 2. The extended Kalman filter (EKF) applies for nonlinear u(y) and h(y), by linearizing the covariance propagation part of the analysis step (Eqn (7)), but not the prediction step, at each point. This is the algorithm sketched above. 3. The ensemble Kalman filter (EnKF) (Evensen, 1994, 2003) is appropriate for high-dimensional problems such as data assimilation into atmospheric and ocean models, where the error covariance matrix for y is too large to store, let alone integrate forward. The EnKF uses stochastic methods based on multiple model runs to propagate the covariance matrix without storing it. Also, the EnKF does not explicitly require the Jacobian matrices u(y) and h(y), which can be difficult to derive analytically and expensive to calculate numerically. Reichle et al., (2002) sum marize the differences between the EKF and the EnKF. 4. The Kalman smoother assimilates multitemporal information to constrain yn at each time point, by running both forward and backward in time (Todling, 2000). It produces an estimate of target variables at time step n based on the entire record, rather than only the record up to time step n. This gives the Kalman smoother the attributes of a nonsequential method, as data at all times are used together. Adjoint methods. These form an additional group of methods applicable to sequential problems. The principle (le Dimet & Talagrand, 1986; Giering 2000) is to update the target variables (including the model state) by using measurements at nearby times such as the interval between steps n and n 1 1, and an estimate of the gradient HyJ obtained by backward integration of an ‘adjoint model’ over that interval. The target variables are effectively the initial state variables for integration of the model from step n to n 1 1. This approach underpins four-dimensional data assimilation (4DVAR) methods for assimilating data into atmospheric and oceanic circulation models on weather and climate time scales (Chen & Lamb, 2000; Park & Zupanski, 2003). Discussion of model–data synthesis methods Differences between nonsequential and sequential strategies. Parameter estimation and data assimilation problems tend to be amenable to solution by nonsequential and sequential search strategies, respectively. However, this is not an absolute correspondence: many problems can be solved using either nonsequential or sequential strategies. The most important advantage for sequential methods is the ability of the optimal state to differ from that embodied in the model equations. This requires that the evolving model state xn be included among the target variables y. In principle, y can also include xn in nonsequential methods but, since all time steps are considered simultaneously, the size of the problem is usually intractable. Sequential methods also have the computational advantages that their size does not grow with the length of the model integration, and that they can easily handle incremental extensions to time series observations. The advantages of nonsequential methods come, naturally, from their ability to treat all data at once. This is a direct advantage in itself. It is, for example, difficult for a sequential method to treat the impact of a datum on a state variable some time in the past, as can occur when, for example, signals are transported through the atmosphere so that the model state at some time is only observed later. This problem is often handled with the Kalman smoother. Model and data error structures. The noise terms in the dynamic and observation models, Eqns (1) and (2), can in principle be quite general in form, including biases, drifts, temporal correlations, extreme outliers and so on. Many extant methods take these noise terms to be Gaussian with zero mean and no temporal correlation, as assumed in Eqns (4)–(10). However, more general error structures are very common, and the development of methods for dealing with such errors is an active area of current research. In the case of biases, drifts and temporal correlations, a promising approach is to introduce extra target variables to represent these features of the model or data error. Evensen (2003) showed how this approach can be used to treat both temporal correlation and bias in the model MODEL –DATA SYNTHESIS IN TERRESTRIAL CARBON OBSERVATION 385 r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

386M.R.RAUPACH etal.error.Wang &Bentley (2002)introduced a targetet al.,1995;Ciais &Meijer,1998; Enting,1999a,b;Rayner et al., 1999; Rayner, 2001; Schimel et al.,2001;variablerepresentingthetemporally correlated part ofthe data error.Enting,2002;Gurney etal.,2002).In summary,atmo-spheric inversions at global scale provide good con-Nonsequential and sequential parameter estimation.straints on total global sources and sinks,butAlthough parameter estimation is typicallycarried out(presently) with very coarse spatial resolution (con-with nonsequential strategies such as least-squarestinental to hemispheric). In addition to global applica-fitting, there can be advantages in using sequentialtions,atmospheric inversionmethods have beenmethods such as Kalman filteringfor parameterapplied regionally (Gloor et al., 2001),in the atmo-estimation.Theapproach is to treat parameters p asspheric boundary layer (Lloyd et al., 1996), and invegetation canopies (Raupach, 2001).components of the target vectory (in addition tothestate variables x), with p governed by the dynamicA third example, combining the previous two, is theequation dp/dt=0（+noise)(Grewal & Andrewsuse of multiple constraints.This involves model-data1993).Thismeansthattheproblem isalmostalwayssynthesis withthe simultaneous useof multiplekindsnonlinearandmustbesolvedwiththeEKForEnKF.ofobservations(forexample,atmosphericcompositionAnnan & Hargreaves (2004) show how this techniquemeasurements,remotesensing,eddy-covariancefluxes,can be used to estimate parameters in the Lorenzvegetation and soil stores, and hydrological data).Thissystem with chaotic dynamics.Apotential advantageofapproach has two advantages: first, different kinds ofthis approach is that parameters can drift through timeobservation constraindifferentprocesses.Forexampletoward new values,in response to observations.Thisatmospheric composition measurements and eddyoffers a meansformodel-datasynthesis to respond tofluxes directly determine net CO2 exchanges (netecosystem exchange, NEE) at large and small spatialexogenous catastrophic events(such asfire,windthrowor clearing)which suddenly changetheparameters in ascales,respectively,while remote sensing providesterrestrial biosphere model, since exogenous changes inindirect constraints on gross exchanges (gross primaryparameters are the usual way that catastrophic eventsproduction, GPP)through indices such as the normal-ized difference vegetation index (NDVI).Second,areincorporated in the absence of a full dynamicmodelfor theprocessesgoverning the catastrophe.different observationshavedifferentresolutions inspace and time.Through assimilation into a terrestrialbiosphere model, the high space-time resolution ofModel-data synthesis:examplesenvironmental remote sensing can add space-timeThe methods outlined above are being applied intextureto estimates of NEE frommethods suchasseveral fields relevant to terrestrial carbon observation.atmospheric inversions or:eddy-covariance fluxes.Some difficulties must also be noted:for example,The first major exampleis parameter estimation.Mostbiogeochemical models contain parameters(p)deter-handling data sources with quite different spatial andmining photosynthetic capacities,light use efficiencies,temporal scales of measurement (discussed further intemperature and nutrient controls on photosynthesis'Scalemismatchesbetweenmeasurementsandmod-and respiration,pool turnovertimes and soon.It isels'), and also with very different sample numbersalmostalways necessaryto chooseptooptimizethefit(remotely sensed data can swamp in situ data withof the model to test data, usually obtained fromrealistic error specifications, becausethe former hasafactor of 103-10°moredata points).multiple study sites. Techniques for doing this rangefrom simple graphical fits('chi-by-eye')to least-squaresApplications of the multiple-constraint conceptfitting procedures based on Eqn (3) or other costincludethecombineduseofatmosphericCOconcen-functions.trationsand surfacedata at continental scales (Wang&A second example is provided by atmosphericBarrett,2003;Wang &McGregor,2003)andglobalinversion methods for inferring the surface-atmospherescales (Kaminski et al.,2001, 2002); use of geneticfluxes of CO2and othertracegases from atmosphericalgorithmsto constrain terrestrial ecosystem modelscompositionobservations.Thedata comefrom globalof the global carbon cycle with multiple ecologicalflask networks and continuous in situ analysers,upper-data (Barrett, 2002); and discriminating vegetation andair measurements from aircraft and tall towers,andsoil sources and sinks in forest canopies with concen-potentially in the future from remotesensing oftration, isotopic and physiological data (Styles et al.,atmospheric composition.The observation model is a2002).model of global atmospheric transport.The basicA fourth example deserves more space than isapproach has been sketched in Eqns (4) and (5). Thereavailable here: the use of data assimilation in atmo-is now a significant literature on this technique (Entingspheric and ocean circulation models.This is now2005Blackwell Publishing Ltd, Global Change Biology,11,378-397

error. Wang & Bentley (2002) introduced a target variable representing the temporally correlated part of the data error. Nonsequential and sequential parameter estimation. Although parameter estimation is typically carried out with nonsequential strategies such as least-squares fitting, there can be advantages in using sequential methods such as Kalman filtering for parameter estimation. The approach is to treat parameters p as components of the target vector y (in addition to the state variables x), with p governed by the dynamic equation dp/dt 5 0 ( 1 noise) (Grewal & Andrews 1993). This means that the problem is almost always nonlinear and must be solved with the EKF or EnKF. Annan & Hargreaves (2004) show how this technique can be used to estimate parameters in the Lorenz system with chaotic dynamics. A potential advantage of this approach is that parameters can drift through time toward new values, in response to observations. This offers a means for model–data synthesis to respond to exogenous catastrophic events (such as fire, windthrow or clearing) which suddenly change the parameters in a terrestrial biosphere model, since exogenous changes in parameters are the usual way that catastrophic events are incorporated in the absence of a full dynamic model for the processes governing the catastrophe. Model–data synthesis: examples The methods outlined above are being applied in several fields relevant to terrestrial carbon observation. The first major example is parameter estimation. Most biogeochemical models contain parameters (p) deter￾mining photosynthetic capacities, light use efficiencies, temperature and nutrient controls on photosynthesis and respiration, pool turnover times and so on. It is almost always necessary to choose p to optimize the fit of the model to test data, usually obtained from multiple study sites. Techniques for doing this range from simple graphical fits (‘chi-by-eye’) to least-squares fitting procedures based on Eqn (3) or other cost functions. A second example is provided by atmospheric inversion methods for inferring the surface-atmosphere fluxes of CO2 and other trace gases from atmospheric composition observations. The data come from global flask networks and continuous in situ analysers, upper￾air measurements from aircraft and tall towers, and potentially in the future from remote sensing of atmospheric composition. The observation model is a model of global atmospheric transport. The basic approach has been sketched in Eqns (4) and (5). There is now a significant literature on this technique (Enting et al., 1995; Ciais & Meijer, 1998; Enting, 1999a,b; Rayner et al., 1999; Rayner, 2001; Schimel et al., 2001; Enting, 2002; Gurney et al., 2002). In summary, atmo￾spheric inversions at global scale provide good con￾straints on total global sources and sinks, but (presently) with very coarse spatial resolution (con￾tinental to hemispheric). In addition to global applica￾tions, atmospheric inversion methods have been applied regionally (Gloor et al., 2001), in the atmo￾spheric boundary layer (Lloyd et al., 1996), and in vegetation canopies (Raupach, 2001). A third example, combining the previous two, is the use of multiple constraints. This involves model–data synthesis with the simultaneous use of multiple kinds of observations (for example, atmospheric composition measurements, remote sensing, eddy-covariance fluxes, vegetation and soil stores, and hydrological data). This approach has two advantages: first, different kinds of observation constrain different processes. For example, atmospheric composition measurements and eddy fluxes directly determine net CO2 exchanges (net ecosystem exchange, NEE) at large and small spatial scales, respectively, while remote sensing provides indirect constraints on gross exchanges (gross primary production, GPP) through indices such as the normal￾ized difference vegetation index (NDVI). Second, different observations have different resolutions in space and time. Through assimilation into a terrestrial biosphere model, the high space-time resolution of environmental remote sensing can add space-time texture to estimates of NEE from methods such as atmospheric inversions or eddy-covariance fluxes. Some difficulties must also be noted: for example, handling data sources with quite different spatial and temporal scales of measurement (discussed further in ‘Scale mismatches between measurements and mod￾els’), and also with very different sample numbers (remotely sensed data can swamp in situ data with realistic error specifications, because the former has a factor of 103 –106 more data points). Applications of the multiple-constraint concept include the combined use of atmospheric CO2 concen￾trations and surface data at continental scales (Wang & Barrett, 2003; Wang & McGregor, 2003) and global scales (Kaminski et al., 2001, 2002); use of genetic algorithms to constrain terrestrial ecosystem models of the global carbon cycle with multiple ecological data (Barrett, 2002); and discriminating vegetation and soil sources and sinks in forest canopies with concen￾tration, isotopic and physiological data (Styles et al., 2002). A fourth example deserves more space than is available here: the use of data assimilation in atmo￾spheric and ocean circulation models. This is now 386 M. R. RAUPACH et al. r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397

MODEL-DATASYNTHESISINTERRESTRIALCARBONOBSERVATION387well-developed and applied routinely in weatherselection of observations from four categories of data:forecasting.A variety of techniques are employed,remotesensingofland surfaceproperties,atmosphericincluding'nudging, three-dimensional and four-composition measurements, direct flux measurements,dimensional variational dataassimilation (3DVARand direct measurements of carbon stores.Theaim is toand 4DVAR)based onadjointmethods,and useofmake estimates of error properties for these categoriesthe ensemble Kalman filter.Recent reviews arepro-ofmeasurement.Thediscussiondoesnotaddressallofvided by Chen &Lamb (2000)and Park &Zupanskithe above issues, largely omitting questions of spatial(2003).and temporal error structure. We present tablesindicating ranges for the diagonal elements [Covz]mm =Cm2 of the error covariance matrix for measurementData characteristics: uncertainty in measurementerror,andthequalitativebehaviourof the correlationsand representationwhich determine the off-diagonal elements.The entriesin thesetables aremostly'expert judgements'by theWe have emphasized that data uncertainties affect notauthors and their colleagues, backed up by quantitativeonly the predicted uncertainty of the eventual result ofevidencewherepossible.Thereis,of course,noclaima model-data synthesis process, but also the predictedthat our estimates are definitive; the intention is ratherbest estimate.This realization raises the challenge ofto indicate the kinds of uncertainty informationevaluatingthe uncertaintyproperties of the mainkindsrequired of observations for model-data synthesisof observationrelevanttoa TCOS,informsdirectlypurposes.The tables characterize measurement errorsusable for model-data synthesis.This is a very largeonly;representation errors, which often exceed mea-goal,which embraces all categories of observationsurementerrors,arediscussedseparatelyinqualitativeidentified at the beginning of Model-data synthesis:terms onlyThe issue of scale mismatches betweenMethods',and also a range of issues:measurements and models, which arises in all cases asThe error magnitudeom for an observation Zma significant contribution to representation error,isinclusive of all error sources (in other words, thetreated generically in a fifth subsection Scale mis-diagonal elements [Covz]mm=Cm? of the covariancematches between measurements and models'.matrix);:The correlations [Covzlmn/(omon)among errors inRemote sensing of land surface propertiesdifferent observations, quantified by the off-diagThe main satellite-borne remotely sensed data on landonal elements ofthecovariancematrix;surface properties comefromtwokinds ofsensor,both. The temporal structure of the errors: whether theypolar-orbiting to provide frequent global coverage:are random in time or temporally correlated,andmoderate-spatial-resolution (~250-1000m)and high-the possiblepresence of unknown long-term driftstemporal-resolution (~1dayrepeatinterval)sensorsorbiases;such as AVHRR and MODIS; and high-spatial resolu-:The spatial structure of errors (random, slowlytion (~10-30 m)and moderate-temporal-resolutionvarying or bias as for temporal structure);(~16dayrepeat interval)sensorssuchasSPOTand.The error distribution: normal (Gaussian), log-LANDSAT. All these sensors provide multi-year re-normal,skewedorthesumofmultipleerrorsourcescords.Onemajor application(amongmany)forthewith different distributions, such as a small Gaus-AVHRR-MODIS family is assessment of vegetationsian noise together with occasional large outliersdynamics with indices such as NDVI (defined asbecause of measurement corruption events;(NIR-Red)/(NIR+Red), where NIR and Red areradiances in the near-infrared and visible red spectral.Possiblemismatches between the spatial and tem-bands) and measures such as surface temperature.poral averaging implicit in the model and theApplications for the SPOT-LANDSAT family includemeasurements(the‘scalingproblem');detection of land cover change and vegetation clearing.The separate contributions to all the above errorand regrowth. In all cases, the measurements are at-properties of measurement error (the distribution ofsensor reflected radiances from the earth in severalthe measurements z around their true values) andspectralbands(5forAVHRR,37forMODIS).representation error (the distribution of theerror inIn usingthese formsof remotesensingdataforthe model representation of the measurement,model-data synthesis applications,three kinds of errorz = h(y).need to be considered:(1) errors associated with theThis challenge is too large to meet fully here. To make ameasurement and spatial attribution of radiances at thestart, we consider (in the next four subsections)asensor; (2) errors in relating radiances at sensor to2005 Blackwell Publishing Ltd, Global Change Biology,11, 378-397

well-developed and applied routinely in weather forecasting. A variety of techniques are employed, including ‘nudging’, three-dimensional and four￾dimensional variational data assimilation (3DVAR and 4DVAR) based on adjoint methods, and use of the ensemble Kalman filter. Recent reviews are pro￾vided by Chen & Lamb (2000) and Park & Zupanski (2003). Data characteristics: uncertainty in measurement and representation We have emphasized that data uncertainties affect not only the predicted uncertainty of the eventual result of a model–data synthesis process, but also the predicted best estimate. This realization raises the challenge of evaluating the uncertainty properties of the main kinds of observation relevant to a TCOS, in forms directly usable for model–data synthesis. This is a very large goal, which embraces all categories of observation identified at the beginning of ‘Model–data synthesis: Methods’, and also a range of issues:
The error magnitude sm for an observation zm, inclusive of all error sources (in other words, the diagonal elements ½Cov z mm ¼ sm2 of the covariance matrix);
The correlations [Cov z]mn/(smsn) among errors in different observations, quantified by the off-diag￾onal elements of the covariance matrix;
The temporal structure of the errors: whether they are random in time or temporally correlated, and the possible presence of unknown long-term drifts or biases;
The spatial structure of errors (random, slowly varying or bias as for temporal structure);
The error distribution: normal (Gaussian), log￾normal, skewed or the sum of multiple error sources with different distributions, such as a small Gaus￾sian noise together with occasional large outliers because of measurement corruption events;
Possible mismatches between the spatial and tem￾poral averaging implicit in the model and the measurements (the ‘scaling problem’);
The separate contributions to all the above error properties of measurement error (the distribution of the measurements z around their true values) and representation error (the distribution of the error in the model representation of the measurement, z 5 h(y)). This challenge is too large to meet fully here. To make a start, we consider (in the next four subsections) a selection of observations from four categories of data: remote sensing of land surface properties, atmospheric composition measurements, direct flux measurements, and direct measurements of carbon stores. The aim is to make estimates of error properties for these categories of measurement. The discussion does not address all of the above issues, largely omitting questions of spatial and temporal error structure. We present tables indicating ranges for the diagonal elements ½Cov z mm ¼ sm2 of the error covariance matrix for measurement error, and the qualitative behaviour of the correlations which determine the off-diagonal elements. The entries in these tables are mostly ‘expert judgements’ by the authors and their colleagues, backed up by quantitative evidence where possible. There is, of course, no claim that our estimates are definitive; the intention is rather to indicate the kinds of uncertainty information required of observations for model–data synthesis purposes. The tables characterize measurement errors only; representation errors, which often exceed mea￾surement errors, are discussed separately in qualitative terms only. The issue of scale mismatches between measurements and models, which arises in all cases as a significant contribution to representation error, is treated generically in a fifth subsection ‘Scale mis￾matches between measurements and models’. Remote sensing of land surface properties The main satellite-borne remotely sensed data on land surface properties come from two kinds of sensor, both polar-orbiting to provide frequent global coverage: moderate-spatial-resolution (  250–1000 m) and high￾temporal-resolution (  1 day repeat interval) sensors such as AVHRR and MODIS; and high-spatial resolu￾tion (  10–30 m) and moderate-temporal-resolution (  16 day repeat interval) sensors such as SPOT and LANDSAT. All these sensors provide multi-year re￾cords. One major application (among many) for the AVHRR-MODIS family is assessment of vegetation dynamics with indices such as NDVI (defined as (NIRRed)/(NIR 1 Red), where NIR and Red are radiances in the near-infrared and visible red spectral bands) and measures such as surface temperature. Applications for the SPOT-LANDSAT family include detection of land cover change and vegetation clearing and regrowth. In all cases, the measurements are at￾sensor reflected radiances from the earth in several spectral bands (5 for AVHRR, 37 for MODIS). In using these forms of remote sensing data for model–data synthesis applications, three kinds of error need to be considered: (1) errors associated with the measurement and spatial attribution of radiances at the sensor; (2) errors in relating radiances at sensor to MODEL –DATA SYNTHESIS IN TERRESTRIAL CARBON OBSERVATION 387 r 2005 Blackwell Publishing Ltd, Global Change Biology, 11, 378–397