Guidelines for Design and Diagnostics of CO 2 Inversions

1. Guidelines for Design and Diagnostics of CO2 Inversions Anna M. Michalak Department of Civil and Environmental Engineering, and Department of Atmospheric, Oceanic and Space Sciences, The University of Michigan

2. Ill-Conditioned Nature of Inverse Problems

3. Environmental Contamination with Unknown Sources Source: http://www.marshfieldclinic.org/nfmc/lab/research_projects.stm#

4. Example – Available Measurements

5. Ownership of Site

6. Source Release Scenarios

7. Blue’s Source Release Scenario

8. Green’s Source Release Scenario

9. Red’s Source Release Scenario

10. Possible Scenarios

11. Propagation of Uncertainty

12. Bayesian Inference Applied to Inverse Modeling for Inferring Historical Forcing Likelihood of forcing given available measurements Posterior probability of historical forcing Prior information about forcing p(y) probabilityofmeasurements y : available observations (n×1) s: discretized historical forcing (m×1)

13. Bayesian Formalism • Use data, y, prior flux estimates, sp, and model (with Green’s function matrix H) to estimate fluxes, s • Estimate obtained by minimizing: • Solution is • Estimates, ŝ have covariance • Residuals:

14. Impact of Aggregation and Independence Assumptions

15. Modeling Tools for NA Carbon Cycle Actual flux history Available data

16. Modeling Tools for NA Carbon Cycle Geostatistical Bayesian 31 data 201 fluxes 31 data 101 fluxes 31 data 41 fluxes 31 data 11 fluxes 31 data 21 fluxes

17. Statistical Diagnostics of Inversions(Ongoing work with Ian Enting)

18. Need for Diagnostics • Wide use of inversion studies • Large set of possible results due to differences in: • transport models • inversion methods (e.g. Bayesian, geostatistical, mass balance) • data choices • meteorological fields • covariance and other parameter choices • TransCom experiments aim to assess variability and derive (relative) consensus • Moving toward operational inversions • Need objective method to evaluate inversions to determine (at a minimum) which inversions are self-consistent

19. Approaches to Inversion Validation • Cumulative plots of residuals (Enting et al., 1995) • Reporting of 2 statistics of residuals from priors and/or observations (e.g. Rayner et al., 1999; Gurney et al., 2002; Peylin et al., 2002; Rödenbeck et al., 2003) • Variance of observation residuals calculated using conditional realizations of a posteriori fluxes (Michalak et al., 2004) • Maximum likelihood approach leading to r2 = 1 (Michalak et al., 2005; Hirsch et al., 2006) • Statistical diagnostics project proposed during TransCom-Tsukuba

20. Bayesian Formalism • Use data, y, prior flux estimates, sp, and model (with Green’s function matrix H) to estimate fluxes, s • Estimate obtained by minimizing: • Solution is • Estimates, ŝ have covariance • Residuals:

21. Testing Residuals • Want to test assumption of zero-mean, multivariate normal with covariance Q and R for and • Unknown strue in r0,sp and r0,y • Sum of squares of normalized residuals from fit, ŝ: should be distributed as 2 with n degrees of freedom, i.e. normalized r1,sp, r1,y are not n + m independent N(0,1) quantities

22. Testing Residuals Residuals from fit, r1,sp and r1,y are correlated because they represent departures from the ‘model’ line.

23. Testing Residuals • Conditional realizations can be generated:where and uk is the k-th realization of vector of N(0,1) values. • Normalized residuals from sc,k : should have covariances R and Q

24. Testing Residuals Residuals from fit, r1,spand r1,y are correlated because they represent departures from the ‘model’ line. Adding perturbations uleads to residuals r2,spand r2,y that are independent

25. Testing Residuals Residuals from fit, r1,sp and r1,y are correlated because they represent departures from the ‘model’ line. Adding perturbations u leads to residuals r2,sp and r2,y that are independent

26. Simple Tests • Do residuals have the specified covariance structure? • Are the residuals unbiased? • Are the residuals normally distributed? • Perform normality tests such as Kolmogorov-Smirnov goodness-of-fit hypothesis test, Lilliefors hypothesis test of composite normality, Filliben normality test, etc.

27. Test Setup • Initial scoping study with 1995 CSIRO setup: • GISS model (8o x 10o) • Cyclo-stationary inversion: constant + fixed seasonal cycle • 12 ocean regions, 12 land regions with (separate) mean plus pre-specified seasonality, 4 regions of (season-dependent) deforestation, fossil source plus explicit CO oxidation • Observations expressed as mean + Fourier components of seasonal cycle • No use of O2 or 13C data • Examined cases: • Reference case, Biased priors, Biased non-fossil priors, Loose data

28. Results Fails KS normality test Fails 2 test Fails unbiasedness test Fails 2 and unbiasedness tests Fails normality and unbiasedness tests Passes all tests

29. Results – Reference Case

30. Results – Biased Priors

31. Results – Biased Non-fossil Priors

32. Results – Loose Data

33. What if Residuals are Correlated? • Q and R are no longer diagonal matrices, but we can use geostatistical methods to calculate equivalent statistics • Do residuals have the specified covariance structure? • Are the residuals unbiased?

34. Conclusions • Residual analysis should be a standard step in validating inversion results • Conditional realizations allow for simple residual tests and tests on subsets of residuals • Diagnostics will not detect errors due to mis-interpretation (CO2 flux ≠ carbon flux ≠ carbon storage rate) • Geostatistics provides a set of tools for dealing with spatially and/or temporally correlated errors and parameters • Many cases suggest that previous studies have used cautious assignments of uncertainty, motivated by risk of unknown correlated errors.

35. Next Steps • Develop additional tests • Analyze residuals from individuals stations / regions • Investigate use of loose priors for “reluctant” Bayesians • Analysis of large ensembles of conditional realizations • Application to existing TransCom inversions

36. Maximum Likelihood Approach for Covariance Parameter Estimation

37. Maximum Likelihood Estimation of Covariance Parameters • Covariance parameters determine: • Relative weight assigned to data and prior information • Posterior covariance / uncertainty estimate • Appropriate estimates of covariance parameters is essential to flux estimation • Lack of objective methods for estimating these parameters: • Described as “greatest single weakness” in some studies (Rayner et al., 1999) • Maximum likelihood approach provides estimates based on available data

38. Model-data Mismatch v. Prior Error Michalak et al. (JGR 2005)

39. Procedure • Define the marginal pdf of covariance parameters (w/out solving inverse problem). Its negative logarithm is: • Minimum of objective function determines best estimate of covariance parameters • Inverse of Hessian of objective function estimates uncertainty of covariance parameters • Currently applied to diagonal matrices, but can directly be used with matrices with off-diagonal terms (i.e. correlated residuals) • Need to define (numerically or analytically) derivative w.r.t. covariance parameters • Demonstrated in the geostatistical framework (Michalak et al., JGR 2004) Michalak et al. (JGR 2005)

40. Stations Used in ML Study Michalak et al. (JGR 2005)

41. Constant Variances R = 1.63 ppm Q = 2.17 GtC/yr Michalak et al. (JGR 2005)

42. Variances Based on Physical Attributes R, MBL = 0.71 ppm R, HI/DES = 1.49 ppm R, CONT = 3.16 ppm Q, LAND = 2.02 GtC/yr Q, OCEAN = 1.07 GtC/yr Michalak et al. (JGR 2005)

43. Variances Based on Inversion Behavior R,1 = 0.58 ppm R,2 = 1.04 ppm R,3 = 1.19 ppm R,4 = 4.36 ppm R,5 = 7.64 ppm Q, LAND = 1.76 GtC/yr Q, OCEAN = 1.21 GtC/yr Michalak et al. (JGR 2005)

44. Variance Based on Auxiliary Information R = CS = 0.11 – 5.96 ppm Q, LAND = 2.08 GtC/yr Q, OCEAN = 0.89 GtC/yr Michalak et al. (JGR 2005)

45. Land Fluxes Michalak et al. (JGR 2005)

46. Ocean Fluxes Michalak et al. (JGR 2005)

47. Land Uncertainty Michalak et al. (JGR 2005)

48. Ocean Uncertainty Michalak et al. (JGR 2005)

49. Ocean Uncertainty Michalak et al. (JGR 2005)

50. Conclusion from ML Method • Data themselves can provide information about model-data mismatch and prior error covariances (in the absence of external information regarding residual covariances) • Covariances R and Q reflect different patterns of residuals • Maximum likelihood approach produces covariance estimates consistent with physical understanding • ML can be applied to more complex covariance structures