fix to prevent overly efficient iron uptake/scavenging and early onset of aggregation/sinking in warm water

[charliestock]

This pull request addresses the issue that Xiao raised at the COBALT doc&dev meeting several weeks ago. The change in the overall simulation globally is surprisingly modest and generates global chlorophyll simulations that are more skillful than the original after the fix and an additional calibration step (frac_mu_stress changed from 0.25 to 0.4). I described the original issues and the fixes in detail below. It would be good if everyone could run a test for their domains using this code and the suggested parameter change to ensure the suggested fix is robust for regional and global applications.

The Issue: Xiao noted that iron concentrations in COBALTv3 were significantly lower than in COBALTv2 in the equatorial Pacific, and that this was impacting the structure of the nutrient/chlorophyll signal associated with the equatorial upwelling. I tracked this down to two issues:

One was a bug introduced with the photoadaptation option: the optimal P_C_max was saved after multiplying by the temperature dependence (see line 3426). The growth had already been calculated correctly at this point, but P_C_max is later used to calculate the iron uptake (the old line 3520) and the onset of stress dependent mortality (the old line 4091). Both of these assumed that P_C_max was at 0 deg. C, resulting in a double counting of temperature for these factors. This resulted in highly effective phytoplankton iron uptake in warm waters, contributing to the very low equatorial iron concentrations Xiao showed in her presentation. It also accelerated the onset of aggregation and sinking for large phytoplankton. It is notable, however, that this issue primarily impacted that the iron concentration at which a given iron uptake was reached in iron limited regions. It had a limited impact on the overall nitrate and chlorophyll distribution. The figure below, for example, compares the chlorophyll before (left) and after (right) the bug fix 10 years into a 0.5 degree global run with the bug corrected by removing the multiplication by expkT on line 3426:

The differences are modest overall. You can see chlorophyll pattern starts to move toward the equator again in the Pacific, but it is still somewhat disorganized. The chlorophyll at high latitudes increases a bit as sinking and aggregation are delayed relative to before, but the relatively small shifts reflect the secondary importance sinking and aggregation as loss terms and the relatively minor changes in the onset of these terms caused by the bug, and suggest that other factors may be contributing.

The second issue was related to luxury iron uptake by directly sinking phytoplankton in COBALTv3. Iron uptake was linked to P_C_max and the iron uptake limitation (see lines 3520-3521 of the old code). There was no light limitation or consideration of the overall growth rate for iron uptake. This allowed phytoplankton to store iron when other nutrients were limiting (as intended), but also allowed sinking phytoplankton to become quite effective nutrient scavengers as they sank. Adding a light limitation for luxury iron uptake recovers a robust equatorial upwelling pattern in the Pacific:

And adding a positive growth condition for iron uptake subtly reinforced it further (not shown).

The only remaining issue to address is that the mean chlorophyll concentration has increased a bit following the two fixes described above. The impact is modest, but can be addressed by re-tuning the onset of the stress-driven terms so that they now initiate closer to the original calibrated run. I thus re-tuned by increasing the threshold at which stress-driven mortality initiates from the current default of 0.25 to 0.4:

#override frac_mu_stress_Sm = 0.4
#override frac_mu_stress_Md = 0.4
#override frac_mu_stress_Lg = 0.4
#override frac_mu_stress_Di = 0.4

This reduces the overall bias relative to the original simulation. The original chlorophyll field (right) and the one resulting from this pull request (left) are compared below:

The Pacific equatorial upwelling looks much better. There is also some reduction in high-latitude chlorophyll in the Arctic associated with the "frac_mu_stress" change that could improve our Arctic and Bering Sea chlorophyll. I have some additional sensitivities ongoing, but please try these setting in your runs and let me know how it looks!

Thanks again to Xiao for spotting this issue.

[yichengt900]

Thanks, @charliestock, for the detailed rundown of the bug fixes! I’d like to encourage everyone to test this PR in their domains and share any feedback or results here. I will conduct a 0.25° global test on my end as well to check the outcomes. CC @liuxiao37k to keep her updated.

If you're trying to submit a run with this bug fix using the FRE workflow, you can add the following lines to the compilation csh section:

          mv CEFI-regional-MOM6 mom6    
          pushd mom6
          (cd src/ocean_BGC && git checkout expkT_iron_bug/clean_code)
          popd

Also do not forget to add the following in your COBALT_override:

truncate -s 0 $work/INPUT/COBALT_override
cat > $work/INPUT/COBALT_override << COBALT_OVERRIDE_EOF
! Put the contents of your COBALT_override file here.
! For parameters which are not specified in the COBALT_input file, just set the new value, ex:  
! do_case2_mod = True
! For parameters which are specified in the COBALT_input file, specify an override, ex:
!#override do_case2_mod = True
#override frac_mu_stress_Sm = 0.4
#override frac_mu_stress_Md = 0.4
#override frac_mu_stress_Lg = 0.4
#override frac_mu_stress_Di = 0.4                     
COBALT_OVERRIDE_EOF

[andrew-c-ross]

I've got NWA12 with this PR running and slowly rolling into /archive/acr/fre/NWA/2024_09/NWA12_COBALT_2024_09_pr101/gfdl.ncrc6-intel23-prod

[andrew-c-ross]

Just realized that I need to override the frac_mu_stress parameters since the defaults for these are unchanged. I will re-do this and overwrite the archive output mentioned above.

[yichengt900]

Ok here are some preliminary results from the global 0.25° 5-year runs:

Surface chlorophyll (log10 scale)

NPP (Net Primary Production)

Surface NO3

It looks like the bug fixes and the parameter adjustments proposed in this PR have recovered the patterns in the Pacific equatorial upwelling regions. Both the patterns and magnitudes of surface chlorophyll/NO3 and NPP now compare ok with satellite-based observations. I'll await results from a longer run and regional domains, but so far the outcome is promising.

[charliestock]

Thanks for adding these comparisons @yichengt900. It is nice to see the improvement in the subtropical gyres from v2 to v3, and I think the equatorial upwelling now looks better as well. The Southern Ocean looks high, though there are numerous papers suggesting the SeaWIFs is a factor of 2 low. There is a meeting in Tasmania in a year to try and resolve that. Could you share the paths to the 0.25 degree run?

[yichengt900]

@charliestock , I just want to clarify that the V2 results shown in the plots are from SPEAR-COBALT runs, so the forcing is different from our ocean-only runs. It's probably not appropriate to compare them directly, and I should remove them (or just ignore the V2 results) to avoid confusion.

At least we can say the bug fix has improved V3's performance in the equatorial upwelling zones. As for the 0.25-degree results, you can find them here: /archive/ynt/fre/OM4_025/2024_10/OM4_025_COBALTv3_jra55_const/gfdl.ncrc6-intel23-prod/pp.

[yichengt900]

Here’s the new surface chlorophyll plot (mean from 1995-1999). The global annual NPP is now approximately 48.5 Pg/yr. The pattern near the equatorial upwelling zones has improved, though chlorophyll levels near coastal regions are now lower compared to the original results.

[andrew-c-ross]

Managed to get 5 years of this run and transferred so far. Here's the 1997--2019 mean chl before this PR, for context:

Here's chl difference with this PR minus before (only for 1993--1997):

The main differences are an increase around the Mississippi, especially in summer (good), a later peak in the northern part of the domain and a decrease overall for the extreme northern part (mixed good and bad?) and a modest increase along the East Coast in summer and fall (also good).

Here is the seasonal chl with the PR, keeping in mind this is only 5 years of data:

[charliestock]

Thanks for posting these @andrew-c-ross . Hopefully the run has progressed a bit further. How do the seasonal cycles look in the Northeast US Shelf for the phenology work that you have been doing? Sometimes I find that the ratio of new to old is a good way to assess the chlorophyll change. Maybe we could add a difference plot calculated in this way?

@yichengt900, I will assess what happened with the extra growth condition today. This is generally of secondary importance.

[gabyneg]

Here is my analysis of a 5 year run for the hawaiian islands domain, before (1st column) and after (2nd column) the COBLAT v3 PR. (You will notice some boundary artifacts that make it difficult see the influence of the changes alone, but I am currently working on improving those boundaries).

Here, I am showing the chlorophyll profile at station aloha, for the observations (blue), COBLATv3 PR (green), and before the PR (orange), and here you can more clearly see that the changes in the PR have significantly improved the depth of the deep chlorophyll maximum, and the chlorophyll is more similar to observations at station aloha.

And lastly, I am showing nitrate before (1st column), after (2nd column) the COBLAT v3 PR, and in the world ocean atlas product (3rd column), where we see a bit more nitrate on this new COBALTv3 PR, but I need to improve the boundaries to be able to understand the influences that those have on the simulation.

[jessluo]

Here's a quick Look at how this PR (without the growth > 0 condition for iron uptake) performs in ESM4.5 (results shown are years 11-15 of the simulation).

Net primary production:

NPP is showing results as expected, with a large increase in the central equatorial Pacific - and corresponding decreases in a horseshoe pattern around it - as well as further decreases in the mid-latitudes.

The total NPP in these simulations are:
ESM4.5-COBALTv3: 55.2 Pg C y-1
ESM4.5-COBALTv3_PR101: 51.3 Pg C y-1

These are pretty close and well within the range of uncertainty.

Surface chlorophyll:

Surface chlorophyll is also exhibiting the increase in the central equatorial Pacific, but also is showing much more of a coastal effect than NPP. There's a shift in the timing of the chlorophyll bloom in the higher latitudes, with a summer decline in the subpolar regions and increases in the subtropical regions. The effect appears to be most pronounced around 60 N in the late summer. It's probably worth keeping an eye on the chlorophyll climatology and bloom timing, but otherwise this change seems to be quite sensible in the coupled model.

[andrew-c-ross]

Mississippi River area hypoxia is looking pretty good with this PR. Still on the low side, but climatology and variability are solid.

[charliestock]

Thanks all for posting sensitivities. It looks as though this bug fix and enhancement had no objectionable outcomes and/or severe changes. I know some of you are still exploring the best parameter choices to associate with this PR, but I think/hope we are close to being able to bring the code changes in while these continue. Let's plan on a hopefully final discussion at tomorrow's COBALT doc&dev meeting before forging ahead.

[andrew-c-ross]

Here's Gulf of Maine surface chl with 3 different values of the frac_mu_stress parameters. The 0.25 and 0.33 cases are only about halfway done, but the climatologies should be fairly robust.

[yichengt900]

Here are the surface chlorophyll-a (chl-a) results for 1998-2002 from NEP10k, incorporating the new bug fix and frac_mu_stress adjustment of 0.4. The top panel shows results before PR101, and the bottom panel shows results with PR101 (frac=0.4). These findings align with what we observed in our global test runs: chl-a shows better agreement in subtropical regions, and there is a significant reduction in chl-a levels in high-latitude coastal and shelf regions during the AMJ and JAS seasons.

BEFORE PR101 (paper results except year from 1998-2002):

PR101 and frac=0.4 (1998-2002):

[yichengt900]

I’ve checked the DCM (2015–2019 annual mean, units: mg m⁻³) along the equatorial transect from the global 1/4° test run, incorporating the recent bug fix and frac_mu=0.4:

The results show that the DCM is located ~30 m, with a chl-a concentration of ~1.1 mg m⁻³ in the eastern equatorial Pacific..

[yichengt900]

Thank you, @charliestock and team, for the excellent discussions. The preliminary results look promising, though a more detailed check (including Fe model-observation comparisons and DCM checks) might be essential. The code changes look good, and I approve this PR. @charliestock, when you get chance could you also add Xiao's SPEAR-BGC results here to ensure we have a comprehensive reference for any future revisit? Since this PR will change the baseline, I’ll handle that separately in a commit.

[charliestock]

As a final step before merging, I wanted to follow up on Xiao's iron comparisons from yesterday. She showed that the surface iron concentrations after this PR increased relative to COBALTv2 and COBALTv3 before the bug fix. The lowest concentrations were ~0.05 nmol L-1:

Whereas COBALTv3 before this PR and COBALTv2 produced values as low as ~0.005 nmol L-1:

I checked the ESM4.1 documentation paper and it looks like the low-end values in COBALTv3 with the bug fix are generally more realistic. The lowest obs look to be around 0.03:

Obviously we'll want to do more quantitative comparisons in the future, but the PR seems like a step in the right direction.

best,
Charlie