Add density-dependent interpolation of BGC tracers to ROMS-Tools BC and IC

[kthyng]

Used UCLA matlab script to start (Add_bgc2ini.m). That code uses nearest neighbor to interpolate in the water column within grid cells, and uses the same for extrapolation. This code uses linear interpolation to interpolate and nearest neighbor to extrapolate.

The interpolation methods give similar results in some areas and different results in others. Here is an example boundary from the end to end notebook domain (Gulf of Mexico eastern boundary), depth-based interpolation:

and density-based interpolation — note the sudden change in behavior midway through the transect:

Example profile from eta_rho=10 where they look similar in the boundary slice images. The reasonable density profile leads to a reasonable resulting alkalinity.

Example profile from eta_rho=29 where they do not look similar in the boundary slice images. The uniform density profile leads to a strange resulting alkalinity.

https://cworthy.atlassian.net/browse/CSD-869

• Tests added
• Passes pre-commit run --all-files
• Changes are documented in docs/releases.md
• New functionality has documentation

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[kthyng]

I set this to be the new default — not sure if that is the desired behavior. I had to update some of the test data accordingly.

[kthyng]

@blsaenz I know you're out this week but when you get a chance can you let me know if you think someone else should instead or also be a reviewer?

[blsaenz]

Looking now!

[blsaenz]

I am curious about these array transpose calls. Looking at the dask dashboard I saw the transpose is 90% of the processing for the roms-tools tidal preparation. I wonder if we could avoid this somehow. I think if it works I'm happy to approve as is, but it would be good to know how we got into the situation of needing to transpose dimensions.

[kthyng]

From using xgcm in the past (xgcm does the vertical interpolation here), I remember always having to transpose afterward back into time x depth x lat/lon proper order, but I can experiment a bit and see if I can avoid it.

[kthyng]

I have just been running this and looking at the dask graph and profile and I am not seeing transpose being a lot of the processing for the example I am running. For the example you are trying it with, do you see a big difference for the density-based interpolation compared to the depth-based interpolation?

As a side note, as I am trying another thing, I do see transpose being a lot of the initial call for the BGC boundary forcing with 2D fill on (it is not present with 1D fill), but I don't think the actual interpolation happens until the save function is called subsequently.

[kthyng]

And I don't think I can avoid needing the transpose after using xgcm.

[blsaenz]

I am a bit confused by the plots. The density-based one seem invalid for a slice with the sharp changes? Does that mean we have extreme density differences in adjacent cells also? I can see some differences near the surface which is where I expected the major changes, but it seems like something is off, maybe in the lateral fill? Can we try with 2D interpolation on for fill?

[kthyng]

I am a bit confused by the plots. The density-based one seem invalid for a slice with the sharp changes? Does that mean we have extreme density differences in adjacent cells also? I can see some differences near the surface which is where I expected the major changes, but it seems like something is off, maybe in the lateral fill? Can we try with 2D interpolation on for fill?

I can try to 2D fill and compare to see if that changes anything, but I think the interpolation is valid — it just gives weird results when the density is mostly uniform. Here are the density vs Alk at the same weird location:

So the alkalinity is being interpolated from the blue x densities to the orange circle densities, but since the orange circles have almost entirely the same density, they get the same alkalinity which is calculated using the x values in the plot. Does that make sense? There may be a better plot to show to illustrate this.

[blsaenz]

Ok, thanks for looking into it. A few of our roms-tools processing steps are pretty opaque, so I’m very glad that you are able to examine in more detail than we have before! On May 20, 2026, at 2:37 PM, Kristen Thyng ***@***.***> wrote: @kthyng commented on this pull request. In roms_tools/setup/utils.py:

+ n_depth = density.sizes[phys_depth_dim]

+ perturbation = xr.DataArray(np.arange(n_depth) * 1e-7, dims=[phys_depth_dim]) + density = density + perturbation + + # xgcm.transform requires a single chunk along the dim being transformed. + density = density.chunk({phys_depth_dim: -1}) + + # Regrid density from physics depth levels to BGC depth levels + ds_phys = xr.Dataset({phys_depth_dim: phys_depth_coord}) + vertical_regrid = VerticalRegrid(ds_phys, source_dim=phys_depth_dim) + source_density = vertical_regrid.apply( + density, + source_depth_coords=phys_depth_coord, + target_depth_coords=bgc_depth_coord, + ) + source_density = transpose_dimensions(source_density) I have just been running this and looking at the dask graph and profile and I am not seeing transpose being a lot of the processing for the example I am running. For the example you are trying it with, do you see a big difference for the density-based interpolation compared to the depth-based interpolation? As a side note, as I am trying another thing, I do see transpose being a lot of the initial call for the BGC boundary forcing with 2D fill on (it is not present with 1D fill), but I don't think the actual interpolation happens until the save function is called subsequently. —Reply to this email directly, view it on GitHub, or unsubscribe.Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you were mentioned.Message ID: ***@***.***>

[kthyng]

Ok actually with 2D fill this boundary plot does look different right in the area that I was showing my test case from (eta_rho=29):

So I will investigate further on how this relates.

[blsaenz]

Now that I see the results I can't believe the UCLA's version could give good results either, interpolating alkalinity or other BGC values directly against density. Is UCLA somehow normalizing the density (source and target), so that the whole profile is always mapped? Perhaps we should try something like that?