4D Variational Data Assimilation Observation Operators

1. 4D Variational Data Assimilation Observation Operators 4D Variational Data Assimilation Observation Operators Hernan G. Arango

2. ROMS 4DVAR ALGORITHMS • Strong Constraint • Conventional (S4DVAR): outer loop,NLM,ADM • Incremental (IS4DVAR): inner and outer loops,NLM,TLM,ADM(Courtier et al., 1994) • IS4DVAR_OLD: inefficient old conjugate gradient algorithm (is4dvar_ocean_old.h, descent.F) • IS4DVAR: new conjugate gradient algorithm, CONGRAD, Fisher 1997 (is4dvar_ocean.h, cgradient.h) • IS4DVAR, LANCZOS: conjugate gradient and Lanczos algorithm, Fisher 1997 (is4dvar_ocean_lanczos.h, cgradient_lanczos.h) • Weak Constraint • Indirect Representer Method (W4DVAR): inner andouter loops,NLM,TLM,RPM,ADM (Egbert et al., 1994; Bennett et al, 1997) • Physical Space Statistical Analysis (W4DPSAS):inner and outer loops,NLM,TLM,ADM (Courtier, 1997)

3. xk = B-1/2(xk + xk-1 – xb) xk = B1/2vk + xk-1 – xb Strong Constraint, Incremental 4DVAR Let’s introduce a new minimization variablev, such that: yielding J(vk) = ½(vk)Tvk + ½(Hxk – dk-1)TO-1(Hxk – dk-1) The gradient ofJin minimization-space, denotedv J, is given by: v J = vk + BT/2HTO-1(Hxk – dk-1) = vk + BT/2x Jo=>vk + W-1/2LT/2GS The background-error covariance matrix can be factored as: B = SCS => S(GL1/2W-1/2)(W-1/2LT/2G)S whereSis the background-error standard deviations,Cis thebackground-error correlations which can be factorized asC = C1/2CT/2,Gis the normalization matrix which ensures that the diagonalelements ofCare equal to unity,Lis a 3D self-adjoint filteringoperator, andWis the grid cell area or volume.

4. Basic IS4DVAR Procedure • Choose an x(0) = xb(0) • Integrate NLROMS t  [0, ] and save x(t) (NLM at OBS) (a) Choose a x(0) (b) Integrate TLROMS t  [0, ] and compute J (TLM at OBS) (c) Integrate ADROMS t  [0, ] to yield (ADM forcing at OBS) (d) Compute (e) Use a descent algorithm to determine a “down gradient” correction to x(0) that will yield a smaller value of J (f) Back to (b) until converged (3) Compute new x(0) = x(0) + x(0) and back to (2) until converged Outer Loop Inner Loop

5. Incremental, Strong Constraint 4DVar (IS4DVAR) • Given a first guess (a forward trajectory)… • And given the available data…

6. Incremental, Strong Constraint 4DVar (IS4DVAR) • Given a first guess (a forward trajectory)… • And given the available data… • IS4DVAR computes the changes (or increments) to the initial conditions so that the forward model fits the observations.

7. 4DVAR Observations NetCDF File • Utility/obs_initial.F • Utility/obs_read.F • Utility/obs_write.F • Utility/obs_scale.F • Utility/obs_depth.F • Utility/extract_obs.F • Adjoint/ad_extract_obs.F • Adjoint/ad_misfit.F

8. Metadata Dimensions: survey Number of different time weight Number of interpolation weight datum Observations counter, unlimited dimension Variables: Nobs(survey) Number of observations per time survey survey_time(survey) Survey time (days) obs_type(datum) State variable ID associated with observation obs_time(datum) Time of observation (days) obs_lon(datum) Longitude of observation (degrees_east) obs_lat(datum) Latitude of observation (degrees_north) obs_depth(datum) Depth of observation (meters or level) obs_Xgrid(datum) X-grid observation location (nondimensional) obs_Ygrid(datum) Y-grid observation location (nondimensional) obs_Zgrid(datum) Z-grid observation location (nondimensional) obs_error(datum) Observation error, assigned weight obs_value(datum) Observation value

9. Observations NetCDF 8 7 5 6 4 3 1 2 (i2,j2,k2) (i1,j1,k1) dimensions: survey= 1 ; tate_variable= 7 ; datum = UNLIMITED ; // (79416 currently) variables: charspherical; spherical:long_name = "grid type logical switch" ; intNobs(survey) ; Nobs:long_name = "number of observations with the same survey time" ; doublesurvey_time(survey) ; survey_time:long_name = "survey time" ; survey_time:units = "days since 2000-01-01 00:00:00" ; survey_time:calendar = "365.25 days in every year" ; doubleobs_variance(state_variable) ; obs_variance:long_name = "global (time and space) observation variance" ; obs_variance:units = "squared state variable units" ; intobs_type(datum) ; obs_type:long_name = "model state variable associated with observation" ; obs_type:units = "nondimensional" ; doubleobs_time(datum) ; obs_time:long_name = "time of observation" ; obs_time:units = "days since 2000-01-01 00:00:00" ; obs_time:calendar = "365.25 days in every year" ; doubleobs_depth(datum) ; obs_depth:long_name = "depth of observation" ; obs_depth:units = "meter" ; doubleobs_Xgrid(datum) ; obs_Xgrid:long_name = "x-grid observation location" ; obs_Xgrid:units = "nondimensional" ; doubleobs_Ygrid(datum) ; obs_Ygrid:long_name = "y-grid observation location" ; obs_Ygrid:units = "nondimensional" ; doubleobs_Zgrid(datum) ; obs_Zgrid:long_name = "z-grid observation location" ; obs_Zgrid:units = "nondimensional" ; doubleobs_error(datum) ; obs_error:long_name = "observation error, assigned weight, inverse variance" ; obs_error:units = "inverse squared state variable units" ; doubleobs_value(datum) ; obs_value:long_name = "observation value" ; obs_value:units = "state variable units" ;

10. Processing • Use hindices, try_range and inside routines to transform (lon,lat) to (,) • Find how many survey times occur within the data set (survey dimension) • Count observations available per survey (Nobs) and assign their times (survey_time) • Sort the observation in ascending time order and observation time for efficiency • Save a copy of the observation file • Several matlab scripts to process observations

11. ROMS GRID • Horizontal curvilinear orthogonal coordinates on an Arakawa C-grid • Terrain-following coordinates on a staggered vertical grid

12. Curvilinear Transformation

13. Staggered C-Grid, RHO-points

14. Staggered C-Grid, U-points

15. Staggered C-Grid, V-points

16. Vieste (Italy) Dubrovnik (Croatia) Longitude Vertical Terrain-following Coordinates Depth (m)

17. Parabolic Splines Reconstruction

18. PARALLEL TILE PARTITIONS 8 x 8 Ny } } Nx

19. East-West MPI Communications Nonlinear With Respect To Tile R With Respect To Tile R Adjoint

20. North-South MPI Communications With Respect to Tile B Nonlinear Adjoint

21. Observations

22. Longterm Ecosystem Observatory New Jersey Shelf Observing System 30km x 30km 1998-2001 Satellites, Aircraft, Surface RADAR, Glider AUVs 300km x 300km Beginning 2001 LEO-15 LEO NJSOS THE STATE UNIVERSITY OF NEW JERSEY RUTGERS Station Field 3km x 3km 1996-Present

23. Assumptions • All scalar observations are assumed to be at RHO-points. • All vector observations are assumed to be rotated to curvilinear grid, if applicable. Vector observations are always measured at the same location. • All observation horizontal locations are assumed to be in fractional curvilinear grid coordinates. • Vertical locations can be in fractional levels (1:N) or actual depths (negative values). • Removal of tidal signal? • Filtering of non-resolved processes?

24. Observation Operators • Currently, all observations must be in terms of model state variables (same units): • 2D configuration: zeta, ubar, vbar • 3D configuration: zeta, u, v, T, S, … • There is not a time interpolation of model solution at observation times: time - 0.5*dt < ObsTime < time + 0.5*dt • There is not observations quality control (screening) inside ROMS, except ObsScale. • No observation constraints are implemented (Satellite SST measurements)

25. Observation Interpolation • Only spatial linear interpolation is coded. • If land/sea masking, the interpolation coefficients are weighted by the mask. • If shallower than z_r(:,:,N), observations are assigned to the surface level. • If deeper than z_r(:,:,1), observations are assigned to bottom level.

26. Recommedations • Create a NetCDF file for each observation type. • Use a processing program to meld NetCDF observation files (nc_4dvar_meld.m). • Keep a master copy of each observation file, in case that you are running your application at different resolutions. • Decimation of observations. Finite volume representation?

27. BACKGOUND ERROR COVARIANCE

28. Model/Background Error Covariance, B • Use a generalized diffusion squared-root operator (symmetric) as in Weaver et al. (2003): B = S C S = S (G L1/2 W-1/2) (W-1/2 LT/2 G) • The normalization matrix,G, ensure that the diagonalelements of the correlation matrix,C, are equal to unity.They are computed using theexact(expensive) orrandomization(cheaper) methods. • The spatial convolution of the self-adjoint filtering operator,L1/2, is split in horizontal and vertical components anddiscretized bothexplicitlyand implicitly. • The model/background standard deviation matrix,S, iscomputed from long (monthly, seasonal) simulations. • The grid cell area or volume matrix,W-1/2, is assumed to betime invariant.

29. Horizontal Vertical (implicit) Vdecay = 100 m Hdecay = 100 km Model/Background Error Correlation (C)

30. SSH Temperature EAC EAC Bottom Level Model/Background Error Correlation Normalization Coefficients (G)

31. job_is4dvar.sh

32. build.sh

33. s4dvar.in