Martian Data Assimilation Using the LETKF

1
Ensemble Data Assimilation and Breeding in the
Ocean, Chesapeake Bay, and Mars
Matthew J. Hoffman Dissertation Defense May 15,
2009 University of Maryland, College
Park Advisors Eugenia Kalnay and James A.
Carton Committee Raghu Murtugudde, Brian Hunt,
Kayo Ide, Michael Evans
2
Outline

• Ocean Instabilities Captured by Breeding
• The Local Ensemble Transform Kalman Filter
  (LETKF)
• Chesapeake Bay Data Assimilation
• Data Assimilation of the Martian Atmosphere

3
Ocean Instabilities Captured by Breeding

• Overview of ocean instabilities
• Overview of the breeding method
• Application to global ocean model
• Development of bred vector energy equations to
  diagnose instability dynamics
• Study of Pacific tropical instabilities

Courtesy NASA Earth Observatory
4
Ocean Instabilities
Ducet et al., 2000

• Flow Instabilities (eddy kinetic energy) are
  prevalent in the upper ocean
• Most occur in strong currentswestern boundary
  currents, Southern Ocean
• Instabilities take place on different timescales
• Tropical Pacific instabilities are some of the
  strongest

5
Pacific Tropical Instabilities
Jesse Allen, NASA Earth Observatory SST from
Advanced Microwave Scattering Radiometer on Aqua

• Pacific Tropical Instability Waves are seen in
  the Pacific equatorial cold tongue
• Periods of 20-30 days, Wavelength of 1000km
• Tropical waves are coupled to the atmospheric
  boundary layer and are important for heat and
  momentum balances
• Masina et al. (1999) argued for baroclinic energy
  conversion dominating
• Qiao and Weisberg (1995)argued for barotropic
  energy conversion dominating

6
Overview of Breeding Method

• Developed by Toth and Kalnay (1993, 1997) to
  estimate the shape of growing errors in a
  non-linear atmospheric model
• Also provides initial conditions for ensemble
  forecasting
• 2 parameters in the methodrescaling size and
  time between rescaling
• Parameters can be tuned to isolate instabilities
  of different time scales
• Yang et al. (2005) used breeding on a coupled GCM
  to identify slow growing ENSO modes

7
Overview of Breeding Method

• A small, random perturbation is added to the
  initial state of the system
• Both the perturbed and unperturbed (control)
  conditions are integrated forward in time
• The control forecast is subtracted from the
  perturbed forecast, yielding the bred vector
• The bred vector is rescaled to its initial size
  and added to the control forecast as a new
  perturbation

8
MOM2 Global Model

• GFDL Modular Ocean Model (MOM) 2b code
• Driven by NCEP reanalysis winds from 1950-1995
• SST and surface salinity from World Ocean Atlas
  1994
• Same setup used by Carton et al. (2000) for SODA
• 1 resolution in longitude, stretched latitude
  grid ranging from 1 in midlatitudes to ½ in
  tropics
• 20 vertical levels 15 meters level thickness
  near surface

9
Bred Vectors

• Bred vectors from1987-1989
• 10 day bred vectors identify many instabilities
  in the ocean
• Instabilities are seen in the Southern Ocean, in
  boundary currents, and in the Tropical Pacific,
  among other locations

10
Tropical Pacific Bred Vectors

• Fall/Winter 1988
• Meridional velocity (shaded) and meridional
  velocity bred vector (contour)
• Bred vectors capture the growing instabilities

11
Tropical Pacific Bred Vectors

• Seasonal cycle is clear
• Speed is 0.46m/s
• Wavelength is 1000km
• 25 day period

12
Tropical Pacific Bred Vectors

• Seasonal cycle is clear
• Speed is 0.46m/s
• Wavelength is 1000km
• 25 day period
• Interannual cycle tied to El Niño-La Niña cycle
  (ENSO)

El Niño
13
Tropical Pacific Bred Vectors

• Seasonal cycle is clear
• Speed is 0.46m/s
• Wavelength is 1000km
• 25 day period
• Interannual cycle tied to El Niño-La Niña cycle
  (ENSO)

La Niña
14
Bred Vector Kinetic Energy

• Momentum Equations
• Kinetic energy defined as
• Bred kinetic energy is

15
Bred Vector Kinetic Energy

• Terms have physical interpretation

16
Bred Vector Kinetic Energy

• Terms have physical interpretation
• Horizontal and vertical divergence of energy
  transport

17
Bred Vector Kinetic Energy

• Terms have physical interpretation
• Horizontal and vertical divergence of energy
  transport
• Work of pressure force

18
Bred Vector Kinetic Energy

• Terms have physical interpretation
• Horizontal and vertical divergence of energy
  transport
• Work of pressure force
• Baroclinic conversion term

19
Bred Vector Kinetic Energy

• Terms have physical interpretation
• Horizontal and vertical divergence of energy
  transport
• Work of pressure force
• Baroclinic conversion term
• Barotropic conversion term

20
Bred Vector Kinetic Energy

• Terms have physical interpretation
• Horizontal and vertical divergence of energy
  transport
• Work of pressure force
• Baroclinic conversion term
• Barotropic conversion term
• Dissipation term

21
Bred Vector Kinetic Energy

• Tropical Pacific shows positive conversion (bred
  potential to bred kinetic)
• Shows Instability Growth
• South Atlantic shows negative conversion (bred
  kinetic to bred potential)
• Stabilizing region

22
Bred Vector Kinetic Energy

• Tropical Pacific shows positive conversion (bred
  potential to bred kinetic)
• Shows Instability Growth
• South Atlantic shows negative conversion (bred
  kinetic to bred potential)
• Stabilizing region

23
Bred Vector Kinetic Energy

• Tropical Pacific shows positive conversion (bred
  potential to bred kinetic)
• Shows Instability Growth
• South Atlantic shows negative conversion (bred
  kinetic to bred potential)
• Stabilizing region

24
Tropical Pacific Instabilities

• Monthly averages over a 30 year period are shown
  for October
• Depth averaged over upper 150m
• Baroclinic conversion is strongest from 3N-5N
• Barotropic conversion is strongest at the equator
• Baroclinic conversion is stronger in this model
• Energy conversion is strongest when bred vectors
  are strongest (La Niña)

25
Tropical Pacific Instabilities
At 3.5N
At 0.25N

• Baroclinic conversion is strongest at coldest
  temperatures (cold tongue)
• Barotropic conversion is strongest at shear
  points (between South Equatorial Current and
  Equatorial Undercurrent)
• Different locations for the different mechanisms

26
Summary

• Breeding is an easy way to identify instabilities
  in a dynamical system
• Breeding energy equations allow bred vectors to
  be used to diagnose the dynamical causes of
  instabilities
• Tropical Pacific instabilities have a baroclinic
  and barotropic component
• Baroclinic component is stronger in this model
  and occurs along the north edge of the cold
  tongue between 3N and 5N
• Barotropic component occurs at the equator
  between the South Equatorial Current and the
  Equatorial Undercurrent

27
Outline

• Ocean Instabilities Captured by Breeding
• The Local Ensemble Transform Kalman Filter
  (LETKF)
• Chesapeake Bay Data Assimilation
• Data Assimilation of the Martian Atmosphere

28
Ensemble Kalman Filter

• Data assimilation combines observations with a
  previous state estimate, called the background,
  based on older observations. The resulting state
  estimate is called the analysis
• Here we use the Local Ensemble Transform Kalman
  Filter (LETKF Hunt et al., 2007)
• The analysis comes from minimizing the cost
  function
• The background covariance, Pb, should change each
  assimilation step based on the current errors
• Computationally, it is not possible to fully
  calculate this covariance every step
• Instead, we use an ensemble method to
  characterize the uncertainty

29
An Ensemble Assimilation Cycle
Observations
x
Model Runs
x
x
x
x
To next background
xa
x
x
LETKF
xb
x
x
x
x
x
Background Ensemble
New Analysis Ensemble
Previous Analysis Ensemble

• Analysis is localanalysis is performed
  independently at each grid point using
  observations in a neighborhood of that point
• In practice the covariance is artificially
  inflated to account for underestimation

30
Chesapeake Bay Data Assimilation

• First phase of implementation of the LETKF on the
  Chesapeake
• Short term goal was studying the impact of the
  LETKF and evaluating the current observing
  system
• Long term goal is an assimilation system as part
  of the Chesapeake Bay forecast system
• Identical Twin Experiments with Random and
  Realistic Data Coverage are presented here

NASA/Goddard Space Flight Center Scientific
Visualization Studio
31
Chesapeake Bay

• Largest estuary in N. America
• Over 1 billion brought in yearly by fishing
  industry
• 300km long, 50km at widest
• Average depth of 6.5m
• Deep, narrow channel in the main stem

NASA/Goddard Space Flight Center Scientific
Visualization Studio
32
Chesapeake Bay Circulation

• Salt water enters Bay in deep channel
• Fresh water enters at surface from rivers
• Salinity distribution has wedge shape
• Tidal amplitude is moderaterange is less than 1m

NASA/Goddard Space Flight Center Scientific
Visualization Studio
33
River Discharge
Susquehanna

• 143 million liters of fresh water per minute
  enter the Bay

Patapsco
Chester
Choptank
Patuxent
Nanticoke
Potomac
Pokomoke
Rappahannock
York
James
NASA/Goddard Space Flight Center Scientific
Visualization Studio
34
River Discharge
Susquehanna

• 143 million liters of fresh water per minute
  enter the Bay
• 50 comes from the Susquehanna River

Patapsco
Chester
Choptank
Patuxent
Nanticoke
Potomac
Pokomoke
Rappahannock
York
James
NASA/Goddard Space Flight Center Scientific
Visualization Studio
35
River Discharge
Susquehanna

• 143 million liters of fresh water per minute
  enter the Bay
• 50 comes from the Susquehanna River
• 18 from the Potomac River

Patapsco
Chester
Choptank
Patuxent
Nanticoke
Potomac
Pokomoke
Rappahannock
York
James
NASA/Goddard Space Flight Center Scientific
Visualization Studio
36
River Discharge
Susquehanna

• 143 million liters of fresh water per minute
  enter the Bay
• 50 comes from the Susquehanna River
• 18 from the Potomac River
• 14 from the James River
• Susquehanna discharge determines flow at the head
  on timescales of 5 days

Patapsco
Chester
Choptank
Patuxent
Nanticoke
Potomac
Pokomoke
Rappahannock
York
James
NASA/Goddard Space Flight Center Scientific
Visualization Studio
37
Chesapeake Bay Model

• Numerics are from the Regional Ocean Modeling
  System (ROMS)
• Curvilinear grid with 100x150x10 resolution
• Same bathymetry and forcing as ChesROMS (Xu et
  al., 2009)only difference is vertical resolution
• Terrain following sigma coordinate

38
Open Boundary Forcing

• 9 tidal constituents from ADCIRC model
• Non-tidal water levels are used from NOAA
  National Ocean Service program
• Identical to Li et al., (2005)
• Salinity and temperature are nudged to
  climatology from WOA01
• Waves propagate through the boundary

39
River and Air Surface Forcing

• Daily freshwater discharges are prescribed for 9
  tributaries from USGS data
• Air-surface boundary is set from North American
  Regional Reanalysis (NARR)
• 3-hourly winds
• Net shortwave and downward longwave radiation
• Temperature
• Relative humidity
• Pressure

40
Identical Twin (IT) Experiments

• The model is run and this run is nature
• Random Gaussian errors with a prescribed standard
  deviation are added to the nature run to create
  observations
• Error standard deviations are 0.5C, 0.6psu, and
  0.05m/s
• Assimilation is started on January 10, 1999 in
  the nature run
• Initial ensemble is created by taking random
  states from the spin up
• Model is run from the background of the initial
  ensemble as the free run forecasti.e. the case
  with no data assimilation
• Analyses are performed every 3 hours

41
10 Data Coverage IT Experiments

• Observations are simulated at 10 of the grid
  points
• Points are randomly chosen and change with
  analysis time
• Observations of temperature, salinity, and
  currents
• Temperature, salinity, and current fields are
  assimilated
• 9 inflation, 5 grid point horizontal
  localization radius, 4 level vertical
  localization, 16 member ensemble
• Drop in free run forecast error indicates the
  importance of forcing

42
Dependence on Data Coverage
Free Run
0.1
0.5
10

• Observations at 10 of the grid points is 4000
  obs. per variable
• This is more than can be expected
• The analysis is degraded by few observations, but
  improvement is still seen
• 0.1 observations is 40 obs. per variable, which
  is close to a real average number

43
Real Observational Data

• Buoy observations are available from the
  Chesapeake Bay Program (CBP) and the Chesapeake
  Bay Observing System (CBOS)

44
Real Observational Data

• Buoy observations are available from the
  Chesapeake Bay Program (CBP) and the Chesapeake
  Bay Observing System (CBOS)
• 6 CBOS stations report profiles of temperature
  and salinity

http//hpl.cbos.org/
45
Real Observational Data

• Buoy observations are available from the
  Chesapeake Bay Program (CBP) and the Chesapeake
  Bay Observing System (CBOS)
• 6 CBOS stations report profiles of temperature
  and salinity
• 120 CBP stations report temperature and salinity
  profiles with 40 CBP stations in the main stem of
  the Bay

46
Real Observational Data

• Buoy observations are available from the
  Chesapeake Bay Program (CBP) and the Chesapeake
  Bay Observing System (CBOS)
• 6 CBOS stations report profiles of temperature
  and salinity
• 120 CBP stations report temperature and salinity
  profiles with 40 CBP stations in the main stem of
  the Bay
• CBOS stations report every 6-30 minutes, CBP
  report every 2 weeks-1 month

47
Simulated Real Observations

• Observations are simulated at real location and
  analysis time
• Temperature and salinity observations are
  assimilated
• 16 ensemble members, 2 inflation, 4 level
  vertical localization, varied horizontal
  localization
• Large improvements are seen with influx of CPD
  data

48
Simulated Real Observations
Free Run SST Error
Analysis SST Error

• Day 5 of the simulation
• Analysis error is lower in Bay
• Analysis error is higher in open oceanno
  observations

49
Simulated Real Observations
Analysis SST Error
Temperature Free Run Error
Temperature Analysis Error
Salinity Analysis Error
Salinity Free Run Error

• Day 5 of the simulation

50
Simulated Real Observations
Analysis SST Error
Temperature Free Run Error
Temperature Analysis Error
Salinity Analysis Error
Salinity Free Run Error

• Day 5 of the simulation

51
Simulated Real Observations
Free Run SST Error
Analysis SST Error

• Day 30 of the simulation
• Analysis error is lower in Bay
• Analysis error is higher in open oceanno
  observations

52
Summary

• The Local Ensemble Transform Kalman Filter
  improves state estimates in the Chesapeake Bay in
  identical twin experiments with randomly
  distributed observations
• Using simulated observations at real locations
  and analysis times, the LETKF improves the
  temperature, salinity, and current fields
• Spatially, the temperature analysis improves most
  in the Bay near the surface
• More errors are seen in the open ocean where
  there are no obs
• Lack of observations over long stretches is an
  issue

53
Martian Data Assimilation Using the LETKF

• First phase of the project, which will be
  continued by Steven Greybush for his dissertation
• Goal To create a Martian climate reanalysis
• Using Thermal Emission Spectrograph (TES) derived
  temperature profiles
• NASA/NOAA Global Circulation Model

Courtesy NASA/JPL-Caltech
54
Mars vs. Earth
55
Mars Topography
Vastitas Borealis
Olympus Mons
Tharsis
Hellas Impact Basin
Read and Lewis, 2004
Courtesy NASA/JPL-Caltech
56
Dominant Unstable Modes
Barnes et al., 1993

• Mars climate attractor has a low dimension
• Atmosphere is dominated by low wavenumber global
  baroclinic modes (m1-3) and diurnal tide (Read
  et al., 2006)
• 80 of total energy is in first 7 EOFs
  (Martinez-Alvarado et al., 2008)

57
NASA/NOAA Mars GCM

• Developed by John Wilson at GFDL
• Uses finite volume numerics in dynamical core
• Latitude-longitude grid
• 60x36 grid points
• (6x5.14 resolution)
• 28 vertical levels
• Uniform dust concentration
• Hybrid vertical coordinate
• s-coordinate near surface

Olympus Mons
Vastitas Borealis
Tharsis
Hellas Impact Basin
Courtesy S. Greybush
58
Full Coverage IT Experiments

• Initial experiments used observations at every
  grid point
• Random Gaussian observation errors with a 1K
  standard deviation
• 16 ensemble members were used
• Covariance inflation was 10
• 1200km horizontal localization radius
• Assimilations performed every 6 hours
• Only temperature observations assimilated
• The assimilation was started on day 10 of the
  nature run
• The initial ensemble is created by taking states
  from the previous 10 days

59
Full Coverage IT Experiments

• Initial ensemble averages out the diurnal cycle
• First analysis recaptures the cycle and greatly
  reduces errors

60
Full Coverage IT Experiments

• Analysis errors are below 1K after 1 day
• Analysis and forecast are below free run for
  majority of the run
• Free run forecast improves in time due to
  importance of forcing

61
Full Coverage IT Experiments

• Zonal wind estimate is improved even without
  observations
• Analysis and forecast are far below free run
• Free run forecast improves in time in zonal wind
  as well

62
Full Coverage IT Experiments

• Level 25 is near the surface and is dominated by
  the diurnal cycle
• Level 25 zonal winds are influenced by topography
• Level 5 has a strong zonal jet in the winter
  hemisphere

63
Full Coverage IT Experiments
Level 25 Temperature RMS Error
Level 5 Temperature RMS Error

• Free run forecast in level 25 drops immediately
  this is due to forcing at the surface
• Free run forecast in level 5 varies more free
  atmosphere
• Analysis is better than free run and below
  observation error at both levels

64
TES Observations

• Observations from 1999-2005, 25,000 polar orbits
• IR radiances with 3x3km footprint
• Vertical temperature profiles and dust opacity
  retrievals are derived
• 19 vertical levels in temperature profiles

Mars Global Surveyor
6 hour observation track
Courtesy NASA/JPL-Caltech
Courtesy S. Greybush
65
Simulated TES Observations

• Observations simulated at closest grid point to
  TES track at closest vertical level
• TES corresponds to observations at 14 vertical
  levels
• Vertical localization radius in the LETKF must be
  increased for the simulated TES experiment
• Average number of observations per 6 hours Full
  coverage 57120, TES 4300
• 1200km horizontal localization radius, 16
  ensemble members, 10 inflation used for both
• Larger temperature observation error was used
  (3K)

Level
1
66
Simulated TES Observations
Temperature
Zonal Wind

• Simulated TES observation experiment asymptotes
  at larger values
• More observations lead to quicker error reduction
  in the analysis
• After 12-15 days analysis error approaches
  asymptote
• Both experiments stay below the observation error
  and the free run forecast

67
Errors in the Analysis
Full coverage analysis error (contour) Zonal
mean Temperature (shaded)
Simulated TES observation analysis error
(contour) Zonal mean Temperature (shaded)
Uppermost Observation 0.11mb

• Largest analysis errors in full coverage
  experiment are in Northern Hemisphere temperature
  front
• From baroclinically unstable waves
• Largest analysis errors using simulated TES
  observations are in upper atmosphere and zonal
  jet
• From lack of observations in upper levels

68
Summary

• The Local Ensemble Transform Kalman Filter
  improves temperature and wind state estimates in
  the MGCM using temperature observations
• In both full coverage and simulated TES
  observation experiments the analysis error stays
  far below the free run forecast
• Largest errors in analysis using full coverage
  observations are from baroclinic waves along the
  Northern Hemisphere temperature front
• Largest errors in analysis using simulated TES
  observations are in the upper atmosphere in the
  zonal jet
• Due to lack of TES observations in above 0.11mb

69
Future Plans

• Complete study of South Atlantic instability
• Chesapeake assimilation in the presence of
  forcing errors
• Real observations in Chesapeake Bay
• Improved localization in both Mars and Chesapeake
  Bay systems
• Effect of assimilation on tracer states (oxygen)

70

• END