Changes in Elevation Fields Between SCHISM Version Tags

[jduckerOWP]

Issue

The schism-dev community has done a fantastic job overall in keeping a record of change logs on the GitHub repository in order for users to understand the development stages of the model itself. The NOAA Office of Water Prediction (OWP) currently utilizes SCHISM v5.10.0 for current operations to calculate the Total Water Level (TWL) fields across US coastlines. The model is compiled using the default settings highlighted under the schism/src/cmake/SCHISM.local.build cmake file to calculate the TWL fields using the baseline hydrodynamic core of the SCHISM code base. As OWP transitions to the NextGen coupled framework, we intend on updating the latest tag version release of the SCHISM code base in order to use the recently developed BMI code. Some work has been recently done on our end in evaluating SCHISM version tags and their sensitivity to changes in elevation fields within the out2d*.nc output files. Below is our full analysis of all SCHISM version tags since the release of v5.10.0 for current OWP operations:

v5.10.0, v5.10.1 -> same elevation fields

v5.11.0, v5.11.1, v5.11.2, v5.13.0, latest master branch -> same elevation fields

Divergence of elevation fields occur between v5.10.1 -> v5.11.0

Steps for Reproducibility

1. Git checkout v5.10.1 and v5.11.0
2. Compile both SCHISM code bases with the default src/cmake/SCHISM.local.build CMake options.
3. Ensure param.nml files mimic each other in overall numerical values (differences between these versions of param.nml files do not impact the SCHISM model settings for this test).
4. Run both SCHISM model versions for some baseline simulation containing meteorological forcings (sflux*.nc file), offshore water level boundaries (elev2d.th.nc file), and a source/sink forcing file (source.nc).
5. diff out2d*.nc files to evaluate a change in their output fields, including elevation fields. The differences grow over the time of the simulation.

Lingering Questions/Concerns
The change logs on the GitHub repository do not appear to reflect some fundamental change to the core SCHISM hydrodyanmic code base to cause such a change in calculations of elevation and velocity fields. Can someone help address what actual code change was implemented to cause a divergence in the baseline out2d*.nc fields overall? Was there some documentation in general that the code changes between these SCHISM versions resulted in an improvement in the calculations of the elevation fields? Ideally, we would just like to obtain some documentation, if possible, on the code changes and potential improvements to the SCHISM hydrodynamics for simulating the elevation fields.

I know this may a heavy lift potentially for looking through the changes in the code base. I've done some digging myself and I've identified a series of code changes that may potentially explain the divergence in the elevation field calculations that I can post in this thread if you would like to continue discussions on this topic here further. I just want to understand from the perspective of the SCHISM developers if they can easily identify for me here a fundamental code change in the SCHISM code base between these versions that would have resulted in changes to output fields within the out2d*.nc files.

Thank you all for your time and support on addressing this SCHISM issue here and I look forward to hearing your response!

[josephzhang8]

Thank you, Jason for the challenge. All of important algorithmic changes are documented inside src/Readme.beta_notes (which are then partially transferred to github change logs). I'm going to take out that change log and ask users to refer to the src file directly.

The most likely cause for the elev difference would be:

(137) cd0b2a3 (14 Aug 2022): [via Qian_pole_treatment] changed nodalvel() under ics=2; commented out a frame transformation
for better vel near poles.

This change was made to improve velocity field near poles, and as a result we had to make 'small' tweaks on how the local frames are used how under lon/lat. Still, the differences at most nodes should be negligible, and my guess is that the differences occurred at some high-elevation watershed nodes. Note that if the node is dry, elevation values do not mean anything, so those differences might not matter after all. Can you please view elev at all nodes?

[jduckerOWP]

@josephzhang8 thank you very much for your insight on the code changes that could have been attributed to perturbations between elevation fields when comparing v5.10.1 and v5.11.0. I've gone ahead and plotted up a few figures here highlighting differences in elevation fields between model output of these two SCHISM versions to just illustrate where those differences reside on the mesh. I'll give you an overview on the work I've done to just shed some light on the analysis I completed here overall.

SCHISM Mesh Extent

• This figure above is to just give you a general visual of the entire extent of the mesh developed for NOAA OWP operations along the US/Gulf coastline.

Filtering Mechanism

• All elevation fields in this analysis have filtered out the dry nodes based on the output of the out2d*.nc files

Hour 1 of SCHISM Simulation Differences

• About 18% of the nodes within the entire mesh had absolute differences in elevation fields greater than 0.0 between the two SCHISM versions of interest here.
• The average absolute difference value was only 0.001 meters, with a maximum absolute difference being capped at 0.005 meters. Overall, this is negligible in the grand scheme of magnitudes for TWL fields. The nodes that showed absolute differences greater than 0.0 though were everywhere throughout the mesh (not along boundaries) just to mention here in general.

Hour 240 of SCHISM Simulation Differences

• Now about 33% of the nodes within the entire mesh had absolute differences in elevation fields greater than 0.0 between the two SCHISM versions of interest here.
• The average absolute difference values now grew to about 0.005 meters, with maximum absolute differences now expanding to 12 meters.
• About 5% of the mesh showed absolute differences of 0.1 meters in elevation fields which extends into the scale of some concern with elevation fields diverging from these two SCHISM source code versions. The plot below shows the locations on the mesh where we've observed between 0.1-to-12-meter absolute differences in elevation fields between the two SCHISM versions:

• The nodes in the image above reside along the US coastline, where TWL fields become important for operational forecasting of surge components in certain scenarios. I guess the overall concern here from this analysis is that do we expect the code update from (137) cd0b2a3 (14 Aug 2022): [via Qian_pole_treatment] to really impact elevation fields this substantially across a given SCHISM mesh? It may indeed be the case from your previous discussion about the implications of this particular work, but I just want to be sure that the behavior here in elevation estimates between the SCHISM versions appears consistent with the frame transformation for better velocity fields near poles.

Let me know if there is any further analysis you would like me to take a deeper dive into here and I'll be happy to help!

[josephzhang8]

Thx @jduckerOWP. The nodes of concerns appear to be all located above shoreline?
If so, this may be related to the change in 1 parameter (meth_sink), which was removed and hardwired in _step as shown below:

! if(meth_sink/=0) then
where(idry_e==1.and.vsource<0.d0) vsource=0.d0
! endif !meth_sink

[josephzhang8]

And in general I don't see any red flags as long as the new results compare well with obs.
The use of slightly different algorithm/local frame can indeed lead to those changes.