Transcript Document

Lagrangian tracer transport models:
Measuring concentrations, assessing emissions
Formulation
Adjoint
Mass Conservation and other complications
S. C. Wofsy, Harvard University, 25 April 2009
(1) Particle dispersion models;
(2) Adjoint implementation at Harvard:
Stochastic Time-Inverted Lagrangian
Transport (STILT) Model;
(3) Application to assessment of sources of
greenhouse gases to the atmosphere.
GV aircraft approaches the
terminator over the Alaska
range in January, 2009, in
HIPPO_1 experiment to define
transport rates and global
sources of greenhouse gases
Global and regional studies of CO2 , other greenhouse
gases, ozone destroying species, etc, use observations of
concentrations at the ground, and from aircraft and
satellites, to infer emission rates from surface sources.
An inverse problem.
Why do we want to know?
• Emissions of these gases drive changes in radiative forcing, ozone
loss, atmospheric chemistry
• Processes regulating emissions, and overall budgets are much less
certain than needed to predict future atmospheric composition
• Natural vs. human—caused balances is debated
• Future monetization of CO2 emissions requires verification of
source reductions
NOAA Greenhouse Gas records
55% of annual
emissions
Basic Lagrangian Equation: Trajectories
The advection of a particle or puff is computed from the
average of the three-dimensional velocity vectors for the
initial-position P(t) and the first-guess position P'(t+∆t). The
velocity vectors are typically linearly interpolated in both
space and time. The first guess is
P'(t+∆t) = P(t) + V(P,t) ∆t
and the final position is
P(t+∆t) = P(t) + 0.5 [ V(P,t) + V(P',t+∆t) ] ∆t,
If V is the mean wind vector from an assimilated meteorological
field, the equation provides a trajectory model. It is fully reversible
in time, because particles (“air parcels”) do not spread. Trajectory
models cannot compute concentrations because the particles have
zero measure.
Puff and Dispersion Models
Puff model: the source releases pollutant puffs at regular
intervals. Each puff is advected according to the trajectory of
its center of mass, as it expands (both horizontally and
vertically) to simulate dispersion in a turbulent atmosphere.
Useful for short term simulations.
Lagrangian particle dispersion model: the source is
simulated by releasing many particles at the same time. A
random component to the motion is added at each step
according to atmospheric turbulence ( unresolved motions).
Useful for longer term simulations.
Both require that stability and mixing coefficients
be computed from meteorological data.
Lagrangian approach to atmospheric chemistry
models : Particle Dispersion Modeling
•Air parcels are represented by particles advected with the
mean wind plus stochastic velocities that simulate turbulent
transport, convection, and other unresolved/"non-conservative"
motions of the atmosphere.
•Particle positions are resolved to arbitrary accuracy, reducing
grid-averaging errors such as representation error (for
observations) and source spreading (for emissions).
•Particles may be followed backwards in time, providing a direct
adjoint model.
Lagrangian particle dispersion models, and the analogous puff
dispersion models, look superficially like trajectory models (same
basic equation), but there are fundamental differences. Reversibility
follows a statistical criterion.
Time reversal and the advectiondiffusion equation for mass transport
n c/t = (Unc) + F
Time reversal velocity reversal: t´  -t ; U´  - U ;
F = <U'c'>, does not change sign !
Diffusion does not reverse—puffs or ensembles of
particles expand, forward or backward in time!
The direction of the velocity in a random walk varies
randomly, thus u'  -u changes nothing.
Dispersion  Random walk
X2/t  n2 steps to move nDL
Receptor-oriented model: timeinverted (adjoint) model
Particles tell where air is arriving
from and link observations with
upstream sources/sinks
Receptor-oriented model: timeinverted (adjoint) model
Backward Lagrangian transport
(x, y, z)
receptor
(xi’,yi’,zi’)
source
Forward = Backward?
(x, y, z)
receptor
(xi’,yi’,zi’)
source
Backward
forward
STILT (Stochastic Time Inverted Lagrangian Transport Model)
u'
stochastic
(turbulent)
advection
u
Receptor
location
mean
advection
u’ depends on:
1) TL (decorrelation timescale)
2) sw (fluctuations in turbulent velocity)
 Based on HYSPLIT (Hybrid Single Particle Lagrangian Integrated Trajectory) model code
[Draxler and Hess, 1998]
 Driven by ETA, AVN (forecasts) or EDAS, GDAS, NGM (assimilations)
Improved turbulence parameterization
(“off-line” model)
•TLw (vertical) and sw after Hanna [1982]
•reflection/transmission scheme at interfaces between high and low turbulence after
Thomson [1997]
Framework: Lagrangian Particle Dispersion
Model Framework: A particle ensemble is released at the receptor
and transported backward in time, tracking air parcels that arrive (in
the forward-time sense) at the receptor at a given time.
The density of STILT particles is used to calculate the influence
function I(xr,tr|x,t) and the footprint f(xr, tr| x,t). I(xr,tr|x,t) and f(xr, tr|
x,t) link concentration measurements at the receptor, C(xr, tr), to the
sum of all upstream contributions:
C(xr, tr) =
tr
)   dt  d 3 x I (x r , tr | x, t ) S (x, t )   d 3 x I ( x r , tr | x, t0 )C ( x, t0 )
t0
V
V
Units: C is in ppm; I is vol-1 (density); S is ppm s-1;
Influence Function, I
I(xr,tr|x,t) is represented by the density r(xr,tr|x,t) of particles at (x,t)
which were transported backward in time from (xr,tr), normalized by
the total particle number Ntot:
I (xr , tr | x, t ) 
r (xr , tr | x, t )
Ntot
1

Ntot
Ntot
 ( x (t )  x)
p
p=1
The fields of I(xr,tr|x,t) and S(x, t), continuous in space and time, are in
practice resolved only at finite resolution with a discrete volume (Dx, Dy,
Dz) and finite time interval (t). Integrating over the finite volume and
time elements to get the time/volume integrated I:
tm t

tm
y j Dy
xi Dx
dt

xi
dx

yj
zk Dz
dy

zk
1
dz I (xr , tr | x, t ) 
Ntot
1

Ntot
tm t

y j Dy
xi Dx
dt
tm

dx
xi
Ntot
 Dt
p 1
p , m ,i , j , k

yj
zk Dz
dy

zk
dz r (xr , tr | x, t )
Source-receptor matrix elements
in I (tr, xr, yr, zr | ts, xs, ys, zs )
(finite small volume, time interval…)
The time- and volume-integrated influence function is simply
quantified by tallying Dtp,m,i,j,k, the total amount of time each
y j Dy p spends
y j D
y
The
Dx particle
zk Dz in a volume element i,j,k over
tm t timestep
xi Dx m.
zk Dz
1 sources at finite temporal
source-receptor
matrix
elements
that
link
dx
dy
dz I (xr , tr | x, t ) 
dt
dx
dy
dz
and spatial resolutions directly to receptor
N concentrations is:

yj

zk
tot
1
I (m´, i ´, j ´, k ´ | m, i, j, k ) 
Ntot

tm


xi
yj
p , m ,i , j , k
/(DV Dt)
Ntot
 Dt
p 1
Model should run equally well forwards or backwards…
now we should check it!

zk
Representing Surface Sources; mass
conservation (Part 1)
 F ( x, y, t )  mair for z  h

S (x, t)   h  r ( x, y, t )
0
for z  h

F is surface flux (e.g. mmole m-2 s-1), h a near-surface mixing depth (m),
r the mean density of dry air in h (kg m-3), mair the molecular mass of dry
_air (kg mole-1), giving S in ppm s-1.
Particles are effectively tagged with a mole fraction of pollutant, and the
mole fraction is conserved as the particle is transported to different
pressure levels. Since S is defined consistently in the interior or at the
boundary, and each receptor is at a single pressure, particles are
counted with the correct weight (mass is conserved) regardless of where
they are emitted or observed.
Footprint function, f: surface influence
S(x,t) is integrated over discrete volume and time elements, defining the
footprint function, f :
y j Dy
tm t
xi Dx
h


mair
DCm,i , j (x r , tr )  
  dt  dx  dy  dz I (x r , tr | x, t )   F ( xi , y j , tm )
 h  r ( xi , y j , tm ) tm

xi
yj
0
 f (x r , tr | xi , y j , tm )  F ( xi , y j , tm )
mair
1
f (xr , tr | xi , y j , tm ) 

h  r ( xi , y j , tm ) Ntot
Ntot
 Dt
p 1
p ,i , j ,k
In our application, f(xr, tr| xi, yj, tm) is in units of [ppm/(mmole/m2/s)], i.e.
given a unit surface flux of 1 mmole/m2/s at (xi, yj, tm) persisting over a
time interval t, f(xr, tr| xi, yj, tm) is the concentration change DC in ppm at
the receptor.
Detailed Linkage between Observation and
Upstream Sources/Sinks (“Footprint”) from STILT
Footprint of Vertical Profile over WLEF from LT1600 Particle Release on June 8th, 1999
-12 hrs
-6 hrs
-24 hrs
log10(footprint)
-8
-6
-4
-2
log10[ppmv/(mmole/m2/s)]
Physical Requirements for Realistic Simulation of
Source-Receptor Relationships
Simulations of transport and upstream influence by particle models
must satisfy the following criteria:
1) maintenance of an initially well-mixed state
2) simulation of the close interaction between wind shear and
vertical turbulence
3) high temporal resolution to resolve the decay in the
autocorrelation of u'
4) consistent representation of particles as air parcels of equal mass
in both the mean and turbulent transport components of the
model.
Turbulent transport: 1st order Markov process
Hanna [1979] showed from Eulerian and Lagrangian
measurements that a Markov assumption is reasonable, i.e. that the
particle velocity vector u can be decomposed into a mean
component and a turbulent component u', with the turbulent
component following the relation:
u' (t+Dt)=R(Dt)u'(t)+u'' (t)
u'' is a random vector, R an autocorrelation coefficient: R(Dt) =
exp(–Dt/TLi), where TLi is the Lagrangian time scale (decorrelation
time) (i=u: horizontal; i=w: vertical). The random velocity u'' is:
u'' = l[1-R2(Dt)]1/2
with l drawn from a Gaussian distribution, and standard deviation
si, is derived from the meteorological field (e.g. from TKE, or
mixed-layer height+surface fluxes of heat, momentum).
Beware of Looking Under the Hood:
Issues Addressed in Time-Reversibility
Comparisons
•Wind coming out of ground
•Failure to incorporate vertical density gradients
•Outliers from random Gaussian generator
•Operator splitting
•Particles trapped in low-turbulence regions
In this house we obey the
Laws of Thermodynamics!
-- Homer Simpson
Problem: Particles become trapped in lowturbulence regions!
(spurious entropy loss)
altitude
Particles accumulate here.
sw
more turbulent
This problem is closely
analogous to the
problem of detailed
balance in the kinetic
theory of gases.
Satisfying the Well-Mixed Criterion
by Treating Turbulence
Profiles as Separate Layers
z
altitude
sw
sw
s w ( zi  )
s w ( zi  ) r ( zi  )
where wt  wi
, =
s w ( zi  )
s w ( zi  ) r ( zi  )
Problem: Particles
become trapped in
regions of spurious
mass loss!
Actual simulation of
particle densities for
a set of winds from
EDAS-80.
Violates 2nd law of
Thermodynamics…
Backward- vs Forward-time
Comparison
NGM Mass-violating Winds
2.0
NGM Mass-violating Winds
for one starting time
1.5
orthogonal
regression
1.0
1st order mass
correction
0.5
1:1 Line
0.0
Particle Density in Source Box [#/km3]
(Backward)
slope=1.852 R2=0.805
0.0
0.2
0.4
0.6
0.8
Particle Density in Receptor Box [#/km3]
Day=15
(Forward)
source
receptor
source
creation of
mass
receptor
source
destruction receptor
of mass
Implications of Mass Violation
•Affects lots of atmospheric models, especially ‘off-line’ models
•Careful attention needed in using wind fields from other models!
•Forward-time not necessarily better than backward-time results
•Resulting time-irreversibility may be an issue for ‘adjoint’ models
as well
Backward- vs Forward-time Comparison
Mass-Conserving
Synthetic Wind FieldWinds
0.4
Mass-conserving, Simplified Winds
0.3
orthogonal
regression
0.2
1:1 Line
0.1
slope=0.98±0.03
R2=0.97
0.0
Particle Density in Source Box [#/km3]
(Backward)
slope=0.979 R2=0.969
0.0
0.1
0.2
0.3
Particle Density in Receptor Box
(Forward)
[#/km3]
0.4
Applications of STILT
1. CH4, N2O, CO at ground stations, from aircraft.
2. CO2 sources and sinks globally
Quantifying Anthropogenic Combustion Contributions
to C Budget
log influence
ppm/(mmol-2s-1)
CO ppb
250
model
obs
200
150
7
8
9
10
Date (August 1999)
log influence
ppm/(mmol-2s-1)
Example: CO at
Harvard Forest, 1999.
Courtesy C. Gerbig,
J. Lin
STILT shows that: Accurate measurements from a
single tall tower "cover" a remarkably large area.
log10
Influence: ppm per unit flux (mmole m-2s-1) at WLEF 396 m
Scott M. Miller (2008) STILT
Methane and Nitrous Oxide
in North America: Using an LPDM to Constrain
Emissions
Eric Kort
[email protected]
Non-CO2 Workshop
October 23, 2008
Case Study- COBRA-NA 2003
•
~300 flasks measured @ NOAA/Boulder, UND Citation II, 23 May to 28 June 2003
1600 1700 1800 1900 2000
CH4
0
50
100 150 200 250 300
Flask number
June 16, 2003
36.72N,96.94W 609 m AGL
WRF
LPDM Model: STILT
Emissions: EDGAR-2000
Met fields: WRF (AER, 35
km, LPDM outputs, Grell-2)
Footprint * Prior
Emission Field
(EDGAR)
*
Result =
Enhancement of
gas at measurement
point due to source
• Methane
Prior Emissions Fields
– AnthropogenicEDGAR32FT2000
– Biogenic- Jed Kaplan
wetland inventory
• Nitrous Oxide
FIRES
– AnthropogenicEDGAR32FT2000
– Anthropogenic & BiogenicGEIA
Measurements- Footprint
Results- Methane
Slope: 0.924 ± 0.13
Scaling Factor: 1.08 ± 0.15
Note: Prior Emissions Field
EDGAR32FT 2000 & JK wetland
1860
1840
1820
1800
1780
Measurement Error:
Atmospheric Variability in a
correlation length
(170 km) derived from
CO measurement &
CO:CH4 ratio:
2s = 23 ppbv
PDF for model CH4 concentrations
for 100 particles from one receptor:
Mean: 1790 ppb; 2s = 38 ppbv
1760
(model, each particle)
CH4 (ppb)
STILT LPDM provides directly applicable
error estimates for Bayesian inverse.
WRF/STILT/EDGAR
model vs data, with gray
and green.
Errors used in fitting are +
38 ppbv for the model,
and + 23 ppbv for the
measurements
Slope = 0.9
sslope = 0.1
EDGAR—2000
confirmed ±10% for
CH4 !
This result pertains to
urban-industrial
sources, which
dominate the flight
region
N2O Observed vs Model (EDGAR—2000 )
Model STILT/ N2O (ppbv)
COBRA-2003
Observedthan
N2O (ppbv)
US sources of N2O are ~2.5x higher
EDGAR Kort et al., 2008
Carbon Monoxide
CO,  3!
EPA Inventory  3!
International dateline (-180 Longitude)
Summary—Lagrangian particle
dispersion models for source assessment
• Lagrangian particle dispersion models provide
powerful tools to utilize atmospheric observations to
study sources and sinks of atmospheric greenhouse
gases.
• Their unique capability to run forward and
backward in a consistent numerical framework
permits checks on model transport not available any
other way.
• Considerable effort is needed to ensure proper
formulation of this type of model. Rigorously massconserving wind fields are required.
• Applications to aircraft measurements over North
America reveal important aspects of greenhouse gas
emissions on continental scale.
R component exercise: make a movie; visualization of a random walk
#generate a random walk movie; N particles in a random walk along the x axis
N=500 ; M = 1000 # 100 particles, up to 1000 time steps
X=rep(0,N) #N particles start at location 0
x=matrix( sample(c(-1,1),N*M,replace=T) , ncol=N) # matrix w/ each row random steps of +/-1
X1=rbind( X, apply(x, 2, cumsum)) ; dimnames(X1)= list(as.character(0:M),as.character(1:N))
# is this a Markov process ? what about staying in place?
#matplot(t(X1)) ##very simple animation
X11(); plot(0:1000,apply(X1,1,var)) ### the variance grows linearly with the number of steps (time)
X11(); plot(0:1000,apply(X1,1,function(x){sqrt(var(x))}))
X11(); plot(sqrt(0:1000),apply(X1,1,function(x){sqrt(var(x))}))
#
#make the movie
for(i in 0:M){
s.n=(paste("000",i,sep="")); ser=substring(s.n,nchar(s.n)-3,nchar(s.n))
if(i%%100==0){yll=range(hist(X1[i+1,],breaks= seq(
- 3.5*sqrt(M),3.5*sqrt(M),length.out=120))$counts)}
png(file=paste("random_walk_movie/Fig.random_walk_",ser,".png",sep="") )
hist(X1[i+1,],breaks= seq(-3.5*sqrt(M),3.5*sqrt(M),length.out=120),xlim=c(-3.5,3.5)*sqrt(M),ylim=yll)
dev.off()
}
30
1000
25
apply(X1, 1, function(x) {
sqrt(var(x))
10
15})
20
800
600
5
400
0
0
200
2
variance
(d
1, var)
apply(X1,
N)
s vs time
(#steps)
0
0
200
400
600
800
1000
200
400
600
0:1000
10
5
0
Frequency
15
Histogram
of X1[i + 1, ]
0:1000
distance
(#steps)
-100
-50
0
distance
50
100
X1[i + 1, ]
distance
800
1000