diff --git a/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90 b/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90
index 05c004902e..4ec25bbadc 100644
--- a/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90
+++ b/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90
@@ -171,8 +171,8 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_soilc, &
     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
@@ -183,7 +183,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_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
@@ -267,8 +267,11 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_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) ! f_a is theta_a/theta_s in Riley et al. (2011)
-            eps =  watsat(c,j)-watfc(c,j) ! Air-filled fraction of total soil volume; theta_a in Riley et al. (2011)
+            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
@@ -280,11 +283,19 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_soilc, &
                   om_frac = 1._r8
                end if
 
+               ! 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 * eps**(10._r8/3._r8) / watsat(c,j)**2 + &
-                    (1._r8-om_frac) * eps**2 * f_a**(3._r8 / bsw(c,j)) ) 
+                    (om_frac        * diffus_millingtonquirk + &
+                    (1._r8-om_frac) * diffus_moldrup ) 
 
                ! calculate anoxic fraction of soils
                ! use rijtema and kroess model after Riley et al., 2000