Assimilating SSH – high res model blowup! Stability issue?

[lmp4]

Having previously successfully assimilated SSTs, last week I moved on to starting to do some test SSH’s assimilation runs for the ROMS I4DVAR system for a 3km model of the Angola/Congo region.

To summarize I have been having a major issue relating to assimilating SSHs into my model;

• Model blows up on initialization of final nonlinear run (after the inner loops) after time step 0.

My observation data is AVISO SSH data, obtained using the matlab script d_ssh_obs.m (located in the matlab trunk of subversion). This script actually required a few simple adjustments for current use (as of 2015) such as datum changes and the conversion to meters which is now not needed. I have performed the spatial and temporal offset so the observations and ICs numbers wise are not far off. Note I have yet to remove the steric signal since these were simply testing the system.

Now for my assimilation here are my observations (top) and ICs (below);

ObsJul06.png (18.52 KiB) Viewed 5505 times

HycomJul06.png (19 KiB) Viewed 5505 times

Now these I am aware that the initial conditions (obtained from hycom climatology) are quite different (especially in the north of my domain) to the observations. Again I have assumed this shouldn’t affect the assimilation process too much since the numbers are still in the same order at least.
I am running the test assimilation with 5 inner loops for speed but at 30 seconds time step to keep the max bara courent number down to around ~0.2 (when assimilating SSTs I found the model kept blowing up with a 60 second timestep, ~0.4 max bara courent number).
However upon the end of the final loop and the forward analysis run; the model blows up after one timestep…. And upon investing the rst file seems to blow up at one point in V at the bottom of the domain… note that none of my observations are rejected.

Rst_SSH_obs_overlayed.png (8.21 KiB) Viewed 5505 times

Here the position of where v (red) has blown up is shown with the observation points superimposed (light orange);
I did this plot since I originally thought maybe a specific observation could be causing it and wanted to get a reference to whether the blow-up and individual observation overlap; however to test this theory anyways I did also try shrinking my observations to a small square in the domain away from the site of blow up and it still blows up in the same place…

Note there is no island at this area and my rxo and rx1 values are good (see log error file), so I’m quite sure it’s not a bathymetry issue. I have also have checked normalization coefficients and they seem good, although I have kept the same correlation length scales for my 3km model as my 5km model, perhaps this could be an issue?

Interestingly this is the same place my model was blowing up when I assimilating SSTs with a 60 second time step, so I thought perhaps it may be the adjoint model requires an even smaller time step? So I tired 20 seconds, max bara courent number ~0.15 and it still blows up. I am sure that a courent number of 0.15 is low enough for the model right?

I have also tried different ICs, with a different month (so different norm coeffs and std dev files) but I still get the same result.

Interesting since I also have a 5km model of the same domain, I tested the SSH assimilation and it worked!

Ok so to sum up this is a brief list of everything I have tried;

1. Made my observations cover a smaller area in my domain;
2. Tried several different ICs, spanning over different months;
3. Changed the timestep from 30 seconds to 20 seconds;
4. Tried to assimilate the observation with my lower resolution 5km model and it worked!

As previously mentioned I am starting to run out of ideas now
although what I havn't tried is;

1.Changing the correlation length scales;
2.Changing ndtfast as well as timestep to even smaller, say 10 seconds;

Since the model blows up with SSTs, with a 60 second timestep in the same area as it is blowing now assimilating SSH with a 30 second timestep, I can only guess that’s it’s a stability issue within my model .

I have been running with a 30.000 timestep size (s) for 3-D equations and 15 ndtfast split, is this an unstable setup, but works for assimilating SSTs?

Why the 5km model can assimilate and the 3km model cannot is strange? The CCP setting in both models are the same so perhaps the change in resolution is too much for assimilation process?

Anyways I’ve attached my .h file, .in, and log error file.

Any suggestions would be greatly appreciated!

[lmp4]

*UPDATE*

It works! Changed the time-step to 10 seconds...

So that's now;

Code: Select all

 Minimum barotropic Courant Number =  2.16784199E-03
 Maximum barotropic Courant Number =  7.56383597E-02
 Maximum Coriolis   Courant Number =  5.22652060E-04

So since running the assimilation driver takes up alot of computational power anyways any advise on how to stabilize this model while increasing the time-step?

I can imagine a 30 inner loop assimilation run over three days would take a very long time with a 10 second time-step :/

[arango]

I hope that you are plotting the cost function to see if 4D-Var is converging. It is not clear to me if you are using 5 or 30 inner loops. If you are using 5 inner loops, I doubt that the 4D-Var algorithm converged and this may explain your blow-up when running the nonlinear model.

If your prior SSH is far away from the AVISO data, the increment is big and may generate surface gravity and internal waves that may shock the system and this may explain why you need to lower the time step for stability. This is all due initial adjustment of posterior fields. Maybe when you have a better background, in latter cycles, the adjustment is very minimal and you can increase the time-step. You can also play with the error hypothesis for SSH and have you estimate the standard deviation for the background error covariance.

[lmp4]

Hi Arango, Thank you for your advice and sorry for the late response.

I was using 5 loops since I wanted a quick test for the system, however after plotting the cost function you are correct in that the 4dvar algorithm was not converging. I must admit I am struggling to find a balance between performing a fast test run for the assimilation process and maintaining stability.

I have also recently found what really was causing the blowup…. Out of everything I didn’t think to double check the standard deviation file. I was under the impression that a dubious std file would produce dubious normalization coeffs, but they seemed fine. Anyways I used the matlab script d_std_unblanaced.m to create the initial condition standard deviation file and it simply went wrong so u and v was in the range of 100-200 m/s for the location of the blowup…. So I am also a little confused about whether to choose unbalanced or balanced standard deviations and the importance of this choice i.e. the effect this has on the 4dvar analysis. I can’t get my head around the difference and the balance operators role for the unbalanced std computation .
Anyways I ran d_std.m instead and my std file looked fine with u and v in sensible ranges.
Hopefully this will fix my problem and I’m currently running another test with 30 loops with a 30 second time step.

I did get another suggestion I have been pondering aswell. Are my gridded observations assimilating as a smooth field or as point sources? Is the observation operator smooth? I think this is referred to as a spatial regularity problem but I'm not sure.

[arango]

The error covariance balanced operator is kind of tricky. It is based on Weaver et al. (2005) method. It imposes a multivariate constraint on the background error covariance such that the unobserved variable(s) information is extracted from the observed (assimilated) data using T-S empirical relationships, hydrostatic balance, and geostrophic balance. The balanced SSH involves the solution of an elliptical equation (bi-conjugate method) or integration of the hydrostatic equation using a level of no-motion. This the tricky part when you have shallow and depth water in your application. You can get unreasonable currents. This maybe is your case. There is lots of literature about 4D-Var.

Reference: Weaver, A.T., C. Deltel, E. Machu, S. Ricci, and N. Daget, 2005: A multivariate balance operator for variational data assimilation, Q. J. R. Meteorol. Soc., 131, 3605-3625.

You also will find lots of information in our papers: Moore et al. (2011a,b,c). The ROMS 4D-Var algorithms are described in detail in these papers.

This is 4D-Var and data are assimilated when and where observed as point data or super-observations. The observation operators in ROMS sample the model solution at the location and time of the observations. The convolution of the error covariance using a diffusion operator smooths the observation influence. We have time covariance but we haven't release that yet. This actually is more complex and requires more computer memory because we need to save the state vector in time.

[lmp4]

I think I'm starting to get my head around the balance operator but I still am a little unsure about the standard deviations scripts and which to choose... Is my following logic correct?

• 1. d_std.m - Computes the 4dvar initial condition standard deviations. If I have chosen this method the balance oporator is computed within the model run itself with the cppdefs

  Code: Select all

  define BALANCE_OPERATOR
  # ifdef BALANCE_OPERATOR
  # define ZETA_ELLIPTIC
  # endif 

  2. d_std_unblanced.m - Computes the unbalanced standard deviations, computing the balance oportaor within the .m script in order to extract the unbalanced field. If I have chosen to compute the unbalanced standard deviations do I still need to compute the balance oporator within the model run? i.e. does the balance operator still needs to be set within the cppdeffs (define BALANCE_OPERATOR)?

Through my brief reading I have learnt that the balanced and unbalanced fields can represent the baroclinic(geostrophic ) and barotropic fields so wouldn't the use of an unbalanced standard deviations (d_std_unbal..m) compared to the regular computed standard deviations (d_std.m) change the 4dvar process?
In any case the fact my domain contains an island that has been smoothed might hinder the elliptic solver so I am inclined to use d_std.m or the integration method as you mentioned.

Moving onto the second point, after going through the ROMS 4dvar tutorial again I now understand how the diffusion operator spreads/smooths the observation information according the correlation length scales. Do you know of a matlab script that would help me visualize the co-variance functions and how the information is spread.

Thanks again for your excellent guidance!