Fix barotropic streamfunction calculation

[xylar]

Instead of computing the BSF using a least-squares solve based on the transport between vertices, the new approach inverts a Poisson equation for the BSF based on the (vertically integrated) vorticity.

I have proven that this produces more accurate results. The vertically integrated vorticity computed from the streamfunciton is within machine precision of the original vertically integrated vorticity. The velocity produced by the streamfunciton has errors that are on the order of 1e-3 compared with the original vertically integrated velocity (likely caused by the field not being perfectly divergence free).

Checklist

• User's Guide has been updated
• Documentation has been built locally and changes look as expected
• Testing comment in the PR documents testing used to verify the changes

Fixes #1068

[xylar]

I wrote out some NetCDF and VTK debugging output with the following code snippet:

        v0 = vertices_on_edge.isel(TWO=0)
        v1 = vertices_on_edge.isel(TWO=1)
        valid = np.logical_and(v0 >= 0, v1 >= 0)
        transport_check = -valid * (bsf_vertex.isel(nVertices=v1) -
                                    bsf_vertex.isel(nVertices=v0))

        velocity_check = 1e6 * transport_check / dv_edge

        vorticity_check = self._compute_vert_integ_vorticity(
            ds_mesh, velocity_check, edge_sign_on_vertex)

        ds_out = xr.Dataset()
        ds_out['barotropicStreamfunction'] = bsf_vertex
        ds_out['verticallyIntegratedVelocity'] = vert_integ_velocity
        ds_out['verticallyIntegratedVorticity'] = vert_integ_vorticity
        ds_out['transportCheck'] = transport_check
        ds_out['verticallyIntegratedVelocityCheck'] = velocity_check
        ds_out['verticallyIntegratedVorticityCheck'] = vorticity_check
        out_filename = f'bsf_{self.subtaskName}.nc'
        ds_out.to_netcdf(out_filename)

        extract_vtk(out_filename, mesh_filename=self.restartFileName,
                    out_dir=f'vtk_bsf_{self.subtaskName}',
                    dimension_list=['vertexDegree=:'],
                    ignore_time=True)

The results indicate that the vertically integrated vorticity is almost identical if computed from the original velocity and from the streamfunction:

The same is not quite true of the vertically integrated normal velocity, which has differences on the order of 1e-3 of the maximum velocity:

This deserves further investigation because it seems like somewhat too large an error to be accounted for by the divergent part of the vertically integrated velocity, which should be very close to zero, especially in a 6-year time average.

[xylar]

Testing

I ran this analysis on a 6-year climatology from a G-case with the SORRMr3 mesh. I had previously noticed that the streamfunction produced for this mesh didn't make sense, particularly in the FRIS area (see #1068). I looked at the transport computed from the BSF compared with the original transport, and they differed by and unacceptably large factors, indicating problems with the methodology even after the fix in #1063.

With the rewrite described above, I'm seeing the results in:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis_testing/BSF/20250125.GMPAS-JRA1p5-DIB-PISMF.TL319_SOwISC12to30E3r3.chrysalis/ocean/index.html#antarctic_extended_bsf

After playing around with colorbars and colormaps quite a bit, I think this linear colormap with an unusual color sequence (the one we often use for resolution plots from Compass) is a better choice than the one from #1063. The fix here increases the maximum value of the BSF (in the middle of the Weddell Gyre) even more than previously, resulting in a completely saturated colormap in Antarctica if we stay with the symmetric-log approach. By effectively having the blue range of the colormap only used in Antarctica, we can have a second linear (divergent) segment that is appropriate for that region.

Thoughts are welcome.

[xylar]

@maltrud and @milenaveneziani, it seems like I messed up again in #1063 and I'll need you to give this another look. Sorry I'm not doing my best work.

This new approach is based on a proper definition of the BSF:

$$\omega H = -\nabla^2 \psi$$

where $\omega H$ is the vertically integrated vorticity and $\nabla^2$ is the Laplacian on the dual mesh (appropriate for vertices).

[xylar]

I want to solve for the divergent contribution to $v H$ and make sure the difference I plotted above is explained by it. If not, I still must be missing something.

[alicebarthel]

Just my 2c but the global plot looks really good to me: good Drake Passage transport, well-defined fronts in Agulhas retroflection and Southern Ocean. I find the new colorbar to be easily readable, well done.

[xylar]

Thanks so much, @alicebarthel!

[xylar]

I also ran the test suite and everything looked fine:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xasay-davis/analysis_testing/chrysalis/fix-bsf-again/

[milenaveneziani]

@xylar : would you be OK with me taking a closer look at this towards the end of this week? My plate is full until at least Th morning, but I am definitely interested/invested in having a correct BSF calculation and thus help with this.

One quick comment I have is about the colorbar: I am personally using the cmocean 'curl' for BSF, but this one looks nice as well. My only suggestion is about flipping the colors. Positive values: orange and reds, negative values: greens and blues (similar to our color association with pos/neg anomalies --> reds/blues). What do you think?

[xylar]

@milenaveneziani, I will wait as long as I can but we are getting close to a new unified release and I need this to be included.

I will see what I can do about the color map but the fact is that we have positive values that are twice as high as negative and green/blue takes up 2/3 of this colormap. So we would need to create a new colormap for what you request.

[xylar]

@milenaveneziani, here's what I came up with for an alternative colormap:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis_testing/BSF/20250125.GMPAS-JRA1p5-DIB-PISMF.TL319_SOwISC12to30E3r3.chrysalis/ocean/index.html

To my eyes, this is not as aesthetic as the one I used above but it does satisfy the desire to have warm colors for positive and cool for negative.

I also tried green instead of blue and purple instead of brown but none of these other combinations looked very good.

Feel free to play around yourself at https://sciviscolor.org/color-moves-app/ and suggest an alternative if you want to take the time. (But I do need to have this done soon for Unified 1.11.0.)

[xylar]

I computed the vertically integrated divergence, and from that I computed the velocity potential. The vertically integrated velocity (VIV) should be the gradient of the potential plus the normal cross the gradient of the barotropic streamfunction. For whatever reason I'm not quite getting this -- the contribution from the velocity potential is actually larger than the difference between the VIV and the piece computed as k cross grad psi. I'd like to understand this better but I'm running out of time I can spare for this work. I think for now I'm content to know that k x grad psi does a very good job now of predicting the VIV field, and the residual is small.

[xylar]

@maltrud, @milenaveneziani and @irenavankova, I've made some more changes to this branch since I asked you to review yesterday, so please take that into account when you review. Sorry for the moving target. I believe I'm done until I get further feedback.

[xylar]

It turns out that zero is not a safe upper bound because this cuts off a few parts of the ocean where the SSH is more than 1 m high (where the middle of the 2-m top layer can be above zero).

[irenavankova]

I made the FRIS focused plots for the 3 versions of the calculations and they are below. I am not sure what to think about the large transport values in the cavity, and all along the coast of Antarctica - what does the nonzero value there imply for the interpretation of clockwise/counterclockwise circulation for positive/negative values?

For comparison, here is a bsf plot from a related regional study of the region, that is easier to interpret:

https://agupubs.onlinelibrary.wiley.com/cms/asset/e4d064df-52a7-4b8f-ba21-2fe278b21c6f/jgrc23944-fig-0008-m.jpg

[xylar]

Thanks for the plots, @irenavankova.

I am not sure what to think about the large transport values in the cavity, and all along the coast of Antarctica - what does the nonzero value there imply for the interpretation of clockwise/counterclockwise circulation for positive/negative values?

There is an arbitrary integration constant in the streamfunction. The only thing that has a physical meaning are differences or gradients in the stream function. So the over all high value (over 150 Sv at the Antarctic coast) tells you that we have that much transport between Antarctica and coastlines to the North (so presumably across Drake Passage). It isn't possible to interpret positive vs. negative signs of the streamfunction as positive vs. negative circulation unless you reference it to the value of the streamfunciton locally (e.g. at the grounding line of FRIS). Unfortunately, that's not possible to do globally.

[xylar]

@irenavankova, following up, it probably makes sense to have a way to plot the streamfunction referenced to the mean value around Antarctica, rather than around the continents north of 45 south. I'll think about how to do that. I think that situation would require a different colormap, too.

[xylar]

@irenavankova, I think we should make the FRIS comparison grid higher resolution. It seems a little blocky in the plots I was looking at when it's 10 km. What do you think?

[irenavankova]

That is fine, as long as there aren't any artifacts popping up from interpolating low res to high res (I have seen that in some field before but here it looks fine)

[xylar]

We always use linear interpolation in MPAS-Analysis and that should always be fine from coarse to fine. The main place it's noticeable is at the coastline (when it doesn't match cartopy) and it isn't a ton better with higher res. But it also improves the streamline a little bit (less jagged).

[xylar]

@irenavankova, I think we might as well add FRIS plots of the BSF to polar regions. Again, what do you think?

[irenavankova]

Sure, I didn't realize that needed to be done explicitly for bsf

[xylar]

Yeah, we don't add any of the projection plots implicitly anywhere. If there are other FRIS plots you want to have on by default, they need to be added as well (but not in this PR).

[irenavankova]

This is looks great, I don't see any issues. Thanks for digging into the inconsistencies and fixing them.

[xylar]

@irenavankova, great, thanks for the review!

[xylar]

@milenaveneziani, please have a look when you have time. It seems the Unified release isn't that eminent -- E3SM-Diags and zstash are both still debugging issues.

[milenaveneziani]

@xylar: this was a lot of work! Thanks for making these big changes: they make sense to me. I also went through the compute_streamfunction code and things look good to me. Great work!

[xylar]

Here are the new plots:

[xylar]

@maltrud, didn't hear from you so will assume you didn't have anything to add.

[xylar]

@milenaveneziani, @alicebarthel and @irenavankova, thanks so much for the input!

[milenaveneziani]

Thanks for posting these!
So, the changes also apply to the subpolar_north_atlantic plot (not just ArcticExtended)? They look great for that region!

[xylar]

Yes, I added the same changes for both regions. But any further changes have to be made in both config sections if we want them to look the same.

[xylar]

Woops! Looks like my last commit didn't get pushed. My Internet connection is bad.

[milenaveneziani]

I think the changes are perfect for the subpolar region. I will test some adjustment to the ArcticExtended region and get back to you, but for now it's great that you added the two sections in the polar_regions config file.

[xylar]

Sounds good.