Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Soil gas diffusivity bugfix #2157

Merged
merged 16 commits into from
Nov 4, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion doc/source/tech_note/Methane/CLM50_Tech_Note_Methane.rst
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,7 @@ molecular free-air diffusion coefficients (:math:`{D}_{0}`
+==========================================================+==========================================================+========================================================+
| Aqueous | 0.9798 + 0.02986\ *T* + 0.0004381\ *T*\ :sup:`2` | 1.172+ 0.03443\ *T* + 0.0005048\ *T*\ :sup:`2` |
+----------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------+
| Gaseous | 0.1875 + 0.0013\ *T* | 0.1759 + 0.0011\ *T* |
| Gaseous | 0.1875 + 0.0013\ *T* | 0.1759 + 0.00117\ *T* |
+----------------------------------------------------------+----------------------------------------------------------+--------------------------------------------------------+

Gaseous diffusivity in soils also depends on the molecular diffusivity,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1199,7 +1199,7 @@ STEM_PROF levdcmp profile for litter C and N
SUPPLEMENT_TO_SMINN_vr levdcmp supplemental N supply gN/m^3/s F
WFPS levdcmp WFPS percent F
anaerobic_frac levdcmp anaerobic_frac m3/m3 F
diffus levdcmp diffusivity m^2/s F
diffus levdcmp diffusivity (from nitrification-denitrification) m^2/s F
fr_WFPS levdcmp fr_WFPS fraction F
n2_n2o_ratio_denit levdcmp n2_n2o_ratio_denit gN/gN F
r_psi levdcmp r_psi m F
Expand Down
44 changes: 34 additions & 10 deletions src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -169,10 +169,11 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
real(r8) :: mu, sigma
real(r8) :: t
real(r8) :: pH(bounds%begc:bounds%endc)
real(r8) :: D0 ! temperature dependence of gaseous diffusion coefficients
!debug-- put these type structure for outing to hist files
real(r8) :: co2diff_con(2) ! diffusion constants for CO2
real(r8) :: eps
real(r8) :: f_a
real(r8) :: fc_air_frac ! Air-filled fraction of soil volume at field capacity
real(r8) :: fc_air_frac_as_frac_porosity ! fc_air_frac as fraction of total porosity
real(r8) :: surface_tension_water ! (J/m^2), Arah and Vinten 1995
real(r8) :: rij_kro_a ! Arah and Vinten 1995
real(r8) :: rij_kro_alpha ! Arah and Vinten 1995
Expand All @@ -183,7 +184,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
real(r8) :: r_max
real(r8) :: r_min(bounds%begc:bounds%endc,1:nlevdecomp)
real(r8) :: ratio_diffusivity_water_gas(bounds%begc:bounds%endc,1:nlevdecomp)
real(r8) :: om_frac
real(r8) :: om_frac, diffus_millingtonquirk, diffus_moldrup
real(r8) :: anaerobic_frac_sat, r_psi_sat, r_min_sat ! scalar values in sat portion for averaging
real(r8) :: organic_max ! organic matter content (kg/m3) where
! soil is assumed to act like peat
Expand Down Expand Up @@ -233,7 +234,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
fmax_denit_carbonsubstrate_vr => soilbiogeochem_nitrogenflux_inst%fmax_denit_carbonsubstrate_vr_col , & ! Output: [real(r8) (:,:) ]
fmax_denit_nitrate_vr => soilbiogeochem_nitrogenflux_inst%fmax_denit_nitrate_vr_col , & ! Output: [real(r8) (:,:) ]
f_denit_base_vr => soilbiogeochem_nitrogenflux_inst%f_denit_base_vr_col , & ! Output: [real(r8) (:,:) ]
diffus => soilbiogeochem_nitrogenflux_inst%diffus_col , & ! Output: [real(r8) (:,:) ] diffusivity (unitless fraction of total diffusivity)
diffus => soilbiogeochem_nitrogenflux_inst%diffus_col , & ! Output: [real(r8) (:,:) ] diffusivity (m2/s)
ratio_k1 => soilbiogeochem_nitrogenflux_inst%ratio_k1_col , & ! Output: [real(r8) (:,:) ]
ratio_no3_co2 => soilbiogeochem_nitrogenflux_inst%ratio_no3_co2_col , & ! Output: [real(r8) (:,:) ]
soil_co2_prod => soilbiogeochem_nitrogenflux_inst%soil_co2_prod_col , & ! Output: [real(r8) (:,:) ] (ug C / g soil / day)
Expand Down Expand Up @@ -267,8 +268,11 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
!---------------- calculate soil anoxia state
! calculate gas diffusivity of soil at field capacity here
! use expression from methane code, but neglect OM for now
f_a = 1._r8 - watfc(c,j) / watsat(c,j)
eps = watsat(c,j)-watfc(c,j) ! Air-filled fraction of total soil volume
fc_air_frac = watsat(c,j)-watfc(c,j) ! theta_a in Riley et al. (2011)
fc_air_frac_as_frac_porosity = 1._r8 - watfc(c,j) / watsat(c,j)
! This calculation of fc_air_frac_as_frac_porosity is algebraically equivalent to
! fc_air_frac/watsat(c,j). In that form, it's easier to see its correspondence
! to theta_a/theta_s in Riley et al. (2011).

! use diffusivity calculation including peat
if (use_lch4) then
Expand All @@ -279,9 +283,20 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
else
om_frac = 1._r8
end if
diffus (c,j) = (d_con_g(2,1) + d_con_g(2,2)*t_soisno(c,j)) * 1.e-4_r8 * &
(om_frac * f_a**(10._r8/3._r8) / watsat(c,j)**2 + &
(1._r8-om_frac) * eps**2 * f_a**(3._r8 / bsw(c,j)) )

! Diffusitivity after Moldrup et al. (2003)
! Eq. 8 in Riley et al. (2011, Biogeosciences)
diffus_moldrup = fc_air_frac**2 * fc_air_frac_as_frac_porosity**(3._r8 / bsw(c,j))

! Diffusivity after Millington & Quirk (1961)
! Eq. 9 in Riley et al. (2011, Biogeosciences)
diffus_millingtonquirk = fc_air_frac**(10._r8/3._r8) / watsat(c,j)**2

! First, get diffusivity as a unitless constant, which is what's needed to
! calculate ratio_k1 below.
diffus (c,j) = &
(om_frac * diffus_millingtonquirk + &
(1._r8-om_frac) * diffus_moldrup )
wwieder marked this conversation as resolved.
Show resolved Hide resolved

! calculate anoxic fraction of soils
! use rijtema and kroess model after Riley et al., 2000
Expand Down Expand Up @@ -373,10 +388,19 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_bgc_soilc, filter_bgc_soilc,
! limit to anoxic fraction of soils
pot_f_denit_vr(c,j) = f_denit_base_vr(c,j) * anaerobic_frac(c,j)

! now calculate the ratio of N2O to N2 from denitrifictaion, following Del Grosso et al., 2000
! now calculate the ratio of N2O to N2 from denitrification, following Del Grosso et al., 2000
! diffusivity constant (figure 6b)
ratio_k1(c,j) = max(1.7_r8, 38.4_r8 - 350._r8 * diffus(c,j))

! Del Grosso et al. (2000) have diffus (their D_FC, "a relative index of gas diffusivity
! through soil assuming a water content of field capacity") as unitless, but diffus history
! field wants m2/s. Here, we use the same theoretical construct as for methane diffusivity
! to convert to m2/s: We multiply by the temperature-dependent free-air diffusion rate.
! NOTE that the coefficients for oxygen are used here; it may be more appropriate to use
! coefficients for the gases being dealt with in this subroutine.
D0 = (d_con_g(2,1) + d_con_g(2,2)*t_soisno(c,j)) * 1.e-4_r8
diffus(c,j) = diffus(c,j) * D0

! ratio function (figure 7c)
if ( soil_co2_prod(c,j) > 1.0e-9_r8 ) then
ratio_no3_co2(c,j) = smin_no3_massdens_vr(c,j) / soil_co2_prod(c,j)
Expand Down
2 changes: 2 additions & 0 deletions src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90
Original file line number Diff line number Diff line change
Expand Up @@ -696,6 +696,8 @@ subroutine InitHistory(this, bounds)
end if

if (use_nitrif_denitrif) then
! NOTE that the calculation for diffusivity here uses coefficients for oxygen.
! It may be more appropriate to use coefficients for N(2)O instead.
this%diffus_col(begc:endc,:) = spval
call hist_addfld_decomp (fname='diffus', units='m^2/s', type2d='levdcmp', &
avgflag='A', long_name='diffusivity', &
Expand Down
Loading