Research Vignette: The TransCom3 Time-Dependent Global CO 2 Flux Inversion … and More

1. Research Vignette: The TransCom3 Time-DependentGlobal CO2 Flux Inversion … and More David F. Baker NCAR 12 July 2007

2. Outline • The TransCom3 CO2 flux inversion inter-comparison project • The fully time-dependent T3 flux inversion • Method (“batch least squares”) • Results • Methods for bigger problems: • Kalman filter (traditional, full rank) • Ensemble filters • Variational data assimilation

3. TransCom3 CO2 Flux Inversions CO2 fluxes for 22 regions, data from 78 sites • Annual-mean inversion, 1992-1996 • Fixed seasonal cycle, no IAV • 22 annual mean fluxes solved for • Gurney, et al Nature, 2002 & GBC, 2003 • Seasonal inversion, 1992-1996 • Seasonal cycle solved for, no IAV • 22*12 monthly fluxes solved for • Gurney, et al, GBC, 2004 • Inter-annual inversion, 1988-2003 • Both seasonal cycle and IAV solved for • 22*12*16 monthly fluxes solved for • Baker, et al, GBC, 2006

4. TransCom3 interannual inversion setup: Solve for monthly CO2 fluxes from 22 regions, 1988-2003, using monthly concentrations from 78 sites

6. Fossil fuel input Atmospheric storage 0.45*FF Land + ocean uptake FF = Atmos + Ocean + Land Bio

8. NH Winter NH Summer CO2 Uptake & Release

9. Form of the Batch Measurement Equations 0 H fluxes Transport basis functions concentrations

11. Batch least-squares or “Bayesian synthesis” inversion ^ Optimal fluxes, x, found by minimizing: where giving

12. Computation of the interannual variability (IAV): Europe Monthly flux [PgC/yr] Deseasonalized flux IAV Estimation error Between-model error

13. Total flux IAV (land+ocean) [PgC/yr]

14. Ocean Flux IAV [PgC/yr] Land flux IAV

15. [PgC/yr]

16. Flux IAV [PgC/yr], 11 land regions

17. Flux IAV [PgC/yr], 11 ocean regions

18. Computational considerations • Transport model runs to generate H: • 22 regions x 16 years x 12 months x 36 months = 152 K tracer-months (if using real winds) • 22 x 12 x 36 = 9.5 K tracer-months (using climatological winds) • Matrix inversion computations: O (N3) • N = 22 regions x 16 years x 12 months = 4.4 K • Matrix storage: O (N*M) --- 66 MB • M = 78 sites x 16 years x 12 months = 15 K

19. Kalman Filter/Smoother 0 H

20. Kalman Filter/Smoother 0 H

21. 1 2 1 1 1 1 1 2 2 2 2 2 Error Time Kalman Filter Equations Measurement update step at time k: Dynamic propagation step from times k to k+1: Put multiple months of flux in state vector xk, method becomes effectively a fixed-lag Kalman smoother • =tangent linear model

22. Inversion methods for the data-rich, fine-scale problem • Kalman filter: some benefit, but long lifetimes for CO2 limit savings • Ensemble KF: full covariance matrix replaced by an approximation derived from an ensemble • Variational data assimilation (4-D Var): an iterative solution replaces the direct matrix inversion; the adjoint model computes gradients efficiently

23. Ensemble Kalman filter • Replace xk, Pk from the full KF with an ensemble of xk,i, i=1,…,Nens • Add dynamic noise consistent with Qk to xk,i when propagating; add measurement noise consistent with Rk to measurements when updating, initial ensemble has a spread consistent with P0 • When needed in KF equations, Pk replaced with • Replace matrix multiplications with sums of dot products • Good for non-Gaussian distributions

24. Kalman filter vs. Kalman smoother For retrospective analyses, a 2-sided smoother gives more accurate estimates than a 1-sided filter. (Gelb, 1974)

25. Estimation as minimization • Solve for x with an approximate, iterative method rather than an exact matrix inversion • Start with guess x0, compute gradient efficiently with an adjoint model, search for minimum along -, compute new  and repeat • Good for non-linear problems; use conjugate gradient or BFGS approaches • Low-rank covariance matrix built up as iterations progress • As with Kalman filter, transport errors can be handled as dynamic noise

26. 4-D Var Iterative Optimization Procedure 0 x0 ° Forward Transport Forward Transport Estimated Fluxes “True” Fluxes 1 2 Modeled Concentrations “True” Concentrations ° x1 ° Measurement Sampling Measurement Sampling x2 Modeled Measurements “True” Measurements Assumed Measurement Errors D/(Error)2 3 x3 ° Weighted Measurement Residuals Flux Update Adjoint Transport Adjoint Fluxes =  Minimum of cost function J

27. 4-D Var Data Assimilation Method • Find optimal fluxes u and initial CO2 field xo to minimize • subject to the dynamical constraint • where • x are state variables (CO2 concentrations), • h is a vector of measurement operators • z are the observations, • R is the covariance matrix for z, • uo is an a priori estimate of the fluxes, • Puo is the covariance matrix for uo, • xo is an a priori estimate of the initial concentrations, • Pxo is the covariance matrix for xo, • is the transition matrix (representing the transport model), and G distributes the current fluxes into the bottom layer of the model

28. 4-D Var Data Assimilation Method Adjoin the dynamical constraints to J using Lagrange multipliers  SettingF/xi = 0 gives an equation for i, the adjoint of xi: The adjoints to the control variables are given byF/uiandF/xooas: The gradient vector is given by F/uiandF/xoo, so the optimal u and xo may then be found using one’s favorite descent method. I use the BFGS method in order to obtain an estimate of the leading terms in the covariance matrix.

29. · 10-8 [ kg CO2 m-2 s-1 ] OSSE Results For Five CO2 Measurement Networks

30. Pros and cons, 4DVar vs. ensemble Kalman filter (EnKF) • 4DVar requires an adjoint model to back-propagate information -- this can be a royal pain to develop! • The EnKF can get around needing an adjoint by using a filter-lag rather than a fixed-interval Kalman smoother. However, the need to propagate multiple time steps in the state makes it less efficient than the 4DVar method • Both give a low-rank estimate of the a posteriori covariance matrix • Both can account for dynamic errors • Both calculate time-evolving correlations between the state and the measurements

31. * * * Adjoint transport model • If number of flux regions > number of measurement sites, then instead of running transport model forward in time forced by fluxes to fill H, run adjoint model backwards in time from measurement sites • What is an adjoint model? • If every step in the model can be represented as a matrix multiplication (= ‘tangent linear model’), then the adjoint model is created by multiplying the transpose of the matrices together in reverse order * measurement sites * FWD ADJ * flux grid