Update constants K1 and K2 for carbon chemistry #446

JorgSchwinger · 2024-12-09T10:24:35Z

This PR introduces updated equilibrium constants K1 and K2 for the carbonate system (Waters et al. 2014), which are valid down to low salinities. The limits for temperature and salinity for the carbon chemistry and gas-echange are adjusted such that salinities down to 5 psu are now taken into account. The limiting is done consistently for the sediment and water column using parameters, which values are set in mo_carchm. This PR (together with #432) solves the problems described in #390.

Changes are mainly observed in the Baltic Sea where low salinities occur (see screenshots below, new results to the left). Note the unstable behavior in the previous results, with very high and very low pCO2 in adjacent grid-cells in the northern Baltic Sea. Global time series otherwise indicate only very minor changes.

@jmaerz , please critically review to make sure that there are no problems related to the extended nitrogen cycle when salinities get as low as 5 psu. The model runs can be found here

/projects/NS2980K/schwinger/NOINYOC_T62_tn21_44 (old)
/projects/NS2980K/schwinger/NOINYOC_T62_tn21_45 (new)

TomasTorsvik

@JorgSchwinger - Thanks, looks fine to me!

TomasTorsvik · 2024-12-09T11:08:17Z

hamocc/mo_carchm.F90

-        ah1 = ah2
+      if (abs( erel ) >= eps) then
+        ! make sure [H+] stays between pH 5 and 11
+        ah1 = max(1.0e-11,min(ah2,1.0e-5))


Suggested change

ah1 = max(1.0e-11,min(ah2,1.0e-5))

ah1 = max(min_ah,min(ah2,max_ah))

Maybe introduce this as (private) static parameters in mo_carchm? Reduce risk of divergence between carchm_solve and carchm_solve_dicsat.

Yes, this would make sense, will change accordingly

TomasTorsvik · 2024-12-09T11:08:44Z

hamocc/mo_carchm.F90

-        ah1 = ah2
+      if (abs( erel ) >= eps) then
+        ! make sure [H+] stays between pH 5 and 11
+        ah1 = max(1.0e-11,min(ah2,1.0e-5))


... and here

jmaerz · 2024-12-11T13:02:52Z

hamocc/mo_carchm.F90

@@ -240,11 +245,9 @@ subroutine carchm(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,prho,pglat,omask,psicomo
              ta   = ocetra(i,j,k,inatalkali) / rrho
              ah1  = nathi(i,j,k)

-              call CARCHM_SOLVE(s,tc,ta,sit,pt,K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p,ah1,ac,niter)
+              call CARCHM_SOLVE(psao(i,j,k),tc,ta,sit,pt,K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p,ah1,ac)


Why psao here (and not just s to be consistent with the call above)?

Yes, this didn't make sense. Fixed.

jmaerz · 2024-12-11T13:09:30Z

hamocc/mo_carchm.F90

-            ! solubility of O2 (Weiss, R.F. 1970, Deep-Sea Res., 17, 721-735) for moist air
-            ! at 1 atm; multiplication with oxyco converts to kmol/m^3/atm
+            ! solubility of O2 (Weiss, R.F. 1970, Deep-Sea Res., 17, 721-735) for moist air at
+            ! 1 atm; multiplication with oxyco converts to kmol/m^3/atm, temp=[-1,40], saln=[0,40]


here and below: I would suggest to write temp_min,temp_max , similar for saln (or at least correct the values below: -2 -> -1, 30 -> 40, ...)

The temperature and salinity ranges given in comments indicate the valid range of the fit. This is usually the range for which measurements have been available, taken from the original publication. I have added a comment in carchm_solve to clarify this

jmaerz · 2024-12-11T13:11:21Z

hamocc/mo_carchm.F90

@@ -413,11 +420,12 @@ subroutine carchm(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,prho,pglat,omask,psicomo
              rpp0  = ppao(i,j)/101325.0

              ! calculate correction for non-ideality of CO2 (fugacity coefficient, Weiss and Price 1980)
-              Bvir  = -1636.75 + 12.0408*tk -0.0327957*tk**2 + 0.0000316528*tk**3
-              delta = 57.7-0.118*tk
+              Bvir  = -1636.75 + 12.0408*tk -0.0327957*tk**2 + 0.0000316528*tk**3   ! temp=[-12,47]


temp ranges in comments?

jmaerz · 2024-12-11T13:11:57Z

hamocc/mo_carchm.F90

@@ -765,20 +770,24 @@ subroutine carchm_kequi(temp,saln,prb,Kh0,K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p,Ksp

    ! Kh0 = [CO2]/ fCO2 (fCO2 = fugacity of CO2 in air)
    ! Weiss (1974), note this does not include a correction for the effect of moist air at 
-    ! air-sea interface [mol/kg/atm]
+    ! air-sea interface [mol/kg/atm], temp=[-1,45], saln=[0,45]


jmaerz · 2024-12-11T13:13:34Z

hamocc/mo_carchm.F90

-    ! Millero p.664 (1995) using Mehrbach et al. data on seawater scale
-    K1 = 10**( -1.0 * ( 3670.7 * invtk - 62.008 + 9.7944 * dlogtk - 0.0118 * s + 0.000116 * s2 ) )
-    K2 = 10**( -1.0 * ( 1394.7 * invtk + 4.777 - 0.0184 * s + 0.000118 * s2 ) )
+    ! Waters et al. (2014) on total pH scale, temp=[0,45], saln=[0,45]


here and below - do the comments provide the valid ranges for the formulations? - or should it be what comes from BLOM?

jmaerz · 2024-12-11T13:20:22Z

hamocc/mo_powach.F90

@@ -432,7 +430,7 @@ subroutine powach(kpie,kpje,kpke,kbnd,prho,omask,psao,ptho,lspin)
            K2p = keqb( 9,i,j)
            K3p = keqb(10,i,j)

-            call carchm_solve(saln,c,alk,sit,pt,K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p,ah1,ac,niter)
+            call carchm_solve(psao(i,j,kbo(i,j)),c,alk,sit,pt,K1,K2,Kb,Kw,Ks1,Kf,Ksi,K1p,K2p,K3p,ah1,ac)


why not remaining consistent with what is calculated in the water column? - do we expect the sediment salinity always being between 5 and 40?

I don't understand this comment. It was inconsistent before (water column saln=[25,40], sediment saln=[0,40]), now it is consistent. I think a range of [5,40] should be good enough.

Hi @JorgSchwinger, ok, I got confused through the doubling above of
s = min(saln_max,max(saln_min,psao(i,j,k)))
which is sometimes also calculated before the call to carchm_solve - so my comment above was also due to this confusion.

Yes, I agree the doubling is not so nice. I finally put it in, because we also use carchm_solve in the sediment and I didn't see a much better solution.

Hi @JorgSchwinger , maybe I am missing something, but we could leave out the
s = min(saln_max,max(saln_min,saln))
in carchm_kequi and carchm_solve since one can calculate it before (as now also done for the water column - would be needed to put in again in mo_powach, which was my initial cause of confusion).

jmaerz · 2024-12-11T13:30:42Z

@JorgSchwinger when you refer to the N-cycle: did you mean it as general comment or particularly with respect to the changes that you were putting in?
Further: could you provide the full reference to the updated constants?
And: any incentives to upgrade the whole carbon chemistry system to mocsy (www.geosci-model-dev.net/8/485/2015/
) in line with preferred CMIP protocols? That would also avoid the salinity issues (https://doi.org/10.5194/gmd-10-2169-2017).

JorgSchwinger · 2024-12-11T14:23:38Z

@JorgSchwinger when you refer to the N-cycle: did you mean it as general comment or particularly with respect to the changes that you were putting in?

I meant in particular with respect to the changes in the salinity range (temperature range is actually decreased from [-3,40] to [-1,40]). There are a few formulas (Schmidt number, etc.) that contain salinity as an parameter. Usually these are quite robust against extrapolation, but I just wanted to make sure you are aware of this.

Waters et al. (2014): http://dx.doi.org/10.1016/j.marchem.2014.07.004

This is essentially what mocsy uses, too (the same set of measurements, a slightly different fit, won't make a difference). So regarding a switch to mocsy, I'm not against it, but right now I don't see a large incentive (a lot of work for getting almost identical results).

jmaerz · 2024-12-11T15:02:15Z

I am not sure, if it is needed, but we previously had the minimum temperature of -3 in mo_carchm congruent with the one used in mo_ocprod. Now this deviates from each other. Not sure, if the congruency is needed, though. If you prefer to do it, I would recommend to follow the approach with temp_max,temp_min as for salinity as suggested by Tomas.

JorgSchwinger · 2024-12-12T09:58:42Z

I am not sure, if it is needed, but we previously had the minimum temperature of -3 in mo_carchm congruent with the one used in mo_ocprod. Now this deviates from each other. Not sure, if the congruency is needed, though. If you prefer to do it, I would recommend to follow the approach with temp_max,temp_min as for salinity as suggested by Tomas.

I don't think we need consistency between mo_carchm and mo_ocprod here. The ranges in temperature and salinity solely ensure that the parameterizations are not used out of their range. Not too far out out of their range one should better say, because we have been (and still are) exceeding the range for which observational data has been available for the functional fits of some of the quantities. I think this is not too worrying, since most of functional fits are quite "friendly", they show a smooth behavior outside the data range, they don't go crazy as some polynomial fits do (maybe an exception would be the Wannikhoff Schmidt numbers, which are polynomials in t).

I reduced temp_min from -3 to -1 because none of the functional fits in mo_carchm is actually based on observations for such low temperatures (some are "valid" down to -1, including the Waninkhoff 2014 stuff). I have not checked any ranges for temperature and salinity in mo_ocprod. I'm pretty sure that there is no data down to -3 for e.g. the Eppley temperate curve(?). So we could also adjust this value, but in general this is independent from the changes inmo_carchm, I'd say.

jmaerz · 2024-12-12T10:04:46Z

Hi @JorgSchwinger , I am fine with the decision - just wanted to raise it to be aware of it. I'll have an additional look into the N-cycle stuff tmw (while I assume that there isn't any issue by your changes).

jmaerz · 2024-12-16T11:21:50Z

hamocc/mo_carchm.F90

@@ -25,9 +25,27 @@ module mo_carchm
  public  :: carchm
  public  :: carchm_solve

-  private :: carchm_kequi
-  private :: carchm_solve_dicsat
+  ! Maximum numnber of iterations for carbon chemistry solver


jmaerz

Hi @JorgSchwinger , I am currently don't see a reason, why the nitrogen cycle formulations should now fail, when the range of applicable temperatures has been rather narrowed compared to before. So thanks for making the model more proof for smaller salinities!
Maybe one thing and if you have time: the compiler should handle it, but in the in the routine a number of float values are written as int (without the dot), something one could change...

JorgSchwinger · 2024-12-16T12:38:07Z

OK, I'll fix the integers before merging.

JorgSchwinger requested review from TomasTorsvik and jmaerz December 9, 2024 10:24

TomasTorsvik approved these changes Dec 9, 2024

View reviewed changes

jmaerz reviewed Dec 11, 2024

View reviewed changes

JorgSchwinger force-pushed the feature-new_cc_constants branch from 1176c99 to c4bf05c Compare December 11, 2024 14:33

jmaerz reviewed Dec 16, 2024

View reviewed changes

jmaerz approved these changes Dec 16, 2024

View reviewed changes

JorgSchwinger added 3 commits December 16, 2024 18:09

Update constants K1 and K2 for carbon chemistry

3ff00f3

Adressed issues raised during review

d4d491e

Fix real numbers without decimal point in mo_carchm

74fd875

JorgSchwinger force-pushed the feature-new_cc_constants branch from c4bf05c to 74fd875 Compare December 17, 2024 08:59

JorgSchwinger merged commit 9e4c3e4 into NorESMhub:master Dec 17, 2024
2 of 4 checks passed

JorgSchwinger deleted the feature-new_cc_constants branch December 17, 2024 09:10

JorgSchwinger mentioned this pull request Dec 17, 2024

Testing with hybrid coordinates show spikes in surface pCO2 - problem in Bothnian Sea? #390

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Update constants K1 and K2 for carbon chemistry #446

Update constants K1 and K2 for carbon chemistry #446

JorgSchwinger commented Dec 9, 2024

TomasTorsvik left a comment

TomasTorsvik Dec 9, 2024

JorgSchwinger Dec 11, 2024

TomasTorsvik Dec 9, 2024

jmaerz Dec 11, 2024

JorgSchwinger Dec 11, 2024

jmaerz Dec 11, 2024

JorgSchwinger Dec 11, 2024

jmaerz Dec 11, 2024

jmaerz Dec 11, 2024

jmaerz Dec 11, 2024

jmaerz Dec 11, 2024

JorgSchwinger Dec 11, 2024

jmaerz Dec 11, 2024

JorgSchwinger Dec 11, 2024

jmaerz Dec 11, 2024

jmaerz commented Dec 11, 2024

JorgSchwinger commented Dec 11, 2024

jmaerz commented Dec 11, 2024

JorgSchwinger commented Dec 12, 2024

jmaerz commented Dec 12, 2024

jmaerz Dec 16, 2024

jmaerz left a comment

JorgSchwinger commented Dec 16, 2024

	ah1 = max(1.0e-11,min(ah2,1.0e-5))
	ah1 = max(min_ah,min(ah2,max_ah))

Update constants K1 and K2 for carbon chemistry #446

Update constants K1 and K2 for carbon chemistry #446

Conversation

JorgSchwinger commented Dec 9, 2024

TomasTorsvik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jmaerz commented Dec 11, 2024

JorgSchwinger commented Dec 11, 2024

jmaerz commented Dec 11, 2024

JorgSchwinger commented Dec 12, 2024

jmaerz commented Dec 12, 2024

Choose a reason for hiding this comment

jmaerz left a comment

Choose a reason for hiding this comment

JorgSchwinger commented Dec 16, 2024