bug: pert_wrf_bc code introduces artifacts to perturbed WRF boundary files

[braczka]

Describe the bug

Potential bug reported by Rui Cheng (UCSD). See dart(at)ucar.edu for more details.

1. Recreate bug: Using the pert_wrf_bc code when running a WRF-DART simulation. The pert_wrf_bc code is invoked both when running init_ensemble_var.csh and driver.csh within the WRF shell scripting.

2. What was the expected outcome? The perturbation of the wrfbdy file should reflect the WRFDA CV3 perturbation which is applied to the wrfinput file, then conveyed to the wrfbdy files. Perturbations are applied to atmospheric variables (e.g. wind, geopotential, temp, precip fields) to enhance ensemble spread during DA.

3. What actually happened? ### Error Message Please provide any error messages.

Artifacts are introduced within the PH (perturbed geopotential) field, that are unrelated to the expected perturbations introduced through perturbation bank. There is no specific error message -- because there is no quality checks performed on the perturbed fields -- however, this could lead to instabilities during subsequent WRF integration (e.g. CFL errors etc).

These artifacts seems to coincide with mountainous terrain. So if the boundary conditions do not coincide with mountainous terrain, this may not be noticeable.

Screenshots

Below I show a screenshot of the wrfbdy fields for U, V, W, PH, T, and QVAPOR for each of the wrf domain boundaries within the WRF-DART tutorial. The wrfbdy 'mean' files generated from the gen_retro_icbc.csh step should be similar (nearly identical) to the mean of the perturbed ensemble of wrfbdy files generated by pert_wrf_bc. If there are signficant pattern differences -- this indicates a bug. Note the artifact for the PH variable (perturbed geopotential) for the southern boundary (BYS_pertmean). Red arrow. The fact that there is a difference between BYS_mean and BYS_pertmean indicates a problem.

Definitions:

BYS_mean: southern boundary generated from gen_retro_icbc_csh BYS_pertmean: mean of perturbed southern boundary ensemble members (pert_wrf_bc)

BYE_mean: northern boundary generated from gen_retro_icbc_csh BYE_pertmean: mean of perturbed northern boundary ensemble members (pert_wrf_bc)

BXS_mean: western boundary generated from gen_retro_icbc_csh BXS_pertmean: mean of perturbed western boundary ensemble members (pert_wrf_bc)

BXE_mean: eastern boundary generated from gen_retro_icbc_csh BXE_pertmean: mean of perturbed eastern boundary ensemble members (pert_wrf_bc)

Version of DART

Latest release.

Have you modified the DART code?

No, but a proposed fix includes replacing the existing lines within pert_wrf_bc with something more related to current WRFv4+ code..

This fix is suggested to work by Rui. I am currently in contact with MMM to verify this fix as well.

Build information

Derecho build with gfortran

[braczka]

Apologies -- Rui Sun (UCSD) first reported this issue.

[braczka]

Update -- I tested the same WRF-DART setup as described above except reverted WRF 4.5.2 to terrain following coordinates by setting &dynamics namelist option to hybrd_opt = 0. The expectation was that this would fix the the perturbation geopotential artifact, however, this did not fix the issue. I need to follow up on why this is the case, but I need to take a short pause on this issue for now.

Rui Sun is also checking this with his setup, and I also have a question into Arthur about how WRF-Chem is influenced by the perturbation issue.

[braczka]

Update -- I tested the same WRF-DART setup as described above except reverted WRF 4.5.2 to terrain following coordinates by setting &dynamics namelist option to hybrd_opt = 0. The expectation was that this would fix the the perturbation geopotential artifact, however, this did not fix the issue. I need to follow up on why this is the case, but I need to take a short pause on this issue for now.

Rui Sun is also checking this with his setup, and I also have a question into Arthur about how WRF-Chem is influenced by the perturbation issue.

Rui Sun followed up on this issue and found that reverting back to the terrain following vertical coordinate system using WRFv4+ resolves this boundary condition perturbation issue for their application. Performing the WRF-DART assimilation with terrain following (hybrid_opt= 0) coordinates with current DART boundary perturbation methodology gives same (scientifically indistinguishable) results as compared to using the default sigma hybrid vertical coordinate (hybrid_opt=2) and local DART edits to BC perturbation code. I am currently following up on this to confirm.

[braczka]

I ran two side by side nested WRF-DART example cases based on support request from Kecheng (Lili Lei group) for a 3 member ensemble -- the first example used sigma hybrid coordinates to generate the initial conditions (real.exe), I then generated and applied perturbations using WRFDA, and then ran a 6 hour forecast. The second case was identical except it used terrain following coordinates (hybrid_opt=0) for all steps. The case where I used the sigma hybrid coordinates developed a strong P (and PH) gradient at the boundary of the parent domain (d01) at the 6 hour forecast, whereas the terrain following coordinates looked normal. See figures below for illustration of spurious boundary gradient (3 member ensemble of d01 at level 1 and 16). Proceeding with side by side nested assimilation example to quantify the impact of this incorrect perturbation setup (strong boundary gradient) impacts skill of assimilation and forecast.

[braczka]

I followed through with a side by side comparison (sigma hybrid (SH) vs. terrain following (TF) coordinate systems) applied to a real case (May 19 2024 Derecho over Kansas) to quantify the impact that these incorrect boundary gradients have on forecast skill. I ran a 20 member ensemble nested WRF-DART case using the WRF shell scripts, generalized for a nested domain setup. The initial perturbations were applied at 2024051900 and the simulations ran for 24 hours total (2024052000). I allowed for 6 hours of spinup followed by the assimilation of NCEP prepbufr observations (e.g radiosonde, ACARS etc) at hour 6 and again at hour 12. Forecasts were then run through hours 18 and 24. Observation space diagnostics were performed at hour 12 and at hour 18.

Key findings:

1. The TF simulation performed significantly better than the SH for all observation space diagnostics in terms of RMSE for temperature, winds and specific humidity at all time periods that I looked at (@ assimilation hour 12, and @ forecast hour 18). See results below.
2. The TF simulation compared more favorably than the SH in terms of forecasting the timing and spatial structure of the derecho storm complex as judged from local radar reflectivity. See results below.

Conclusions:

1. The results suggest that the unintended perturbations from the boundary gradients when using SH coordinates with pert_wrf_bc significantly degrades the WRF forecast skill.
2. The performance of the SH coordinate system, although poor compared to the TF, still exceeds that of a free WRF run. Given that the SH simulations completed successfully and show improvement compared to 'free' forecasts, it is likely only advanced WRF users would detect any reduced performance.
3. WRF-DART code should restrict use to only TF coordinate setup. Should not allow for use with SH coordinates given relative poor performance.

Caveats/disclaimer:

1. I specifically chose a domain setup that coincided with mountainous terrain (on western boundary) because this tends to amplify boundary gradients leading to unintended perturbations. The WRF user guide specifically recommends against generating domain boundaries that coincide with mountainous terrain. Using boundaries that did not coincide with mountainous terrain would likely improve the performance of the SH relative to the TF coordinate system -- perhaps even make forecast skill indistinguishable -- I did not test for this.
2. My WRF model setup was not optimized in terms of dynamics and microphysics schemes nor in terms of nested domain positioning and resolution ratio relative to parent domains (2:1). However -- the key point is that the WRF model setup was identical in every way except for the coordinate system, thus was a valid test for the impact of the boundary gradient on forecast skill.

Observation space diagnostics:

Storm Cell Position/Timing: