bug in iron flux computation

[matt-long]

There is a bug in the iron flux computation in ecosys_forcing.F90. The intended implementation was for black carbon (BC) to have a constant iron fraction and constant bioavailable fraction (of 100%).

The flux computation is here.

We have

forcing_field%field_0d(:,:,iblock) = atm_fe_bioavail_frac(:,:) * &
     (iron_frac_in_atm_fine_dust * atm_fine_dust_flux(:,:,iblock) + &
      iron_frac_in_atm_coarse_dust * atm_coarse_dust_flux(:,:,iblock) + &
      iron_frac_in_atm_bc * atm_black_carbon_flux(:,:,iblock))

! add component from seaice

seaice_fe_bioavail_frac(:,:) = atm_fe_bioavail_frac(:,:)

forcing_field%field_0d(:,:,iblock) = forcing_field%field_0d(:,:,iblock) + seaice_fe_bioavail_frac(:,:) * &
     (iron_frac_in_seaice_dust * seaice_dust_flux(:,:,iblock) + &
      iron_frac_in_seaice_bc * seaice_black_carbon_flux(:,:,iblock))

We should have

atm_bc_fe_bioavail_frac = 1.0_r8
seaice_bc_fe_bioavail_frac = 1.0_r8
...

forcing_field%field_0d(:,:,iblock) = atm_fe_bioavail_frac(:,:) * &
     (iron_frac_in_atm_fine_dust * atm_fine_dust_flux(:,:,iblock) + &
      iron_frac_in_atm_coarse_dust * atm_coarse_dust_flux(:,:,iblock)) + &
      atm_bc_fe_bioavail_frac  * iron_frac_in_atm_bc * atm_black_carbon_flux(:,:,iblock)

! add component from seaice

seaice_fe_bioavail_frac(:,:) = atm_fe_bioavail_frac(:,:)

forcing_field%field_0d(:,:,iblock) = forcing_field%field_0d(:,:,iblock) + seaice_fe_bioavail_frac(:,:) * &
     (iron_frac_in_seaice_dust * seaice_dust_flux(:,:,iblock)) + &
      seaice_bc_fe_bioavail_frac * iron_frac_in_seaice_bc * seaice_black_carbon_flux(:,:,iblock)

The introduction of atm_bc_fe_bioavail_frac and seaice_bc_fe_bioavail_frac parameters is for clarity; this makes it explicit that we are assuming Fe in BC is 100% bioavailable.

I don't think this should be fixed on the CMIP branch as will dramatically change the forcing. I don't see any reason, however, that it should not be fixed moving forward with new runs, including ocean-ice runs.

cc @mnlevy1981, @klindsay28, @kristenkrumhardt

[matt-long]

We should consider a refinement here.

In an email, Douglas Hamilton said,

Iron from combustion is related to black carbon from combustion. But it is not iron in BC. BC is just carbon. Total iron mass emitted from combustion is approximately on parity with black carbon mass emitted. The soluble amount of this iron has been been assumed to be 6%. Therefore soluble iron from combustion is: BC*0.06.

atm_bc_fe_bioavail_frac = 0.06_r8
atm_fe_to_bc_ratio = 1.0_r8
seaice_fe_to_bc_ratio = 1.0_r8
...

forcing_field%field_0d(:,:,iblock) = atm_fe_bioavail_frac(:,:) * &
     (iron_frac_in_atm_fine_dust * atm_fine_dust_flux(:,:,iblock) + &
      iron_frac_in_atm_coarse_dust * atm_coarse_dust_flux(:,:,iblock)) + &
      atm_bc_fe_bioavail_frac  * atm_fe_to_bc_ratio * atm_black_carbon_flux(:,:,iblock)

! add component from seaice

seaice_fe_bioavail_frac(:,:) = atm_fe_bioavail_frac(:,:)

forcing_field%field_0d(:,:,iblock) = forcing_field%field_0d(:,:,iblock) + seaice_fe_bioavail_frac(:,:) * &
     (iron_frac_in_seaice_dust * seaice_dust_flux(:,:,iblock)) + &
      seaice_bc_fe_bioavail_frac * seaice_fe_to_bc_ratio * seaice_black_carbon_flux(:,:,iblock)

Which is mathematically the same, but with different parameter names for conceptual clarity.

[mnlevy1981]

There's some active discussion about this in an on-going email thread -- once consensus is reached there, I'll update this ticket and start working on the agreed-upon solution.

[matt-long]

I think we have some consensus.

From Keith Moore:

Your last version of the code on github is correct, and the suggested parameter name changes are good.

Suggested parameter values:

iron_frac_in_atm_fine_dust = 0.035
iron_frac_in_atm_coarse_dust = 0.035
atm_bc_fe_bioavail_frac = 0.06_r8
atm_fe_to_bc_ratio = 1.0_r8
dust_ratio_thres = 60.0
fe_bioavail_frac_offset = 0.006

I think these are all consistent with our defaults, except the last one, fe_bioavail_frac_offset, which is currently set to 0.01.

[mnlevy1981]

I think these are all consistent with our defaults, except the last one, fe_bioavail_frac_offset, which is currently set to 0.01.

The defaults for dust_ratio_thres and fe_bioavail_frac_offset are a little complicated:

<dust_ratio_thres>55.0</dust_ratio_thres>
<dust_ratio_thres ocn_bgc_config="cesm2.0">60.0</dust_ratio_thres>
<dust_ratio_thres ocn_coupling="partial">66.2818</dust_ratio_thres>

<fe_bioavail_frac_offset>0.01</fe_bioavail_frac_offset>
<fe_bioavail_frac_offset ocn_bgc_config="cesm2.0">0.01</fe_bioavail_frac_offset>
<fe_bioavail_frac_offset ocn_coupling="partial">0.0131875</fe_bioavail_frac_offset>

Note that the default value of ocn_bgc_config is latest; cesm2.0 is the pre-CMIP6 tuning value. If we are changing the default value of these two variables, is it only for the fully-coupled compset or should it change for partially-coupled as well?

[matt-long]

We may have to revisit the partially-coupled values. I think we set those differently to account for nonlinearities in the parameterization, since we force the ocean-ic model with time-mean dust fields.

[mnlevy1981]

I'll update https://github.com/marbl-ecosys/marbl-forcing/blob/master/Fe_aeolian_dep/iron_flux_parameterization.ipynb and recalculate the partially-coupled values

[mnlevy1981]

Given the notebook output

------------------------------------------------------------
orig values
dust_ratio_thres = 60.0
fe_bioavail_frac_offset = 0.01
dust_ratio_to_fe_bioavail_frac = 0.0058823529411764705
------------------------------------------------------------
new values
dust_ratio_thres = 69.00594084513892
fe_bioavail_frac_offset = 0.014675596610025312
dust_ratio_to_fe_bioavail_frac = 0.002729899489804088
------------------------------------------------------------

I have updated defaults (when ocn_coupling="partial"):

 dust_ratio_thres = 69.00594
 dust_ratio_to_fe_bioavail_frac_r = 366.314
 fe_bioavail_frac_offset = 0.0146756

The test suites are currently being re-run, but I don't anticipate any issues. I should note that I added the following lines to the namelist default files:

<dust_ratio_thres ocn_bgc_config="cesm2.0" ocn_coupling="partial">60</dust_ratio_thres>
<fe_bioavail_frac_offset ocn_bgc_config="cesm2.0" ocn_coupling="partial">0.01</fe_bioavail_frac_offset>
<dust_ratio_to_fe_bioavail_frac_r ocn_bgc_config="cesm2.0" ocn_coupling="partial">170.0</dust_ratio_to_fe_bioavail_frac_r>

because I was unsure of how build-namelist would choose between

<dust_ratio_thres ocn_bgc_config="cesm2.0">60.0</dust_ratio_thres>
<dust_ratio_thres ocn_coupling="partial">69.00594</dust_ratio_thres>

if those two could both be considered the "best" match.