Assimilation of surface observations in LETKF: GSI vs. JEDI

[SamuelDegelia-NOAA]

Extensive work is ongoing to validate H(x) and observation errors in FV3-JEDI compared to GSI (e.g., #80). However, we also need to evaluate the solver component for both EnVar and EnKF. EnVar comparisons are contained in #46 whereas this issue is created to compare LETKF between FV3-JEDI and GSI. I will follow up with additional updates as I perform these comparisons.

Some notes: the default solver in GSI is the EnSRF but there is also an option for LETKF (LETKF_FLAG = T) so I will use it for comparisons. In FV3-JEDI, there is also an option for solver: GSI LETKF to emulate the single precision fortran implementation using LAPACK. I also plan to use @delippi's GSI branch to output GSI geovals and input those into FV3-JEDI. That way we can ensure that H(x) is equivalent between GSI and FV3-JEDI and so we can isolate comparisons of the solver only.

[delippi]

Hi @SamuelDegelia-NOAA, I'm glad to see this. I need to work on putting my changes to ufo's vertInterp operator into a branch. In the mean time, you could use my RDASApp build that I made for the demo. That build already has those changes to directly apply GSI geovals to vertInterp. You can also use the GSI build that I have there that outputs the geovals that you need.

You could also make this comparison following almost exactly what we did in the demo yesterday. Start by replicating an exact one-to-one comparison of the hofx's and oberrors (you said you would turn off oberror inflation in GSI and then just remove that part of the JEDI yaml to inflate errors). I assume you are going to, but I just wanted to say to use the mesonet obs since I've already validated those yamls. And it should also go without saying that a good first step is always a single ob test!

I think all of that should be pretty straight forward to do. I look forward to seeing the results!

My demo directory on Hera is here: /scratch2/NCEPDEV/fv3-cam/Donald.E.Lippi/RRFSv2_demo.

[SamuelDegelia-NOAA]

@delippi Thank you for the feedback! Sounds like we are on the same page here. My plan is to use your RDASApp and GSI builds to ensure that H(x) is an exact fit between GSI and JEDI. I will just copy these builds from your installation on Hera since that will be easiest. I am also starting with a single observation test from the mesonet obs. For the observation error, I decided to just write a quick python script that replaces the output observation error in the diag file with the input error so that I do not have to track down every inflation routine in GSI. Before I do any comparison, I'll make sure that the obs errors match in addition to H(x).

I am still in the early staged of this validation (just got the GSI-LETKF up and running for our test case) but will hopefully have some results soon.

[SamuelDegelia-NOAA]

As an update, I have been unable to get the GSI-based LETKF working for a single observation test. Assimilating the full observation suite works, but a single observation produces very small increments (~10^-5). The EnSRF does not have this problem. I also tried using some older test case data I had and it has the same problem. It seems that the GSI implementation of LETKF has some issues with a single observation test. So I am now proceeding with this validation using EnSRF which will introduce some (likely small) differences from the JEDI-based LETKF.

Below is an initial result comparing the assimilation of a single mesonet temperature observation between GSI and JEDI. The o-b, observation error, and background spread are the same. However there are some differences in the extent and magnitude of the analysis increments. The increments in FV3-JEDI are slightly larger and more widespread compared to GSI (despite using the same localization cutoff). I will provide additional updates as I continue to investigate these differences.

[SamuelDegelia-NOAA]

Increment profiles at the observation location reveal some larger differences around 700 hPa. The rest of the increment profile looks pretty good. Note that these experiments do not use any vertical localization since this functionality does not exist in OOPS (except for GETKF).

[delippi]

Hi @SamuelDegelia-NOAA, this is some good investigative work. I'm thinking that for single ob test that EnSRF and LETKF should be formally equivalent and provide identical results, do you agree?

[SamuelDegelia-NOAA]

Theoretically they should be very close for a single observation test. But I do not think they will be exactly the same due to differences in localization methods. The EnSRF applies localization to B while the LETKF applies localization via R. To my understanding, the latter method inflates obs errors based on the distance between the current model grid point and obs within the radius. But I should confirm that this is how things work in JEDI/OOPS.

This does make it hard to confirm whether the differences above are due to localization differences or something else. But I would lean to there being another difference that I am missing. At the very least, I would think that the increment and o-a would be the equivalent between LETKF and EnSRF. Also, since vertical localization is turned off, I would expect the profiles in my most recent plot to match up more closely than they do.

[SamuelDegelia-NOAA]

I reran the single observation experiments with very large localization radius (10,000 km) and the increments are now much closer. See below. The maximum increment and the increment in obs-space are also much closer now. This implies that the differences in my previous plots were largely due to different localization methods between EnSRF and LETKF. However, there are still some larger differences at 700 hPa that I cannot explain yet.

I should mention @delippi that I was not able to get the GSI geovals working for the ensemble members. The code added to UFO looks like it only reads the geovals in from the IODA file which currently only has one set of geovals in it (for the ensemble mean in my case). So the H(x) for the ensemble mean matches but the members are very slightly different. They H(x) for the members only differ by about 0.01 K but that might be enough to change the covariance structure and cause some slight differences in the increment structure like in the profile plot below. My next step is to compare the covariances computed using GSI and JEDI geovals to see if there are some large differences at 700 hPa.

[SamuelDegelia-NOAA]

For the sake of Github images being hard to view and annotate, I have moved most of this analysis to a Google Doc found here.