ROMS 4DVAR Data Assimilation - myroms/roms GitHub Wiki

Incremental Analysis Update

Since no dynamical balance constraint is imposed on the analysis increments computed during 4D-Var, the resulting geostrophic adjustment process at the beginning of the assimilation window can lead to the generation of large amplitude inertia-gravity waves, so-called initialization shocks. These waves typically dissipate within 12-24 hours and are not usually a concern for the physical circulations that develop at the end of the assimilation window. However, these inertia-gravity waves can be very detrimental for coupled biogeochemical configurations of ROMS or those that include sediment processes. Therefore, it is desirable to suppress the generation of inertia-gravity waves during the assimilation window to the maximum extent possible.

Initialization shock is a severe issue in numerical weather prediction (NWP), and mitigation attempts have a long history (e.g., Machenhaur,1977). Two common approaches used to minimize the influence of initialization shock in modern NWP systems are the application of a digital filter (Gauthier and ThΓ©paut, 2001) and incremental analysis updates (IAU; Bloom et al., 1996), although under certain conditions, the two approaches are equivalent (Polavarapu et al., 2004).

IAU is straightforward to implement and has been included as an option in ROMS 4D-Var. Following the usual convention, we will denote by $𝒙(𝑑_k)$ the ocean state vector at time $𝑑_k$. Integration of the ROMS nonlinear model can then be represented as $𝒙(𝑑_{k+1}) = β„³(𝒙(𝑑_k))$ where $β„³$ represents the nonlinear model operators. Let $[𝑑_kβˆ’\tau, 𝑑_k]$ represent the 4D-Var data assimilation window of length $\tau$, where $𝑑_k$ is the current time. The usual practice in 4D-Var is to compute the analysis increment $\Delta 𝒙(𝑑_k βˆ’ \tau)$ at the beginning of the analysis window and then integrate $𝒙(𝑑_k βˆ’ \tau) + \Delta 𝒙(𝑑_k βˆ’ \tau)$ forward in time to the end of the assimilation window using the nonlinear model $β„³$. This occurs at the end of each 4D-Var outer loop. However, since no explicit dynamical balance constraints are imposed on $\Delta 𝒙(𝑑_k βˆ’ \tau)$, the ensuing geostrophic adjustment process during the assimilation window $[𝑑_k βˆ’ \tau, 𝑑_k]$ can lead to the generation of inertia-gravity waves as described above. Instead of adding the increment $\Delta 𝒙(𝑑_k βˆ’ \tau)$ to the initial condition $𝒙(𝑑_k βˆ’ \tau)$, the IAU approach imposes the increment as a forcing term in the nonlinear model so that:

$$𝒙(𝑑_{k+1}) = β„³ 𝒙(t_k)) + 𝒇(𝑑)$$

where $𝒇(𝑑) = \Delta 𝒙(𝑑_k βˆ’ \tau)\Delta 𝑑/𝑇$ for $𝑑_k βˆ’ \tau \leq 𝑑 \leq 𝑑_kβˆ’ \tau + 𝑇$ and $𝒇(𝑑) = 0$ for $𝑑 > 𝑑_k βˆ’ \tau + 𝑇$. In this way, the nonlinear model at the end of each outer loop can be gently steered toward the desired analysis increment over some time interval $[𝑑_k βˆ’ \tau, 𝑑_k βˆ’ \tau + 𝑇]$ during the assimilation window with the result that the initialization shock is greatly reduced. The size of the forcing term $𝒇(𝑑)$ each timestep $\Delta 𝑑$ depends on the choice of $𝑇$, which in turn influences the filtering properties of the IAU procedure. The IAU procedure, as implemented in ROMS, is illustrated in the figure below. Since the forcing amplitude is generally small, the response function resulting from IAU can be conveniently quantified using the tangent linear approximation described by Bloom et al. (1996). The connection between the IAU and the digital filter approach was explored by Polavarapu et al. (2004).

cropped

Implementation:

In ROMS 4D-Var, the IAU is activated by the parameter timeIAU in s4dvar.in, where time_IAU represents $T$ in the equation for $𝒇(𝑑)$ introduced in above. If timeIAU=0, then the IAU is not applied, and the increment $\Delta 𝒙(𝑑_k βˆ’ \tau)$ is added to the initial condition at the beginning of the analysis window as normal.

[!NOTE]

  • The IAU is only available for dual 4D-Var using the RBL4DVAR and SPLIT_RBL4DVAR algorithms. IAU is not currently implemented in the primal formulation of 4D-Var (I4DVAR).
  • If timeIAU $\neq$ 0, the IAU is activated in each outer-loop.
  • If running SPLIT_RBL4DVAR, the analysis phase can be run separately with different values of timeIAU but only if timeIAU was initially chosen as non-zero in the background and increment phases. This allows users to explore different choices of timeIAU without rerunning the inner loops.
  • Thanks to Andy Moore for implementing this option and writing the documentation. Notice that the PDF file uses timIAU instead of timeIAU in the code. ROMS_IAU_v2.pdf

References:

  • Bloom, S.C., L.L. Takacs, A.M. da Silva and D. Ledvina, 1996: Data assimilation using incremental analysis updates, Mon. Wea. Rev., 124, 1256-1271, doi: 10.1175/1520-0493(1996)124<1256:DAUIAU>2.0.CO;2.
  • Gauthier, P., and J.-N. ThΓ©paut, 2001: Impact of the digital filter as a weak constraint in the preoperational 4DVAR assimilation system of MΓ©tΓ©o-France, Mon. Wea. Rev., 129, 2089–2102, doi: 10.1175/1520-0493(2001)129<2089:IOTDFA>2.0.CO;2.
  • Machenhaur, B., 1977: On the dynamics of gravity oscillations in a shallow water model with application to non-linear initialization, Beit. Phys. Atmos., 50, 253-271.
  • Polavarparu, S., R. Shuzan, A.M. Clayton, D. Sankey and Y. Rochon, 2004: On the relationship between incremental analysis updating and incremental digital filtering, Mon. Wea. Rev., 132, 2495-2502, doi: 10.1175/1520-0493(2004)132<2495:OTRBIA>2.0.CO;2.

Modeling Background Standard Deviation

In 4D-Var, the background (prior) error covariance, B, is a large matrix that cannot be computed or stored directly. Still, its effects can be modeled using spatial correlations, C, and spatial convolutions via diffusion operators. To convert correlations into error covariance, we multiply by a diagonal matrix of the background error standard deviations, $\bf\Sigma$. Hence,

$$ \mathbf{B = K \Sigma C \Sigma^{T} K^{T}}$$

Here, K is a balance operator if BALANCE_OPERATOR is activated in ROMS. It allows the information on unobserved state variables to be extracted from directly observed quantities by imposing linear balance relationships between temperature and other state variables using T-S empirical formulas, the linear equation of state, hydrostatic balance, and geostrophic balance. Please check Moore et al. (2011) for more information.

This update includes an approach to computing the standard deviation, $\bf\Sigma$, directly from the background state field as an alternative to climatological values read from the input NetCDF files. It follows the work of Mogensen et al. (2012) by assuming the background errors are proportional to the vertical derivatives of the background field. The background field's practical standard deviation shape has a vertical profile similar to its vertical derivative. For example, the difference between background (prior) temperature and actual temperature is due to a vertical displacement.

In the past, the standard deviation was read from input NetCDF files and computed from a long simulation of the ROMS application. The climatological $\bf\Sigma$ can be categorized as monthly, seasonal, or annual values.

The STD_MODEL option activates the modeling of the standard deviation from the background (prior), which is done in the new module background_std.F. If you want more information about this capability in ROMS, please check Moore et al. (2020).

New parameters are added to the 4D-Var input script s4dvar.in to constrain the standard deviation profile at the mixed-layer depth and deep ocean:

! Modeled standard deviation (STD) of Background Error Covariance parameters.
!
! The Mogensen et al. (2012) formulation assumes that the background errors
! are proportional to vertical derivatives of the state vector field. Its
! error has a similar field profile shape, but the difference with its
! ture error value is due to a vertical displacement.
!
! If COMPUTE_MLD is activated, the mixed-layer depth is computed using the
! criterion from Kara et al. (2000). Otherwise, it will be set to a uniform
! value provided below.

Sigma_max(isFsur) == 0.025d0           ! free surface  maximum STD value

Sigma_max(isUvel) == 0.06d0            ! U-velocity maximum STD value
 Sigma_ml(isUvel) == 0.05d0            ! U-velocity minimum STD at mixed layer
 Sigma_do(isUvel) == 0.02d0            ! U-velocity minimum STD in deep ocean
 Sigma_dz(isUvel) == 500.0d0           ! U-velocity vertical displacement

Sigma_max(isVvel) == 0.06d0            ! V-velocity maximum STD
 Sigma_ml(isVvel) == 0.05d0            ! V-velocity minimum STD at mixed layer
 Sigma_do(isVvel) == 0.02d0            ! V-velocity minimum STD in deep ocean
 Sigma_dz(isVvel) == 500.0d0           ! V-velocity vertical displacement

Sigma_max(isTvar) == 0.33d0   0.056d0  ! 1:NT tracers maximum STD
 Sigma_ml(isTvar) == 0.05d0   0.05d0   ! 1:NT tracers minimum STD at the mixed layer
 Sigma_do(isTvar) == 0.02d0   0.0028d0 ! 1:NT tracers minimum STD in the deep ocean
 Sigma_dz(isTvar) == 40.0d0   40.0d0   ! 1:NT tracer vertical displacement

      mld_uniform == -75.0d0           ! Uniform mixed layer depth value

Notice we have the option COMPUTE_MLD to compute the mixed-layer depth using the approach of Kara et al. (2000) or set a constant value with the mld_uniform input parameter.

Therefore, the difference between the background values of a 3D state variable $\phi_{b}$ and the true value $\phi_{t}$ is due to a vertical displacement $\delta{z}$ of the profile. Thus, to a first-order approximation $$\phi_{t} \approx \phi_{b}(z+\delta{z}) + (\partial{\phi_{b}}/\partial{z}) $$

The error in $\phi_{b}$ is approximated by $(\partial{\phi_{b}}/\partial{z})\delta{z}$, and it is assumed that the standard deviation $\bf\Sigma$ is approximated by $\sigma_{\phi} \sim \mid(\partial{\phi_{b}}/\partial{z})\delta{z} \mid$, with the folllowing constraints:

$$ \sigma_{\phi} = \begin{cases} \hbox{\bf max}(\widehat{\sigma_{\phi}}, \enspace\sigma_{\phi}^{ml}), &\text{if} \quad z \geq -D_{ml} \ \hbox{\bf max}(\widehat{\sigma_{\phi}}, \enspace\sigma_{\phi}^{do}), &\text{if}\quad z < -D_{ml}\end{cases} $$

where

$$ \widehat{\sigma_{\phi}} = \hbox{\bf min}(\mid(\partial{\phi_{b}}/\partial{z})\mid, \enspace\sigma_{\phi}^{max}) $$

Here, $\sigma_{\phi}^{max}$ is the maximum value of $\sigma_{\phi}$, prescribed as Sigma_max(:), while $\sigma_{\phi}^{ml}$ and $\sigma_{\phi}^{do}$ are the minimum values in the mixed layer and deep ocean, prescribed as Sigma_ml(:) and Sigma_do(:), respectively. The vertical displacemen $\delta{z}$ is specified in Sigma_dz(:).

Two new routines are added, def_std.F and wrt_std.F, to write the standard deviation modeled from the background into an output NetCDF file. It will be needed in the split 4D-Var algorithms and for postprocessing. The new filename is also specified in s4dvar.in:

! If computing the standard deviation from the background (prior) state
! vector as an alternative to climatological values read from the
! input NetCDF file, enter the output standard deviation file name,
! [1:Ngrids].

       STDnameC == roms_std_c.nc

We foresee enhancing this capability in the future. We only model the standard deviation for adjusting the initial state vector (zeta, u, v, T, and S). We don't have options for background error on the model (weak constraint), lateral model boundary conditions, surface tracer fluxes, and surface momentum stress. Thus, those standard deviation values are still read from input files.

Notice that the standard deviation structure for I/O management increased its inner dimension from 4 to 5 in mod_iounits.F:

       IF (.not.allocated(STD)) THEN
         allocate ( STD(5,Ngrids) )
       END IF

References:

  • Kara A., P. Rochford, and E. Hulburt, 2000: An optimal definition for ocean mixed layer depth, ''J. Geophys. Res.'', 105, NoC7, pp 16, 803-16, 821, doi: 10.1029/2000JC900072.

  • Mogensen, K., M.A. Balmaseda, and A.T. Weaver, 2012: The NEMOVAR ocean data assimilation system implemented in the ECMWF ocean analysis for system 4, ''ECMWF Tech. Memorandum 668'', 59.

  • Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, and J. Zavala-Garay, 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems, Part I - System overview and formulation, ''Prog. Oceanogr.'', 91, 34-49, https://[doi: 10.1016/j.pocean.2011.05.004](https://www.sciencedirect.com/science/article/abs/pii/S0079661111000516).

  • Moore, A., J. Zavala-Garay, H.G. Arango, C.A. Edwards, J. Anderson, and T. Hoar, 2020: Regional and basin scale applications of ensemble adjustment Kalman filter and 4D-Var ocean data assimilation systems, ''Progress in Oceanography'', 189, 102450, doi: 10.1016/j.pocean.2020.102450.

Split Mixed-Resolution 4D-Var

The split 4D-Var scheme allows mixed-resolution applications in which the background and analysis phases are run at a higher resolution (outer loops) than the increment phase doing the minimization (inner loops) with the tangent linear model (TLM) and adjoint model (ADM) kernels.

  • The coarse grid trajectory used to linearize the TLM and ADM kernels in the inner loops is extracted by a factor of two decimation using the CPP option GRID_EXTRACT when running the higher-resolution background phase.

  • The mixed resolution algorithm improves this application's computational efficiency by over 88 percent compared to the 3km non-splitted 4D-Var in double precision.

  • In some cases, mixed-precision outer loops (double) and inner loops (single) are possible but not recommended because they affect the stability of the tangent linear and adjoint trajectories. It improves the efficiency by an additional 3 percent. However, stability of the solution takes precedence.

  • Both input and output files are processed with the PIO-NetCDF library, accelerating the computations.

  • Grid decimation is only possible if the parent grid size satisfies MOD(Lm+1, 2) = 0 and MOD(Mm+1, 2) = 0. Please check its documentation for more information. Currently, we only support ExtractFlag=2 for decimation in the mixed-resolution split 4D-Var scheme because the land/sea masking complicates extraction at larger factors.

  • The coarse-to-fine increments interpolation is done level-by-level in new ROMS modules roms_interp.F and state_regrid.F, which uses Fortran 2003 CLASS objects. It is very efficient in parallel.

  • A new CPP option SPLIT_EXECUTABLE is added to avoid allocation of the control vectors and adjoint state arrays to save computer memory in the computation of the high-resolution background and analysis phases, which may limit running larger applications because they do not fit into computer memory.

  • A new generic script, script submit_split_rbl4dvar.sh, is provided to run the split 4D-Var data assimilation system. Configuring the system is tricky; the devil is in the details, which require full attention! It is recommended for users with experience in 4D-Var and not for novice users learning variational data assimilation.

  • The mixed precision requires two ROMS executables: romsM_nl (outer loops) and romsM_da (inner loops). Special build scripts, build_splt.csh and build_split.sh are provided to generate the excutables:

    • To create the romsM_nl executable, use: build_split.sh -nl -pio -j 10
    • To create the romsM_da executable, use: build_split.sh -da -pio -j 10

  • The ROMS standard input files roms_nl.in and roms_da.in are created from the template scripts.

  • The split sheme facilitates running ROMS 4D-Var data assimilation cycles more efficiently. Here is an efficiency table from a desktop Linux box (16 CPUs and 1 GPU with 384 CUDA cores, NVIDIA) on 12 CPUs with a 3x4 partition. I compiled with ifort using Spack-Stack 1.9.

  • A test case is provided for a regional application of the U.S. East Coast (USEC), where the outer loops are run at 3km grid resolution (561x243x50) resolution (CFF-USEC3). In contrast, the inner loops (Increment phase) are run at a coarser 6km grid (281x122x50) resolution to accelerate the computations. For detailed information and Instructions about the mixed-resolution 4D-Var test, please check the provided test case:

    https://github.com/myroms/roms_test/blob/main/USEC/RBL4DVAR_mixres/Readme.md

    The USEC test case shows that 3km and 6km resolution increments are indistinguishable. The figures below show the increments for the 4D-Var Cycle 1, Aug 27 - Aug 30, 2019, ERA-5 forcing case. From top to bottom, there are 3km and 6km increments for free surface, potential temperature, salinity, u-velocity, and v-velocity at 20m depth.

3km Increments at z=20m 6km Increments at z=20m

Similarly, the 3km and 6km increments for free surface, potential temperature, salinity, u-velocity, and v-velocity at 50m depth are indistinguishable, illustrating the robustness of the mixed-resolution 4D-Var in the USEC application in a highly dynamically active region that includes the Gulf Stream.

3km Increments at z=50m 6km Increments at z=50m