Merge pull request #88 from EcohydrologyTeam/KW_Phosphorus

Phosphorus, N2, POM, and PX

examples/dev_sandbox/DOX_sat_test.py

@@ -25,15 +25,15 @@ def DOs_atm_alpha(
 
 def DOX_sat(
     TwaterK: xr.DataArray,
-    pressure_atm: xr.DataArray,
+    pressure_mb: xr.DataArray,
     pwv: xr.DataArray,
     DOs_atm_alpha: xr.DataArray
 ) -> xr.DataArray:
     """Calculate DO saturation value
 
     Args:
         TwaterK: Water temperature kelvin
-        pressure_atm: Atmospheric pressure (atm)
+        pressure_mb: Atmospheric pressure (mb)
         pwv: Patrial pressure of water vapor (atm)
         DOs_atm_alpha: DO saturation atmospheric correction coefficient
     """
@@ -44,8 +44,8 @@ def DOX_sat(
 
     DOX_sat_uncorrected = -139.34410 + ( 1.575701E05 / TwaterK ) - ( 6.642308E07 / (TwaterK**2.0) ) + ( 1.243800E10 / (TwaterK**3.0) ) - ( 8.621949E11 / (TwaterK**4.0) )
 
-    DOX_sat_corrected = DOX_sat_uncorrected * pressure_atm * \
-        (1 - pwv / pressure_atm) * (1 - DOs_atm_alpha * pressure_atm) / \
+    DOX_sat_corrected = DOX_sat_uncorrected * pressure_mb * \
+        (1 - pwv / pressure_mb) * (1 - DOs_atm_alpha * pressure_mb) / \
         ((1 - pwv) * (1 - DOs_atm_alpha))
 
     print(no_exp)

examples/dev_sandbox/prof_nsm.py

@@ -163,7 +163,7 @@
     'topwidth': 1,
     'slope': 2,
     'shear_velocity': 4,
-    'pressure_atm': 2,
+    'pressure_mb': 2026.5,
     'wind_speed': 4,
     'q_solar': 4,
     'Solid': 1, 

src/clearwater_modules/nsm1/constants.py

@@ -98,7 +98,7 @@ class BalgaeStaticVariables(TypedDict):
     krb_20=0.2,
     kdb_20=0.3,
     mub_max_theta = 1.047,
-    krb_theta = 1.047,
+    krb_theta = 1.06,
     kdb_theta = 1.047,
     b_growth_rate_option=1,
     b_light_limitation_option=1,
@@ -130,11 +130,11 @@ class NitrogenStaticVariables(TypedDict):
     kdnit_20=0.002,
     rnh4_20=0,
     vno3_20=0,
-    knit_theta= 1.047, ## Check values RAS/Kelsey's
-    kon_theta= 1.047,
-    kdnit_theta= 1.047,
+    knit_theta= 1.083,
+    kon_theta= 1.074,
+    kdnit_theta= 1.08,
     rnh4_theta= 1.047,
-    vno3_theta= 1.047,
+    vno3_theta= 1.045,
     KsOxdn=0.1,
     PN=0.5,
     PNb=0.5
@@ -156,10 +156,10 @@ class CarbonStaticVariables(TypedDict):
 DEFAULT_CARBON = CarbonStaticVariables(
     f_pocp = 0.9,
     kdoc_20= 0.01,
-    kdoc_theta = 1.047,
     f_pocb=0.9,
     kpoc_20= 0.005,
     kpoc_theta = 1.047,
+    kdoc_theta = 1.047,
     KsOxmc=1.0,
     pCO2 = 383.0,
     FCO2 = 0.2,
@@ -200,10 +200,12 @@ class N2StaticVariables(TypedDict):
 
 class POMStaticVariables(TypedDict):
     kpom_20: float
+    h2: float
     kpom_theta: float
 
 DEFAULT_POM = POMStaticVariables(
     kpom_20 = 0.1,
+    h2=0.1,
     kpom_theta = 1.047
 )
 
@@ -216,7 +218,7 @@ class PathogenStaticVariables(TypedDict):
 
 DEFAULT_PATHOGEN = PathogenStaticVariables(
     kdx_20=0.8,
-    kdx_theta = 1.047,
+    kdx_theta = 1.07,
     apx=1,
     vx=1
 )
@@ -232,7 +234,7 @@ class PhosphorusStaticVariables(TypedDict):
     kop_20 = 0.1,
     rpo4_20 =0,
     kop_theta = 1.047,
-    rpo4_theta = 1.047,
+    rpo4_theta = 1.074,
     kdpo4 = 0.0,
 )
 
@@ -299,7 +301,7 @@ class GlobalVars(TypedDict):
     topwidth: float
     slope: float
     shear_velocity: float
-    pressure_atm: float
+    pressure_mb: float
     wind_speed: float
     q_solar: float
     Solid: int
@@ -319,24 +321,24 @@ class GlobalVars(TypedDict):
     vs = 999,
     SOD_20 = 999,
     SOD_theta = 999,
+    theta=1.047,
     vb = 0.01,
     fcom = 0.4,
     kaw_20_user = 999,
     kah_20_user = 999,
-    kaw_theta = 1.047,
-    kah_theta = 1.047,
-    hydraulic_reaeration_option = 2,
-    wind_reaeration_option = 2,  
+    kaw_theta = 1.024,
+    kah_theta = 1.024,
+    hydraulic_reaeration_option = 1,
+    wind_reaeration_option = 1,  
     dt = 1,    #TODO Dynamic or static?
     depth = 1.5,         #TODO Dynamic or static?
     TwaterC = 20,
-    theta = 1.047,
     velocity = 1,
     flow = 2,
     topwidth = 1,
     slope = 2,
     shear_velocity = 4,
-    pressure_atm = 2,
+    pressure_mb = 2026.5,
     wind_speed = 4,
     q_solar = 500,
     Solid = 1,

src/clearwater_modules/nsm1/dynamic_variables.py

@@ -1371,14 +1371,6 @@ class Variable(base.Variable):
     process=processes.KHN2_tc
 )
 
-Variable(
-    name='P_wv',
-    long_name='Partial pressure water vapor',
-    units='atm',
-    description='Partial pressure water vapor',
-    use='dynamic',
-    process=processes.P_wv
-)
 
 Variable(
     name='N2sat',

src/clearwater_modules/nsm1/processes.py

@@ -7,13 +7,7 @@
 ############################################ From shared processes
 
 def celsius_to_kelvin(tempc: xr.DataArray) -> xr.DataArray:
-    return tempc + 273.16
-
-
-
-def kelvin_to_celsius(tempk: xr.DataArray) -> xr.DataArray:
-    return tempk - 273.16
-
+    return tempc + 273.15
 
 def arrhenius_correction(
     TwaterC: xr.DataArray,
@@ -1512,6 +1506,8 @@ def NH4_ApGrowth(
 def NH4_AbRespiration(
     use_Balgae: bool,
     rnb: xr.DataArray,
+    Fb: xr.DataArray,
+    depth: xr.DataArray,
     AbRespiration: xr.DataArray,
 
 ) -> xr.DataArray:
@@ -1521,10 +1517,12 @@ def NH4_AbRespiration(
         use_Balgae: true/false to use benthic algae module (unitless),
         rnb: xr.DataArray,
         AbRespiration: Benthic algal respiration rate (g/m^2/d),
+        depth: water depth (m),
+        Fb: Fraction of bottom area for benthic algae (unitless),
     """
     # TODO changed the calculation for respiration from the inital FORTRAN due to conflict with the reference guide
 
-    return xr.where(use_Balgae, rnb * AbRespiration, 0.0 )
+    return xr.where(use_Balgae, (rnb * AbRespiration*Fb)/depth, 0.0 )
 
 def NH4_AbGrowth(
     use_Balgae: bool,
@@ -2022,17 +2020,21 @@ def DIP_ApGrowth(
 def DIP_AbRespiration(
     rpb: xr.DataArray,
     AbRespiration: xr.DataArray,
-    use_Balgae: bool
+    use_Balgae: bool,
+    Fb: xr.DataArray,
+    depth: xr.DataArray
 
 ) -> xr.DataArray :
     """Calculate DIP_AbRespiration: Dissolved inorganic phosphorus released for benthic algal respiration (mg-P/L/d).
 
     Args:
         rpb: Benthic algal P : Benthic algal dry ratio (mg-P/mg-D)
         AbRespiration: Benthic algal respiration rate (g/m^2/d)
-        use_Blgae: true/false to use benthic algae module (t/f)        
+        use_Blgae: true/false to use benthic algae module (t/f) 
+        Fb: Fraction of bottom area available for benthic algal (unitless)
+        depth: water depth (m)       
     """     
-    return xr.where(use_Balgae, rpb * AbRespiration,0)
+    return xr.where(use_Balgae, rpb * Fb * AbRespiration /depth,0)
 
 def DIP_AbGrowth(
     rpb: xr.DataArray,
@@ -2199,7 +2201,7 @@ def POM_algal_settling(
     Ap: xr.DataArray,
     vsap: xr.DataArray,
     rda: xr.DataArray,
-    depth: xr.DataArray,
+    h2: xr.DataArray,
     use_Algae: xr.DataArray
 ) -> xr.DataArray:
     """Calculates the particulate organic matter concentration change due to algal mortality
@@ -2208,10 +2210,10 @@ def POM_algal_settling(
         Ap: Algae concentration (mg/L)
         vsap: Algal settling velocity (m/d)
         rda: Ratio of algal biomass to chlorophyll-a 
-        depth: Depth of water in computation cell (m)
+        h2: active sediment layer thickness (m)
         use_Algae: Option to consider algal kinetics  
     """
-    da: xr.DataArray = xr.where(use_Algae == True, vsap * Ap * rda / depth, 0)
+    da: xr.DataArray = xr.where(use_Algae == True, vsap * Ap * rda / h2, 0)
 
     return da
 
@@ -2234,7 +2236,7 @@ def POM_dissolution(
 def POM_POC_settling(
     POC: xr.DataArray,
     vsoc: xr.DataArray,
-    depth: xr.DataArray,
+    h2: xr.DataArray,
     fcom: xr.DataArray,
     use_POC: xr.DataArray
 ) -> xr.DataArray:
@@ -2243,11 +2245,11 @@ def POM_POC_settling(
     Args:
         POC: Concentration of particulate organic carbon (mg/L)
         vsoc: POC settling velocity (m/d)
-        depth: Depth of water (m)
+        h2: active sediment layer thickness (m)
         fcom: Fraction of carbon in organic matter (mg-C/mg-D) 
         use_POC: Option to consider particulate organic carbon
     """
-    da: xr.DataArray = xr.where(use_POC == True, vsoc * POC / depth / fcom, 0)
+    da: xr.DataArray = xr.where(use_POC == True, vsoc * POC / h2 / fcom, 0)
 
     return da
 
@@ -2257,7 +2259,7 @@ def POM_benthic_algae_mortality(
     kdb_tc: xr.DataArray,
     Fb: xr.DataArray,
     Fw: xr.DataArray,
-    depth: xr.DataArray,
+    h2: xr.DataArray,
     use_Balgae: xr.DataArray
 ) -> xr.DataArray:
     """Calculates particulate organic matter concentration change due to benthic algae mortality
@@ -2267,10 +2269,10 @@ def POM_benthic_algae_mortality(
         kdb_tc: Benthic algae death rate (1/d)
         Fb: Fraction of bottom area available for benthic algae growth
         Fw: Fraction of benthic algae mortality into water column
-        depth: Depth of water in computation cell (m)
+        h2: active sediment layer thickness (m)
         use_Balgae: Option for considering benthic algae in DOC budget (boolean)
     """
-    da: xr.DataArray = xr.where(use_Balgae == True, Ab * kdb_tc * Fb * (1 - Fw) / depth, 0)
+    da: xr.DataArray = xr.where(use_Balgae == True, Ab * kdb_tc * Fb * (1 - Fw) / h2, 0)
 
     return da
 
@@ -2279,16 +2281,16 @@ def POM_benthic_algae_mortality(
 def POM_burial(
     vb: xr.DataArray,
     POM: xr.DataArray,
-    depth: xr.DataArray
+    h2: xr.DataArray
 ) -> xr.DataArray:
     """Calculates particulate organic matter concentration change due to POM burial in the sediments
     
     Args:
         vb: Velocity of burial (m/d)
         POM: POM concentration (mg/L)
-        depth: Depth of water in computation cell (m)
+        h2: active sediment layer thickness (m)
     """
-    return vb * POM / depth
+    return vb * POM / h2 #note removed 365 from FORTRAN
 
 
 
@@ -2898,19 +2900,19 @@ def DOs_atm_alpha(
 
 def DOX_sat(
     TwaterK: xr.DataArray,
-    pressure_atm: xr.DataArray,
+    pressure_mb: xr.DataArray,
     pwv: xr.DataArray,
     DOs_atm_alpha: xr.DataArray
 ) -> xr.DataArray:
     """Calculate DO saturation value
 
     Args:
         TwaterK: Water temperature kelvin
-        pressure_atm: Atmospheric pressure (atm)
+        pressure_mb: Atmospheric pressure (mb)
         pwv: Patrial pressure of water vapor (atm)
         DOs_atm_alpha: DO saturation atmospheric correction coefficient
     """
-    pressure_atm = pressure_atm * 0.000986923
+    pressure_atm = pressure_mb * 0.000986923
 
     DOX_sat_uncorrected = np.exp(-139.34410 + ( 1.575701E05 / TwaterK ) - ( 6.642308E07 / (TwaterK**2.0) ) + ( 1.243800E10 / (TwaterK**3.0) ) - ( 8.621949E11 / (TwaterK**4.0) ))
 
@@ -3465,36 +3467,21 @@ def KHN2_tc(
     return 0.00065 * np.exp(1300.0 * (1.0 / TwaterK - 1 / 298.15))   
 
 
-def P_wv(
-    TwaterK : xr.DataArray,
-) -> xr.DataArray :
-
-    """Calculate partial pressure water vapor (atm)
-
-    Constant values found in documentation
-
-    Args:
-        TwaterK: water temperature kelvin (K)
-
-    """
-    return np.exp(11.8571  - (3840.70 / TwaterK) - (216961.0 / (TwaterK**2)))
-
-
 def N2sat(
     KHN2_tc : xr.DataArray,
-    pressure_atm: xr.DataArray,
-    P_wv: xr.DataArray
+    pressure_mb: xr.DataArray,
+    pwv: xr.DataArray
 ) -> xr.DataArray:
 
     """Calculate N2 at saturation f(Twater and atm pressure) (mg-N/L)
 
     Args:
         KHN2_tc: Henry's law constant (mol/L/atm)
-        pressure_atm: atmosphric pressure in atm (atm)
-        P_wv: Partial pressure of water vapor (atm)
+        pressure_mb: atmosphric pressure in mb (mb)
+        pwv: Partial pressure of water vapor (atm)
     """
 
-    N2sat = 2.8E+4 * KHN2_tc * 0.79 * (pressure_atm - P_wv)  
+    N2sat = 2.8E+4 * KHN2_tc * 0.79 * (pressure_mb*0.000986923 - pwv)  
     N2sat = xr.where(N2sat < 0.0,0.000001,N2sat) #Trap saturation concentration to ensure never negative
 
     return N2sat