Geostatistics Principles of spatial analysis Anna M Michalak

Geostatistics: Principles of spatial analysis Anna M. Michalak Department of Civil and Environmental Engineering Department of Atmospheric, Oceanic and Space Sciences The University of Michigan

Key Points n If the parameter(s) that you are modeling exhibits spatial (and/or temporal) autocorrelation, this feature must be taken into account to avoid biased solutions n Spatial (and/or temporal) autocorrelation can be used as a source of information in helping to constrain parameter distributions n The field of geostatistics provides a framework for addressing the above two issues A. M. Michalak (amichala@umich. edu)

Outline n n n n Motivation for geostatistical tools What is geostatistics? Traditional applications Application to OCO sampling design Introduction to inverse modeling Application to groundwater contamination Application to CO 2 flux estimation A. M. Michalak (amichala@umich. edu)

What is Geostatistics? n A short answer: q An interpolation and extrapolation toolkit n A more sophisticated answer: q All of the above for modeling spatial relationship of available data and building from such a model (e. g. kriging, stochastic simulation, …) n Formal definition q Analysis and prediction of spatial or temporal phenomena (e. g. pollutant concentrations, soil porosities, elevations, etc. ) A. M. Michalak (amichala@umich. edu)

Spatial Correlation n Measurements in close proximity to each other generally exhibit less variability than measurements taken farther apart. n Assuming independence, spatially-correlated data may lead to: 1. Biased estimates of model parameters 2. Biased statistical testing of model parameters n Spatial correlation can be accounted for by using geostatistical techniques A. M. Michalak (amichala@umich. edu)

Parameter Bias Example map of an alpine basin Q: What is the mean snow depth in the watershed? snow depth measurements kriging estimate of mean snow depth (assumes spatial correlation) mean of snow depth measurements (assumes spatial independence) A. M. Michalak (amichala@umich. edu)

Example cont… H 0 is TRUE 5% H 0 rejected A. M. Michalak (amichala@umich. edu) 5% H 0 Rejected! H 0 H 5% 0 Not Rejected

Variogram z(x)Model = m(x) n Used to describe spatial correlation 1 2 3 4 A. M. Michalak (amichala@umich. edu) + e(x)

Geostatistics in Practice n Main uses: q Data integration q Numerical models for prediction q Numerical assessment (model) of uncertainty A. M. Michalak (amichala@umich. edu)

Caveats DOESN’T Provide practical solution to real problems Fully automate estimation process Honor data Replace good or additional data Expand from data Create data Integrate data Provide causal / physical relationships Save time & effort Geostatistics is a set of decision-making tools A. M. Michalak (amichala@umich. edu)

Steps in Geostatistical Study n Exploratory Data Analysis (EDA) q Data cleaning q Consistency of data q Identification of populations n Spatial Continuity Analysis q Experimental q Analysis, interpretation q Quantitative n Estimation q Uncertainty assessment q Account for spatial correlation q Integrate hard and soft information n Simulation q Alternative images of the field q Reproduce field heterogeneity q Honor all available information A. M. Michalak (amichala@umich. edu)

Go to Matlab… A. M. Michalak (amichala@umich. edu)

OCO Satellite q Planned launch in September 2008 q Will provide global columnintegrated CO 2 measurements q 1 ppm measurement accuracy at a 1000 km scale. A. M. Michalak (amichala@umich. edu)

OCO Measurements q 1 ppm measurement accuracy at a 1000 km scale. q Processing all spectral radiances to XCO 2 is computationally prohibitive. q Limit Sampling to optimal locations A. M. Michalak (amichala@umich. edu)

OCO Subsampling Strategy n Objective: q Determine optimal sampling locations as a function of time and space that allow for the interpolation of XCO 2 at unsampled locations with estimation error within a set threshold n Recent work: q Define modeled XCO 2 spatial variability using CASA-MATCH data (Olsen and Randerson 2004) subsampled at 1 pm local time q Preliminary approach for identifying optimal sampling locations A. M. Michalak (amichala@umich. edu)

Sample Modeled XCO 2 Data April July August October A. M. Michalak (amichala@umich. edu)

Optimal Sampling Locations q Optimal sampling locations = potential sampling locations that will achieve a set estimation error threshold at unsampled locations q Estimation error = estimation standard deviation at unsampled locations q Geostatistical interpolation tools: n Use spatial correlation as a basis of estimation n Provide best linear unbiased estimates n Quantify associated estimation error A. M. Michalak (amichala@umich. edu)

Spatial correlation (Variogram model) h 6 h 2 h 1 h 4 h 3 Semivariance, γ(h) h 5 2 1 4 3 6 5 Separation Distance, h A. M. Michalak (amichala@umich. edu)

Global Spatial Variability ½ variance Correlation Length A. M. Michalak (amichala@umich. edu)

Global Spatial Variability A. M. Michalak (amichala@umich. edu)

Local Variability (2000 km radius) 5. 5 degrees 2000 km A. M. Michalak (amichala@umich. edu)

XCO 2 Variance and Correlation Length April Correlation length (km) A. M. Michalak (amichala@umich. edu) Variance (ppm 2)

Distance to Achieve 1 ppm Uncertainty (h 0) h 0 =? Vmax=1 ppm A. M. Michalak (amichala@umich. edu) n h 0 = max distance from the interpolation point to sample for 1 ppm error n h 0 depends on spatial variability near interpolation point n Interpolation at each grid point on a 5. 5 o by 5. 5 o global grid

Maximum Sampling Interval h 0 - April Maximum sampling interval (km) A. M. Michalak (amichala@umich. edu)

Regular Grid Sampling Uncertainty July A. M. Michalak (amichala@umich. edu) April

Optimal Sampling Locations and Associated Uncertainties July A. M. Michalak (amichala@umich. edu) April

Sampling Constraints A. M. Michalak (amichala@umich. edu) n n n n Aerosols Clouds Satellite track Maximum (sub)sampling rate Albedo Measurement error Temporal aggregation Others?

Conclusions from OCO Study n n XCO 2 exhibits strong spatial correlation XCO 2 covariance structure is variable in space and time Uniform sampling will not achieve uniform/acceptable interpolation uncertainty Geostatistical tools can be used to incorporate the variability in the XCO 2 covariance structure into a subsampling protocol A. M. Michalak (amichala@umich. edu)

Inverse Modeling A. M. Michalak (amichala@umich. edu)

Inverse models n Geostatistical inverse modeling objective function: H = transport information s = unknown fluxes y = CO 2 measurements R = model-data mismatch covariance Q = spatial/temporal covariance of flux deviations from trend X and = model of the trend Deterministic component A. M. Michalak (amichala@umich. edu) Stochastic component

Bayesian Inference Applied to Inverse Modeling for Inferring Historical Forcing Posterior probability of historical forcing Likelihood of forcing given available measurements y : available observations (n× 1) s : discretized historical forcing (m× 1) A. M. Michalak (amichala@umich. edu) Prior information about forcing p(y) probability of measurements

Dover Air Force Base Case Study n Dover Air Force Base located in Delaware, U. S. A. n Unconfined aquifer underlain by two-layer aquitard n Aquitard cores used to infer PCE and TCE contamination history in aquifer n Solute transport controlled by diffusive process: A. M. Michalak (amichala@umich. edu)

TCE at Location PPC 11 Measured TCE concentration as a function of depth A. M. Michalak (amichala@umich. edu) Time variation of boundary condition

TCE at Location PPC 13 Measured TCE concentration as a function of depth A. M. Michalak (amichala@umich. edu) Time variation of boundary condition

Sources of Atmospheric CO 2 Information North American Carbon Program A. M. Michalak (amichala@umich. edu)

What Surface Fluxes to Atmospheric Samples See? Longitude Source: Arlyn Andrews, NOAA-GMD A. M. Michalak (amichala@umich. edu) Height Above Ground Level (km) Latitude 24 June 2000: Particle Trajectories -24 hours -48 hours -72 hours -96 hours -120 hours Longitude

Large Regions Inversion Trans. Com 3 Sites & Basis Regions Trans. Com, Gurney et al. (2003) A. M. Michalak (amichala@umich. edu)

Study Goals 1. Estimate carbon fluxes at fine spatial resolution (3. 75 o x 5. 0 o) 2. Avoid use of prior flux estimates 3. Incorporate and quantify effect of available auxiliary data Questions: n What will be the effect on estimated fluxes and their uncertainties? n Is there sufficient information in the atmospheric measurements to “see” the relationship between auxiliary data and fluxes? A. M. Michalak (amichala@umich. edu)

Auxiliary Data and Carbon Flux Processes: Other: Spatial trends Anthropogenic Flux: Fossil fuel combustion (sine latitude, absolute value latitude) (GDP density, population) Environmental parameters: (precipitation, %land use, Palmer drought index) Oceanic Flux: Gas transfer Terrestrial Flux: Photosynthesis (sea surface (FPAR, LAI, NDVI) temperature, air temperature) Respiration (temperature) Image Source: NCAR A. M. Michalak (amichala@umich. edu)

Sample Auxiliary Data A. M. Michalak (amichala@umich. edu)

Global Inversion Setup n Monthly fluxes for 1997 to 2001 at 3. 75 o x 5. 0 o resolution (s) n Atmospheric data from NOAA/ESRL cooperative air sampling network (y) n TM 3 gridscale basis functions (H) n Select subset of auxiliary variables (X) n Quantify spatial covariance (Q) n Perform inversion to obtain: ^ q Influence of auxiliary variables on fluxes (β) q Flux best estimates (ŝ) q Estimates of uncertainty for s and β A. M. Michalak (amichala@umich. edu)

Final Set of Auxiliary Variables Combined physical understanding with results of VRT to choose final set of auxiliary variables: • GDP Density • Leaf Area Index (LAI) • Fraction of photosynthetically active radiation (FPAR) • Percent forest / shrub • Precipitation Variable GDP LAI FPAR Precip. F/S b 1. 1 -5. 4 4. 7 0. 6 -0. 6 b + 2 s 0. 6 -6. 4 3. 7 0. 2 -1. 0 b - 2 s 1. 5 -4. 4 5. 7 1. 0 -0. 2 |b/s| 4. 9 10. 8 9. 1 2. 8 3. 0 A. M. Michalak (amichala@umich. edu)

Building up the Best Estimate A. M. Michalak (amichala@umich. edu)

Location of 22 Transcom Regions A. M. Michalak (amichala@umich. edu)

Conclusions - Methodology n Geostatistical inverse modeling avoids the use of prior flux estimates n Covariance structure of flux residuals and model-data mismatch can be quantified using atmospheric data n Benefit of auxiliary data can be quantified n Fluxes and the influence of auxiliary data are estimated concurrently (w/ uncertainties) n Approaches maximizes the use of information while minimizing assumptions n Geostatistical inverse modeling not constrained by prior estimates q Provides independent validation of bottom-up estimates in wellconstrained regions q Approach well suited to show inter-annual variability q Provides accurate measure of uncertainty A. M. Michalak (amichala@umich. edu)

Key Points n If the parameter(s) that you are modeling exhibits spatial (and/or temporal) autocorrelation, this feature must be taken into account to avoid biased solutions n Spatial (and/or temporal) autocorrelation can be used as a source of information in helping to constrain parameter distributions n The field of geostatistics provides a framework for addressing the above two issues A. M. Michalak (amichala@umich. edu)

Acknowledgments n Collaborators: q Pieter Tans, Adam Hirsch, Lori Bruhwiler, Kevin Schaefer, Wouter Peters, Andy Jacobson NOAA/CMDL q Alanood Alkhaled, Sharon Gourdji, Charles Humphriss, Meng Ying Li, Miranda Malkin, Kim Mueller, and Shahar Shlomi, UM q Bhaswar Sen and Charles Miller, JPL q Kevin Gurney, Purdue U. q Peter Kitanidis, Stanford U. n Funding sources: q q q Elizabeth C. Crosby Research Award University Corporation for Atmospheric Research (UCAR) National Oceanic and Atmospheric Administration (NOAA) National Aeronautic and Space Administration (NASA) and Jet Propulsions Laboratory (JPL) National Science Foundation (NSF) Michigan Space Grant Consortium (MSGC) n Data providers: q q q NOAA / CMDL cooperative air sampling network Seth Olsen (LANL) and Jim Randerson (UCI) Christian Rödenbeck, MPIB Kevin Schaefer, NOAA / ESRL NOAA CDC NASA, EROS USGS, CEISIN, Global Precipitation Climatology Centre, UCAR A. M. Michalak (amichala@umich. edu)

QUESTIONS? Anna. Michalak@umich. edu http: //www-personal. engin. umich. edu/~amichala/ A. M. Michalak (amichala@umich. edu)