How to do regional surface salinity restoring (nudging) in POP/CESM2.1.5? (target climatology + mask file format)

[CCl]

What version of the code are you using?
CESM version: CESM2.1.5


Have you made any changes to files in the source tree?
No source-code changes (only case-level settings / user_nl_pop / input files).


Describe every step you took leading up to the problem:
1. Create a fully coupled PI control case (any standard PI compset is fine for this question), then run long enough to create a climatology for surface salinity (e.g., 100-yr monthly climatology).

2. Build a target restoring file from POP history output:
- Start from a climatology of salinity (3D).
- Apply a regional perturbation only within a limited Southern Hemisphere mid/high-latitude band (example: a latitude band like ~<lat1>–<lat2>), e.g. target = climatology + 0.5 (psu or g/kg; only surface) inside the region, unchanged elsewhere.

3. Prepare a regional mask file so that restoring is applied only in the selected region (land points excluded, ocean points outside region excluded).

4. Enable salinity restoring/nudging in POP via user_nl_pop (I think it is handled by forcing_s_interior / pop_in namelist).
- I am not fully sure about the correct namelist variables and the expected netCDF variable names / dimensions (details below).

-----------------------my user_nl_pop ---------------------------------
&forcing_s_interior_nml
s_interior_formulation = 'restoring'
s_interior_data_type = 'annual'
s_interior_file_fmt = 'nc'
s_interior_filename = '/path/to/your/annual_SALINITY_3D_clim.nc'
s_interior_surface_restore = .true.
s_interior_restore_max_level = 1
s_interior_variable_restore = .true.
s_interior_restore_file_fmt = 'nc'
s_interior_restore_filename = '/path/to/your/S_RESTORE_MAX_LEVEL_and_S_RESTORE_RTAU_2D.nc'
/
----------------------------------------------------------------------------

I am attempting to apply regional restoring of surface salinity in POP (CESM2.1.5) using forcing_s_interior_nml, with restoring restricted to the top model level only (s_interior_surface_restore=.true. and s_interior_restore_max_level=1) and with a user-generated target climatology.

My main issue is that the model crashes shortly after enabling salinity restoring. The run aborts with a MARBL/CO2SYS pH solver error shortly after enabling restoring; I suspect salinity restoring may be triggering unphysical states, but I am not yet sure. I am not sure whether this indicates:

1. a misunderstanding on my side about how forcing_s_interior restoring should be configured for surface-only and regional-only salinity restoring, or
2. POP/MARBL (and possibly sea ice) being very sensitive to salinity restoring if the input target/mask has any issues (units, missing values, land points, sharp boundaries, etc.).

Additional context:

1. I also tested a configuration with only ocean + sea ice coupling (to simplify debugging), but the run still fails after enabling salinity restoring, and in some attempts the sea ice component (CICE) aborts soon after the run begins.
2. Without salinity restoring (same case otherwise), the model runs normally.

What I would like to confirm / ask:
1. Is forcing_s_interior_nml the recommended/supported pathway to do surface-only (k=1) and regional-only salinity restoring in POP?
2. Are there known best practices / pitfalls to avoid pH/chemistry failures when restoring salinity (e.g., required handling of land points, fill values, units, tapering at the boundary, recommended restoring strength)?
3. What are the key sanity checks you recommend to verify that the restoring target and mask are not introducing unphysical salinity tendencies (e.g., min/max checks, ensuring restoring is strictly zero outside the region/over land, etc.)?


The error I had met:
----
ERROR: bounding bracket for pH solution not found
co2calc.F90: drtsafe_row
----

 

[mlevy]

I suspect salinity restoring may be triggering unphysical states,

Can you provide more of your log output? Specifically, I'm looking for something like the output in this comment; MARBL should report on the values of DIC (dic), alkalinity (ta), phosphorus (pt), silicon (sit), temperature (temp), and salinity (salt), and it will be clear if something is very wrong.

Another possibility is that physical values are reasonable, but your restoring is somehow leading to more ice melt and some of the poleward cells that are typically ice-covered are now open ocean -- this can lead to CFL violations that are resolved by taking a smaller time step. I think this is unlikely though, because we usually see those issues in the northern hemisphere (we hide the north pole in Greenland and as the coastal sea-ice melts it exposes very small grid cells; in the southern hemisphere Antarctica is large enough that the small grid cells are all land rather than hiding under sea ice)

 

[CCl]

Hi Michael, thanks for taking a look.
Here is the MARBL output around the crash (this is right after enabling salinity restoring):

(Task 60, block 1) MARBL ERROR (marbl_co2calc_mod:drtsafe): bounding bracket for pH solution not found
(Task 60, block 1) MARBL ERROR (marbl_co2calc_mod:drtsafe): ... dic = 0.1732991E+005
(Task 60, block 1) MARBL ERROR (marbl_co2calc_mod:drtsafe): ... ta = 0.1863307E+005
(Task 60, block 1) MARBL ERROR (marbl_co2calc_mod:drtsafe): ... temp = -0.1969854E+001
(Task 60, block 1) MARBL ERROR (marbl_co2calc_mod:drtsafe): ... salt = 0.2343485E+003

So it does look like something is going very wrong (salt ~234, and DIC/TA are also extremely large), which makes me suspect the restoring is driving the model into an unphysical state.

What’s confusing is that I checked my restoring target file beforehand and the salinity values look reasonable (max is only ~40 g/kg). Do you have any idea why MARBL/POP would end up seeing such a large salinity if the input target field itself is not that high? Could this be due to a units / scaling issue, or the way the restoring term is applied?

Related question: for the restoring timescale *_restore_tau, what units does POP expect? I set it to 365 assuming “365 days”, but is it actually interpreted as seconds? If so, that would make the restoring far too strong, which might explain the blow-up — but it seems a bit surprising, so I wanted to confirm.

Thanks again for your time and guidance.

 

[CCl]

One more detail: the failing location appears to be in the Arctic, around (lon,lat) ≈ (316–318°, 84–85°N), roughly i,j ≈ (162–167, 375–378). If that’s a region where small/open-ocean cells can be exposed during sea-ice retreat, I wonder if that might also be relevant.

 

[mlevy]

Do you have any idea why MARBL/POP would end up seeing such a large salinity if the input target field itself is not that high? Could this be due to a units / scaling issue, or the way the restoring term is applied?

Unfortunately the inline documentation for the restoring module is pretty sparse, so it's unclear what units you should be using. I believe internally POP uses msu (or g/g) for salinity rather than psu (or g/kg), so it's expecting values around 0.035 rather than values around 35. You could try setting s_interior_data_renorm=0.001, which would would scale the values in your forcing file by 1/1000. (Sorry, this is all stuff I should know rather than guess at, but it's been years since I've had to dig into the POP code.)

Related question: for the restoring timescale *_restore_tau, what units does POP expect? I set it to 365 assuming “365 days”, but is it actually interpreted as seconds? If so, that would make the restoring far too strong, which might explain the blow-up — but it seems a bit surprising, so I wanted to confirm.

I see a line s_interior_restore_rtau = c1/(seconds_in_day*s_interior_restore_tau), so you are correct that s_interior_restore_tau has units days and then it is converted to seconds in the body of the code.

One more detail: the failing location appears to be in the Arctic, around (lon,lat) ≈ (316–318°, 84–85°N), roughly i,j ≈ (162–167, 375–378). If that’s a region where small/open-ocean cells can be exposed during sea-ice retreat, I wonder if that might also be relevant.

This definitely looks like the CFL issue, but I bet it's related to the units issue in the salinity restoring; I suspect setting s_interior_data_renorm will fix it. If the model still crashes (but with reasonable salinity and DIC values) we can talk about reducing the timestep.

 

[CCl]

Unfortunately the inline documentation for the restoring module is pretty sparse, so it's unclear what units you should be using. I believe internally POP uses msu (or g/g) for salinity rather than psu (or g/kg), so it's expecting values around 0.035 rather than values around 35. You could try setting s_interior_data_renorm=0.001, which would would scale the values in your forcing file by 1/1000. (Sorry, this is all stuff I should know rather than guess at, but it's been years since I've had to dig into the POP code.)



I see a line s_interior_restore_rtau = c1/(seconds_in_day*s_interior_restore_tau), so you are correct that s_interior_restore_tau has units days and then it is converted to seconds in the body of the code.


This definitely looks like the CFL issue, but I bet it's related to the units issue in the salinity restoring; I suspect setting s_interior_data_renorm will fix it. If the model still crashes (but with reasonable salinity and DIC values) we can talk about reducing the timestep.

Hi Michael — thank you so much for the suggestion! Setting s_interior_data_renorm=0.001 did the trick: the salinity values are now in a reasonable range and the run no longer blows up. This completely solved the issue on my end and saved me a lot of time. Really appreciate you taking the time to look into it and share your insight !