Transcript Document

Carbon and water flux measurements and
applications to regional biogeochemistry Eddy covariance revisited
David Hollinger
USDA Forest Service
Northern Research Station
Durham, NH 03824 USA
NCAR ASP 2007 Colloquium
Regional Biogeochemistry: Needs and Methodologies
June 4-15, 2007
Boulder, CO
National Center for Atmospheric Research
Advanced Study Program
A biased history of Eddy Covariance –
Approx number of research
groups using EC
~500
~100
~20
~5
Limited by
technique
development
- not clear
how to make
it work
1970s
Limited by
instrument
access –
The basic
approach
works, but
instruments
(fast IRGAs,
sonics) not
widely
available
1980s
Technique
explodes –
Companies
offer very
Limited by
easy,
ease of use turnkey
–Instruments systems
become
available,
but the
technique is
still
cumbersome
and requires
a lot of effort
to learn
1990s
Saturation –
The
challenge
becomes
finding
important,
unanswered
science
questions
Why ?
2000s
1. There’s always a need for top quality data !
2. To use the data effectively, must understand the
assumptions
Organization:
• Turbulence and Eddies
• Eddy Covariance Methodology
– History (Baldocchi 2003, GCB 9:479)
– Theory
– Practice
– Using other people’s data – how do you
know the data are any good?
• Uncertainty
• Crucial assumptions
• Methods
• Key problems
Philosophy:
• All models are wrong
– Useful simplifications
• All data are corrupt
– Noise (random component)
– Bias (fixed component, compared to an
independent measure)
– Variable bias
• But, we need to move forward…
Advantages of eddy flux data for evaluating
and parameterizing models
• Same or finer timestep as models
• Simultaneous measurements of CO2, water, &
energy fluxes
• Environmental data also recorded
– PFD, temperature, soil moisture, etc.
• Often abundant site data
– LAI, leaf %N, C pool data, etc.
• Area averaging, whole ecosystem net flux
Eddy flux data useful for the study of
biogeochemical processes are collected
worldwide and readily available
• 294 active FLUXNET sites worldwide
• 57% with woody vegetation
• Available data on-line:
– AmeriFlux
380 site-years
– EUROFLUX
51 site-yearsa
– Fluxnet-Canada, Asiaflux, Ozflux, etc.
• In a few years we’ve gone from a data poor to a
data rich environment. How do we get to a
knowledge rich environment?
a. Additional data available to network participants or by special arrangement
Turbulence and Eddies
“…an irregular condition of flow in which the
various quantities show a random variation
with time and space coordinates, so that
statistically distinct average values can be
discerned” (Hinze, 1975)
Turbulence can be visualized as eddies (local
swirling motions). Over forests, eddies may
form as coherent “transverse rollers”
(Raupach et al. 1996)
Raupach, M. et al. 1996. Boundary-Layer Met. 78:351-82.
Turbulence and Eddies
.
Finnigan, J. 2000. Turbulence in plant canopies. Annual Review of Fluid Mechanics 32: 519-571
• Eddies are carried by the mean flow (turbulence is not local)
• Eddies overlap in space, with large eddies carrying smaller
eddies.
• Turbulent kinetic energy transfers from large eddies to
smaller eddies in a cascading process and is dissipated as
heat energy once molecular viscosity becomes important.
“Big whorls have little whorls that feed on
their velocity,
And little whorls have lesser whorls, and so
on to viscosity.”
-- Lewis F. Richardson
The supply of energy from and to Atmospheric Eddies' (1920)
Proceeding of the Royal Society A 97 354-73.
Great fleas have little fleas upon their backs to bite 'em, And little fleas have lesser fleas, and so ad infinitum. And
the great fleas themselves, in turn, have greater fleas to go on; While these again have greater still, and greater
still, and so on. -- Augustus De Morgan
Turbulence and Eddies
.
Finnigan, J. 2000. Turbulence in plant canopies. Annual Review of Fluid Mechanics 32: 519-571
The size of an eddy () is its “turbulent length scale”,
above a forest canopy this may be 10’s-1000’s of meters
• For Monin-Obukhov (surface layer) scaling,   z/u
z is measurement height, u is mean windspeed
(Bussinger et al. 1971, J. Atmos. Sci 28:181; Wyngaard et al. 1971, J.
Atmos. Sci 28:1171)
• For Raupach-Finnigan (mixing layer) scaling,   h/u
h is vegetation height
Turbulence integrates sources/sinks and
connects them to the measurement system:
Population
>>106
10’s – 100’s of meters
Transport
uncertainty
+ Measurement
uncertainty
Eddy Method: Simplification of
the Continuity Equation
• Fluxes consist of mean and fluctuating
components (Reynold’s decomposition)
• The mean components (advection) = 0
• The only relevant fluctuating component
is with the vertical wind
Baldocchi et al. (1988) Ecology 69:1331-40.
Aubinet et al. (2000) Adv. Ecol. Res., 30:113-175
Massman and Lee (2002) Ag and For Met 113:121-44
c uc
c  wc c
F u

w


x
x
z
z
t
z
Overbar signifies time average,
primes deviations from mean
x
u is horizontal wind,
w is vertical wind
Crucial assumption: F = eddy flux + Dstorage
c uc
c  wc c
F u

w


x
x
z
z
t
z
Overbar signifies time average,
primes deviations from mean
x
u is horizontal wind,
w is vertical wind
The covariance
• F = w’c’
• Cov(w,c) = E{[ w - E(w) ][ c - E(c) ]}
– E() expected value (mean)
– Insensitive to uncorrelated noise
Cospectra show the distribution of the flux in the
frequency domain:
Co w't' (fx)/w't'
1
0.2
0.1
0.02
0.01
0.001
(z-d) = 3.8 m (Schidler)
(z-d) = 27 (Morgan Monroe)
0.01
0.1
1
10
fx
Tall tower, low wind
Short tower, high wind
Methods
(Massman & Lee 2002, Ag. For. Met. 113:121.
Papale et al. 2006, Biogeosciences 3:571)
•
•
•
•
Measure vertical wind & scalar (lag scalar)
No net w; rotate coordinates as necessary
Correct for spectral deficiencies in measurement
Correct for density fluctuations (WPL, selfheating)
• Reject bad data (variances too high/low, bad
calibrations,u* threshold, etc.)
• Fill in gaps for annual sums
• Software:
– MacMillen (ATDD), EddySol, EddyRe, TK2, etc.
ATI & Gill
Crucial Assumption – Site:
No horizontal or vertical advection
• Homogeneous surface
roughness
exchange characteristics
• Flat
• Extensive
old rule: 100h
• No local sources of CO2
avoid buildings and roads
Meteorological conditions
affect turbulence:
• a u* threshold ensures
there will be sufficient
mixing
T
Ri  g
z

z
T u
z

Stable - bad

2
T
Unstable - good
Neutral - ok
z
z
T
T
Crucial Assumptions Meteorology
• Turbulence is present
 no measurement
without transport
- u* (mixing) threshold (Goulden, 1996)
• Conditions are “stationary” (statistical
properties not changing or “not changing too
quickly”)
– Stationarity Tests (Foken & Wichura 1996, Ag For Met
78:83)
Measurement Diagnostics
• Does the energy balance?
Rn = H + LE + G + DS
• Is there flux divergence?
Are fluxes the same at different heights?
• Do the power spectra and cospectra look
o.k.? (raw data)
• Does the flux plateau with u*?
• Is the intercept of the light response similar to
the mean nocturnal respiration.
Nocturnal respiration (g m-2 y-1 )
Relative sample size
and r2
600
1996
1997
1998
1999
500
400
(a)
300
0
0.1
0.2
0.3
0.4
u* threshold
1
n (1996)
r2 (1996)
0.8
0.6
0.4
(b)
0.2
0
0.1
0.2
u* threshold
0.3
0.4
Compare nocturnal data with parameters
estimated from daytime data:
NEE = PmaxI / (Km+I) + Rd
Another solution:
Are the data any good?
• Use data from long established sites.
• See what sites keep appearing in
synthesis papers.
• Avoid sites with “mountain”, “ridge”, etc.
in the name.
What is uncertainty and why do we care?
• All measurements are
corrupted by random
error:
We want the flux, F, but we
measure F + d + e
where d is the random error
and e is a bias
• The random error is
characterized by its PDF;
the first moment (s)
referred
to as
“uncertainty”
-10
-5
0
5
10
•
Uncertainty information is
required to properly fit
models to data; e.g. for
estimating the
parameters of canopy
processes models.
 Maximum Likelihood
 Data assimilation
(Kalman filter)
 Bayesian methods
•
Uncertainty useful for
assessing model fits,
annual uncertainties, risk
analysis, etc.
Maximum likelihood – “given the data, what are the
most likely model coefficients (parameters)?”
• Determined by minimizing the difference between
data and model via a likelihood (cost) function:
• The likelihood function requires uncertainty
information (s) from the data
N
 
2
i
N
 
2
i
( yi  y( xi ))2
s
2
For Gaussian data
i
yi  y( xi )
2i
For double exponential data
Where do we get uncertainty info?
What is uncertainty info?
"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."
– Raupach et al. (2005). Model data synthesis in
terrestrial carbon observation: methods, data
requirements and data uncertainty specifications.
Global Change Biology 11:378-97.
Flux Uncertainties are nonGaussian with nonconstant variance
•
•
•
1.
2.
3.
Simultaneous measurements
at 2 towers (Howland)1
Single tower next day
comparisons (Howland,
Harvard, Duke, Lethbridge,
WLEF, Nebraska)2
Data-model residuals3
=28.7
=31.6
Hollinger & Richardson (2005) Tree Physiology 25: 873-885.
Richardson et al. (2006) Ag. Forest Met. 136: 1-18.
Hagen et al. (2006), Journal of Geophysical Research 111,
D08S03.
• a double exponential PDF
better represents the random
error distribution of eddy fluxes
2   s , where  
yˆ  yi
n
=2.7
Flux uncertainty
increases linearly
with flux magnitude
~proportional to flux
(for least squares uncertainty
is constant)
Forests
FC > 0
2 = 0.62 + 0.63*FC
FC < 0
2 = 1.42 - 0.19*FC
Richardson et al. (2006) Ag. Forest Met. 136: 1-18.
Maximum likelihood – “given the data, what are the
most likely model coefficients (parameters)?”
• Probability
• Determined by minimizing the difference between
data and model via a likelihood function:
• The likelihood function depends on the uncertainty
characteristics (PDF) of the data
N
 
2
i
N
 
2
i
( yi  y( xi ))2
s
For Gaussian data
2
i
yi  y( xi )
2i
For double exponential data
Choices affect the outcome:
Weighting of data model deviations
Rd = Rref e(q/(T-T0))
Richardson & Hollinger 2005
Ag. For Met 131:191
Inverse modeling - choosing the
“best” respiration model
(Richardson et al. 2006)
• Linear
q1  q2T
• Q10
q1  q 2
(T Tref )/10

• temperature varying Q10

• time varying Q10
• Arrhenius
N
2  
i
• Lloyd and Taylor
s2
i
• Fourier
• Neural Networks
(T Tref ) /10
q1  q2  q3T 
(T Tref ) /10
q1q2  q3 sin(JD )  q4 cos(JD )
 
q 2  283.15 
q1  exp
1 

283
.
15

T






( yi  y( xi ))2


q2
q1  exp

T  273.15 q3 
• Lloyd and Taylor (restr)
  308.56 

 T  46.02 
q1  exp
Ranking respiration models based on
eddy flux or automated chamber data
Results: Q10 and restricted Lloyd & Taylor
performed poorly
Harvard
Harvard
0.23
0.21
0.20
0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.18
0.18
0.18
0.18
0.18
0.17
NA
NA
NN-SAJ
NN-SA
NN-A
Q10-vTime
P olynomial
Q10-A
Lloyd&Taylor
NN-S
Linear
Logistic
Q10-vTemp
Arrhenius
Q10-S
Fourier-4
L&T-Rest
Fourier-2
Fourier-1
NN-SAJW
Q10-Gresp
0.19 NN-SA
0.19 NN-A
Howland
HowlandHowl
below
canopy
Howl an d-Mai n
an d-S u bcan opy
0.38
0.38
0.36
0.35
0.34
0.34
0.34
0.34
0.32
0.32
0.32
0.30
0.30
0.30
0.30
0.30
0.29
0.26
0.25
NN-SAJW
NN-SAJ
NN-SA
Fourier-4
Fourier-2
Q10-vTime
P olynomial
NN-S
Lloyd&Taylor
Q10-vTemp
Logistic
NN-A
Linear
Q10-Gresp
Arrhenius
Q10-S
Q10-A
Fourier-1
L&T-Rest
0.34 NN-SA
0.33 NN-SAJ
0.21
0.21
0.20
0.20
0.19
0.19
0.19
0.19
0.19
0.19
0.19
0.18
0.18
0.18
0.18
0.17
0.17
0.16
0.15
NN-SAJW
NN-SAJ
NN-SA
Fourier-4
P olynomial
Q10-vTime
NN-S
Fourier-2
Q10-vTemp
Logistic
Lloyd&Taylor
Q10-Gresp
Arrhenius
Linear
Q10-S
L&T-Rest
Fourier-1
NN-A
Q10-A
0.18 Q10-vTemp
0.18 Lloyd&Taylor
Howland
soil
Howl an d-Au toch am be
0.54
0.53
0.50
0.47
0.46
0.46
0.46
0.46
0.46
0.46
0.46
0.46
0.43
0.42
0.41
0.38
0.36
0.26
0.23
NN-SAJW
NN-SAJ
NN-SA
Q10-vTime
Q10-Gresp
P olynomial
NN-S
Q10-vTemp
Logistic
Arrhenius
Q10-S
Lloyd&Taylor
Linear
Fourier-4
Fourier-2
L&T-Rest
Fourier-1
NN-A
Q10-A
0.46 NN-SA
0.44 Q10-S
The better the model fit, the higher the estimated annual
respiration – bias from using different models?
Cost functions?
• There can be many different cost functions,
each yielding different parameters and
integral estimates
• Be humble
Words to remember with Eddy Covariance?
• Be aware of the assumptions
– Turbulence present (stationary)
– No advection
• Reject bad data - Use good data
•Large corrections are
necessary due to
density effects of
sensible & latent heat
fluxesa
•FCO2 is thus affected
by measurement error
and uncertainty in H
and LE
• Corrections for self
heating also necessary
FCO2 Correction ( mol m-2 s-1)
A Problem with
openpath sensors:
20
H
LE
15
10
5
0
-5
-100
0
100
200
300
400
Energy flux (W m-2 )
a. Webb, E.K., G.I. Pearman, and R. Leuning. 1980. Correction of flux measurements for density effects due to heat and water
vapor transfer. Quart. J. Met. Soc. 106: 85-100
500
Flux Correction ( mol m-2 s-1)
WPL (density)
Corrections:
30
Coniferous Forest (2001)
20
10
0
-10
-20
~ 0.025 (0.325 LE+1.684 H)
~
0
100
200
300
JD
FCO2 ( mol m-2 s-1)
20
10
0
-10
-20
-30
0
100
200
JD
300
Self heating
Corrections:
(Burba et al. AMS 2006)
0
0
0.2
0.4
0.6
0.8
1
1.2
u*
FCO2 ( mol m-2 s-1)
(c)
10
slightly stable (4)
stable (5)
slightly unstable (3)
5
very stable (6)
unstable (2)
0
very unstable (1)
0
0.2
0.4
0.6
u*
0.8
1
1.2
Inverse modeling: A simple canopy model example
driven by PFD, (Tair +Tsoil)/2, & DW
PFD
NEE  Amax
e
( K m  PFD)
NEE =
Amax –maximum PS
Km – half saturation
Topt – temperature optimum
– temp shape factor
– DW shape factor
Ra – base respiration
E0 – activation energy
X
+
 T Topt 


  
4
e
 w w 
 i a 
  
2
 Ra  e
  E0 


 T  30. 
X
All 2004 flux
data from
Howland &
Bartlett
PFD
NEE  Amax
e
( K m  PFD)
 T Topt 


  
4
e
 w w 
 i a 
  
2
 Ra  e
  E0 


 T  30. 
data
N
 
2
i
( yi  y( xi ))
2
s2
i
uncertainty
By guided trial and error, find the parameters that minimize the
sum of the weighted deviations between the model and
measurements (2
The time to do this is linear in N and factorial in parameter
number
Example results: The 95% confidence intervals of the
parameter sets for Howland and Bartlett are far apart
PFD
NEE  Amax
e
( K m  PFD)
 T Topt 


  
4
e
 w w 
 i a 
  
2
 Ra  e
  E0 


 T  30. 
Functional responses for Howland and Bartlett based on
most likely parameter estimates are distinct
Determining parameter
covariances help to
understand model
behavior:
PFD
NEE  Amax
e
( K m  PFD)
 T Topt 


  
4
e
 w w 
 i a 
  
2
 Ra  e
Bartlett
Howland
  E0 


 T  30. 
Determining measurement uncertainty by
multiple measurements of the same thing
•
Compare independent measurements of the same thing (x
and y) - The surface must be homogenous
•
Calculate mean and standard deviation of the difference
between the towers: qi = xi - yi (~ 0) sq  0
y1
If x and y are independent
and have the same random
measurement uncertainty:
y2
x1
x2
q  x y
s q  s x2  s y2
assume s x  s y
Howland
sx 
sq
2
then
Openpath gas analyzers have superior
high frequency response and thus
smaller corrections:
CO2 flux correction
1.5
1.4
1.3
z = 10 m
1.2
z = 20 m
1.1
z = 40 m
1
openpath
0
2
4
6
8
u (m s-1 )
on
1.5
1.4
z = 20 m
10
Corrections have Implications for energy
balance closure & u* corrections!
How do we minimize the likelihood function?
• Traditional non-linear approaches (LM)
– Available in SAS and other statistical packages
– Often fail if parameters poorly constrained or model
contains local minima
• If functions have local minima, we need a global
method
– Available in MATLAB or program form
– Simulated annealing (SIMANN, Goffe et al. 1994 )
– Genetic algorithms (SRES, Runarrson & Yao 2000)
(See Press et al. “Numerical recipes”, chapters 10 & 15)
Fitted Pmax ( mol CO2 m-2 s-1)
Independent data sets from separate towers
yielded similar parameter variation:
0
Main tower
West tower
-10
-20
-30
r = 0.97, P<0.001
1999
2000
Season
2001
Respiration flux ( mol m-2 s-1)
Lloyd & Taylor Respiration
8
Measured respiration
2nd-order Fourier
Lloyd & Taylor
6
Simple Respiration models:
4
2
0
J FMAM J J A SOND
8
a4 cos(2d )  a5 sin(2d )
8 7
9
6
4
10
6
5
11
2
4
12
231
0
0
2
4
6
• Fourier series
Unbiased, good fit
No need for other data
Not mechanistic
R  a1  a2 cos(d )  a3 sin(d ) 
8
2nd-order Fourier Respiration
Lloyd and Taylor (1994) Functional Ecology 8:315-23
• Lloyd and Taylor
Possibly biased, good fit
Mechanistic
Need complete temperature data
  Eo


R  Ae

(T To ) 
History of Model Inversions with flux data
• First to use eddy flux data to determine parameters
of canopy models
– Wang et al. (2001)
• 10 parameters (Leuning & Wang 1998), 3 weeks data (pasture) X
6 sites
• 3-4 parameters constrained
– Schulz et al. (2001)
• 21 parameters, 2 weeks flux data (wheat), 3 parameters
• Bayesian inversions – prior distribution specified
– Braswell et al. (2005)
• SIPnET model, 23 parameters, 10 years day/night data Harvard
• Constrained 12 parameters, 7 edge hitting
– Knorr and Kattge (2005)
• Bethy model, 24 parameters, 7 days of forest flux data
• Constrained 5 parameters, evaluated parameter covariances