Assimilation of METAR observations: GSI vs JEDI

[spanNOAA]

I am working on assimilating METAR observations 187 and 287 in FV3 JEDI EnVar and am observing some unusual results when comparing them to the results from the GSI EnVar analysis.

Below are the model level 64 and level 1 wind, temperature and specific humidity increments comparing between GSI (left) and JEDI EnVar using the same background and the same ensemble. The GSI is configured with horizontal localization of 110 km and vertical localization of 3 levels. And the JEDI is configured with horizontal localization of 110 km and vertical localization of 3 (but I'm not pretty sure the unit of the vertical localization, cause the JEDI documentation says "the vertical lengthscale is in vertical units of model's choice").

The major surface increment differences (@ level 64) between GSI and JEDI are in the observation locations. It is evident that most of the METAR observations used in JEDI are within CONUS, whereas many METAR observations over the marine areas are assimilated in GSI. Another notable difference is that, even with the same horizontal localization, the JEDI increment spreads less horizontally compared to the GSI increment. To determine the difference in length scales between the two systems, a single observation experiment will be necessary in the future.

The third issue concerns vertical localization. Below is a comparison of wind, temperature, and specific humidity increments at model level 1 (top of the model). Humidity has reasonable increments at the top of the model in both systems, which are less than 1‰ of the surface increments. However, wind and temperature increments have similar magnitudes to their surface increments in both GSI and JEDI. This issue likely requires further investigation through a single observation experiment. Additionally, more results from HofX will be posted in the near future.

If anyone has insights regarding these issues, please let me know.

[ShunLiu-NOAA]

@spanNOAA Thank you for sharing this. Could you also rerun this test with RDAS test case? Also you may consider starting a single observation test with RDAS test case for 187 and 287.

[delippi]

@spanNOAA, if you switch to the RDAS case as @ShunLiu-NOAA suggests, that is currently only working on Hera. I'll need to add support for Orion which I will do asap. I can also help you get set up with a single observation test with ob types *87. I've been wanting to add single ob generation/testing support to RDASApp as well, but still not sure about the best way to do it. I've already started working on some generic documentation in the wiki that might be a good place for you to start.

[SamuelDegelia-NOAA]

@spanNOAA Thanks for sharing these results! I agree that a single observation test with the RDAS case would be helpful. Regarding your question about vertical localization: the units for this can be set in the geometry section of the yaml you use for BUMP. See below:

geometry:
  nml_file: "./namelist.atmosphere_15km"
  streams_file: "./streams.atmosphere_15km"
  deallocate non-da fields: true
  bump vunit: "height" # vertical unit for vertical localization [modellevel, height]

I believe the default is in model levels, so if you do not set this then it should be acting as you described in your post. But you can try to change this to see.

Another point is to make sure your resolution is set high enough for running BUMP to generate the localization files for EnVar. I found in #39 that a low resolution can add artifacts to my EnVar updates in JEDI. It is recommended to use resolution: 8.0 if possible.

[spanNOAA]

@ShunLiu-NOAA @delippi Thanks for your advice. I'll discuss with Ming how to set up a single observation test.

@SamuelDegelia-NOAA Yes, I'm currently using an 8.0 resolution. For a 13 km model resolution, JEDI can't use a 110 km horizontal localization with an 8.0 resolution and the default max horizontal grid size. Using a 110 km horizontal localization with the default max horizontal grid size causes the resolution to automatically decrease to below 3, which results in a crash. Therefore, I manually increased the max horizontal grid size to 120,000 to maintain the resolution at 8 and prevent BUMP NICAS from automatically decreasing it. For the vertical localization, I didn't set 'bump vunit', so it should be in model levels. Now, I have to figure out why the model top shows such a large increment in both systems when surface observations are assimilated with a vertical localization of 3 model layers.

[SamuelDegelia-NOAA]

@spanNOAA The wiki pages here and here have some instructions on running single observation tests. The yaml files in RDASApp/rrfs-test/testinput/msonet_singleob* are used in those wikis and are useful examples.

Additionally, you can do online single observation tests (without needing a IODA file). Here is an example for a single temperature observation at 3 km:

  observations:
    observers:
    - obs space:
        name: Aircraft
        simulated variables: [airTemperature]
        observed variables: [airTemperature]
        obsdatain:
          engine:
            type: GenList
            lats: [27.44]
            lons: [278.61]
            vert coord type: height
            vert coords: [3000]
            dateTimes: [0]
            epoch: "seconds since 2022-05-26T19:00:00Z"
            obs errors: [1.0]
            obs values: [287]
        obsdataout:
            engine:
               type: H5File
               obsfile: ./aircraft_hofx_2022052619.nc4
               allow overwrite: true
      obs operator:
        name: VertInterp
        vertical coordinate: height
        observation vertical coordinate: height
        interpolation method: log-linear

[delippi]

@spanNOAA, if you switch to the RDAS case as @ShunLiu-NOAA suggests, that is currently only working on Hera. I'll need to add support for Orion which I will do asap. I can also help you get set up with a single observation test with ob types *87. I've been wanting to add single ob generation/testing support to RDASApp as well, but still not sure about the best way to do it. I've already started working on some generic documentation in the wiki that might be a good place for you to start.

@spanNOAA, I just wanted to let you know that the GSI test case support was added for Orion this morning. Running hofx validation wiki was also updated and should be able to set you set up with running the same case that we have been using (2022052619).

[SamuelDegelia-NOAA]

I ran a single observation test using the RDASApp test case for FV3-JEDI to check if we were experiencing similar problems with vertical localization that @spanNOAA presented today. I used the test yaml RDASApp/rrfs-test/testinput/msonet_singleob_airTemperature_fv3jedi.yaml that assimilates the single mesonet surface observation in RDASApp/rrfs-test/obs/msonet_singleob_airTemperature.nc4. Also, I modified the yaml to use pure EnVar instead of the hybrid (similar to the configuration described by @spanNOAA.

Increment profile at observation location:

Overall it seems that vertical localization here is working correctly. I believe the ensemble localization files from bump are using 0.3 ln(p) as the localization radius. @TingLei-NOAA Would we happen to have the yaml files used for BUMP for the fv3-jedi case in RDASApp? Sharing that might help @spanNOAA track down his problem.

[TingLei-NOAA]

@SamuelDegelia-NOAA Thanks for those clarification. I will upload the bump yaml file (used to create the localization bump/nicas files) tomorrow when Hera come back.

[spanNOAA]

@SamuelDegelia-NOAA This result is very convincing. I'll try to generate ensemble localization files separately from the analysis process and see how that works out. @TingLei-NOAA Please let me know when the bump YAML file is available in the repository. Thanks.

[TingLei-NOAA]

@spanNOAA I will double-check. But I realize that would be the same yaml I shared with you before (maybe of the name : d-bump.yaml). Anyway, I will make a PR to add it in the repo to allow all interested developers have access to it.

[ShunLiu-NOAA]

I will double-check. But I realize that would be the same yaml I shared with you before (maybe of the name : d-bump.yaml). Anyway, I will make a PR to add it in the repo to allow all interested developers have access to it. @SamuelDegelia-NOAA Thank you for the fast test.

[HuiLiu-NOAA]

@spanNOAA : Could you share the run directory of this case? I would like to take a look of the details of your settings. Thanks.

[TingLei-NOAA]

@spanNOAA I added the yaml used to create the bump loc files used in this fv3-jedi case https://github.com/NOAA-EMC/RDASApp/pull/53. It should be the same file with the name d-bump.yaml I shared with you before. Hope this example would help for you and other rdasapp developers.

[spanNOAA]

@ShunLiu-NOAA You can find it on Orion: /work/noaa/wrfruc/span/jedi_summer_retro/jedi_prepbufr_analysis/2023061012_singleobs.

[delippi]

@ShunLiu-NOAA You can find it on Orion: /work/noaa/wrfruc/span/jedi_summer_retro/jedi_prepbufr_analysis/2023061012_singleobs.

Thanks @spanNOAA, I was interested too. First, it looks like your prepbufr.surface.nc link is broken and some of your permissions are locked for users outside of wrfruc. I'm curious what kind of single ob this is? Is it a temperature ob? It isn't clear from the yaml and I can't inspect the IODA file to see. It might be good to simplify the yaml to only have the filters/functions for the single ob type, but not completely necessary--just easier to debug things. For example, I have a separate yaml for airTemperature and specificHumidity and do the testing one at a time.

I noticed a couple of things in your current yaml related to ob error assignment:

• airTemperature ob error = 2.2582 K where it looks like it should be 1.2595 K based on entry 187 second column of the GSI error table. You also use ObsFunction/ObsErrorModelStepwiseLinear for ob errors that are constant with height. You could do something like this to make it simpler:

       # ------------------
       # airTemperature
       # ------------------
       # Initial error assignment
       - filter: Perform Action
         filter variables:
         - name: airTemperature
         where:
         - variable: ObsType/airTemperature
           is_in: 187
         action:
           name: assign error
           error parameter: 1.2595  # K
         defer to post: true

• specificHumidity ob error = 0.5912 where it looks like it should be 0.027059 based on entry 187 third column of the GSI error table. I don't think you're using humidity ob, but I don't think ObsFunction/ObsErrorModelHumidity is what you want either, strangely. You just want it something like this. Make sure the units are correct for RH ob error! Look at the comment after error parameter below.

       # ------------------
       # specificHumidity
       # ------------------
       # Initial error assignment
       - filter: Perform Action
         filter variables:
         - name: specificHumidity
         where:
         - variable: ObsType/specificHumidity
           is_in: 187
         action:
           name: assign error
           error parameter: 0.027059  # GSI errtable is units of RH %/10 so 0.27059 --> 0.027059
         defer to post: true

[spanNOAA]

@delippi I’ve changed the permissions for the obs file. Please check if it’s working for you now.

Regarding the humidity section, it’s not present in my YAML file. You might have been looking at the wrong YAML file. The one I’m using is 3denvar_adpsfc_t.yaml.

[guoqing-noaa]

Thanks @delippi for suggestions!

Where did you get those obseror values? Here is what I found in fix/gsi/errtable.rrfs:

187 OBSERVATION TYPE
  0.11000E+04 0.22585E+01 0.59120E+00 0.10000E+10 0.53890E+00 0.10000E+10
  0.10500E+04 0.22585E+01 0.59120E+00 0.10000E+10 0.53890E+00 0.10000E+10
  0.10000E+04 0.22585E+01 0.59120E+00 0.10000E+10 0.53890E+00 0.10000E+10
  0.95000E+03 0.22585E+01 0.59120E+00 0.10000E+10 0.53890E+00 0.10000E+10
  0.90000E+03 0.22585E+01 0.59120E+00 0.10000E+10 0.53890E+00 0.10000E+10

The important thing is to use the same obs error for both GSI and JEDI.

[HuiLiu-NOAA]

@spanNOAA : Have you tried different values of vertical length-scale and saw any changes/sensitivity in the analysis increment at the upper levels? I suspect the vertical localization may not work for some reasons.

[guoqing-noaa]

Per discussions in this PR https://github.com/NOAA-EMC/RDASApp/pull/53, it is possible that 0.3 instead of 3.0 should be used for vertical localization. @spanNOAA Could you try a test using 0.3 instead of 3? Thanks!

[delippi]

@delippi I’ve changed the permissions for the obs file. Please check if it’s working for you now.

Regarding the humidity section, it’s not present in my YAML file. You might have been looking at the wrong YAML file. The one I’m using is 3denvar_adpsfc_t.yaml.

@spanNOAA, yes I was looking at the wrong yaml! Sorry! But I think @guoqing-noaa said it right here:

The important thing is to use the same obs error for both GSI and JEDI.

It doesn't really matter at this point of testing what the exact value is just that we are using the same between the two systems. And it looks like you are using 2.2585 in both. I think I see what the issue is though. You only specify the initial ob error in JEDI, but there is no ob error inflation. It looks like the GSI:Erriv_Input=0.4427717 and the GSI:Errinv_Final=0.8855435. So what is happening is the GSI is deflating the ob error which I think would match what you were showing in your ppt. You should be able to add ob error inflation something like this. Double check the final ob error in the JEDI hofx file and make sure this works as expected.

         # Ajusted error after initial assignment (qcmod.f90)
         - filter: Perform Action
           filter variables:
           - name: airTemperature
           where:
           - variable: ObsType/airTemperature
             is_in: 187
           action:
             name: inflate error
             inflation variable:
               name: ObsFunction/ObsErrorFactorConventional
               options:
                 test QCflag: PreQC
                 test QCthreshold: 3
                 inflate variables: [airTemperature]
                 pressure: MetaData/pressure
                 distance threshold: -1.
           defer to post: true

         # error inflation based on pressure check (setupt.f90)
         - filter: Perform Action
           filter variables:
           - name: airTemperature
           where:
           - variable: ObsType/airTemperature
             is_in: 187
           action:
             name: inflate error
             #inflation factor: 1.09757
             inflation variable:
               name: ObsFunction/ObsErrorFactorPressureCheck
               options:
                 variable: airTemperature
                 inflation factor: 8.0
                 geovar_sfc_geomz: surface_geometric_height
           defer to post: true

I think your ombg/hofxs also look a little "too" different. I still suggest using the same case as we are using at EMC. I have some pretty useful plotting tools for easily plotting and checking these sorts of things. I'm more than happy to help you get set up with that.

[spanNOAA]

inflation

@delippi That's a great point. I'll double-check the final observation error and see how things work when both systems have the same final error.

[guoqing-noaa]

@delippi Thanks for identifying the final obs error used.

@spanNOAA Please set l_gsd_terrain_match_surfTobsto falsein gsiparm.anl. The obs error and observation will be changed when matching the surface T observation elevation with the model terrain.

[spanNOAA]

@delippi I’ve changed the permissions for the obs file. Please check if it’s working for you now.

Regarding the humidity section, it’s not present in my YAML file. You might have been looking at the wrong YAML file. The one I’m using is 3denvar_adpsfc_t.yaml.

@spanNOAA, yes I was looking at the wrong yaml! Sorry! But I think @guoqing-noaa said it right here:

The important thing is to use the same obs error for both GSI and JEDI.

It doesn't really matter at this point of testing what the exact value is just that we are using the same between the two systems. And it looks like you are using 2.2585 in both. I think I see what the issue is though. You only specify the initial ob error in JEDI, but there is no ob error inflation. It looks like the GSI:Erriv_Input=0.4427717 and the GSI:Errinv_Final=0.8855435. So what is happening is the GSI is deflating the ob error which I think would match what you were showing in your ppt. You should be able to add ob error inflation something like this. Double check the final ob error in the JEDI hofx file and make sure this works as expected.

         # Ajusted error after initial assignment (qcmod.f90)
         - filter: Perform Action
           filter variables:
           - name: airTemperature
           where:
           - variable: ObsType/airTemperature
             is_in: 187
           action:
             name: inflate error
             inflation variable:
               name: ObsFunction/ObsErrorFactorConventional
               options:
                 test QCflag: PreQC
                 test QCthreshold: 3
                 inflate variables: [airTemperature]
                 pressure: MetaData/pressure
                 distance threshold: -1.
           defer to post: true

         # error inflation based on pressure check (setupt.f90)
         - filter: Perform Action
           filter variables:
           - name: airTemperature
           where:
           - variable: ObsType/airTemperature
             is_in: 187
           action:
             name: inflate error
             #inflation factor: 1.09757
             inflation variable:
               name: ObsFunction/ObsErrorFactorPressureCheck
               options:
                 variable: airTemperature
                 inflation factor: 8.0
                 geovar_sfc_geomz: surface_geometric_height
           defer to post: true

I think your ombg/hofxs also look a little "too" different. I still suggest using the same case as we are using at EMC. I have some pretty useful plotting tools for easily plotting and checking these sorts of things. I'm more than happy to help you get set up with that.

With these filters applied, I got errors indicating that VADER cannot generate geopotential_height and surface_geometric_height. Do you have any insight into what might be causing this issue?

[delippi]

@spanNOAA, I'm not too sure without digging into your particular case more--it could be something related to how your particular case was create--I can't say for sure.

I really would like to urge you to try this wiki which will set up the case we are using at EMC. I'm almost certain that once you regenerate your single ob you can just change the bufr type from 188 to 187 and the initial ob error in the yaml that I provided called testinput/msonet_airTemperature.yaml and you shouldn't be having any of these problems. It'll also make digging into these sorts of issue a little easier.

Basically all you have to do is:

cd RDASApp/rrfs-test/scripts
bash setup_experiment.sh # Make sure to edit USER INPUT block first

And you should be good to go with running the test case.

[spanNOAA]

Based on the discussion in #67 , three more single observation experiments were conducted and are shown below. GSI h110km: Number of values falling between 0.050857356826984675 (maxval/e) and 0.13824462890625 (maxval): 12 JEDI h110km: Number of values falling between 0.05738280633564391 and 0.1559826397281654: 6 JEDI h311km: Number of values falling between 0.05985338604609267 and 0.1626983716608379: 9 JEDI h401km: Number of values falling between 0.05958699176913942 and 0.16197423693859037: 11

[guoqing-noaa]

@spanNOAA Thanks a lot for those experiments and plots. I will copy your photos and results to issue #67 and continue some discussions there.