Merge pull request #28 soil surface evaporation

soil surface evaporation

Author: vcorcoro

src/physiology/gasexchange/gasexchange.jl

@@ -10,9 +10,11 @@ include("../../atmosphere/atmosphere.jl")
 @system GasExchange(Atmosphere, BoundaryLayer, Calendar, Stomata, IntercellularSpace, Irradiance, EnergyBalance, C3) begin
     PPFD: photosynthetic_photon_flux_density ~ track(u"μmol/m^2/s" #= Quanta =#, override)
     LAI: leaf_area_index ~ track(override)
-    drought_factor ~ track(override)
+    water_stress ~ track(override)
+    s: senescence_reduction_factor ~ track(override)
 
     A_net_total(A_net, LAI): net_photosynthesis_total => A_net * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)
     A_gross_total(A_gross, LAI): gross_photosynthesis_total => A_gross * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)
     E_total(E, LAI): transpiration_total => E * LAI ~ track(u"mmol/m^2/s" #= H2O =#)
-end
+    R_sw_total(R_sw, LAI): shortwave_radiation_absorbed_total => R_sw * LAI ~ track(u"W/m^2")
+end

src/physiology/gasexchange/stomata.jl

@@ -14,7 +14,8 @@
 end
 
 @system Stomata(StomataBase) begin
-    drought_factor ~ hold()
+    water_stress ~ hold()
+    s: senescence_reduction_factor ~ hold()
     g0 => 0.01 ~ preserve(u"mol/m^2/s/bar" #= H2O =#, parameter)
     g1 => 9.670307198008624 ~ preserve(parameter)
 
@@ -28,10 +29,9 @@ end
         g0 + (g1 * (1-s) * m * (A_net * hs / Cs))
     end ~ track(u"mol/m^2/s/bar" #= H2O =#, min=g0)
 
-    m(drought_factor): transpiration_reduction_factor => begin
+    m(water_stress): transpiration_reduction_factor => begin
         #TODO: implement soil water module
-        drought_factor
+        water_stress
     end ~ track
 
-    s: senescence_reduction_factor ~ track(override)
 end

src/physiology/nitrogen.jl

@@ -73,20 +73,10 @@
         f
     end ~ track(max=1)
 
-    # Relative drought factor from CROPGRO. Used for N_uptake_conversion_factor.
-    "Relative drought factor"
-    drought_factor(SW, minSW, field_capacity, soil_saturation, WP) => begin
-        if SW > field_capacity
-            1.0 - (SW - field_capacity) / (soil_saturation - field_capacity)
-        else
-            1 * ((SW - WP) / (field_capacity - WP))
-        end
-    end ~ track(min=0.1, max=1) 
-
     # Nitrogen uptake conversion factor.
     # How much kg/ha of nitrogen for mg/cm of nitrogen (root)?
-    N_uptake_conversion_factor(RLD, drought_factor, soil_depth) => begin
-        RLD * sqrt(drought_factor) * soil_depth 
+    N_uptake_conversion_factor(RLD, water_stress, soil_depth) => begin
+        RLD * sqrt(water_stress) * soil_depth 
     end ~ track(u"kg*cm/mg/ha")
 
     "Amount of NO3 that stays in soil"
@@ -320,4 +310,4 @@
     # LFSENWT(SENRTE, NMINELF) => SENRTE * NMINELF / 0.16 ~ track(u"g/m^2/hr")
 
     # STSENWT(LFSENWT, petiole_to_leaf) => LFSENWT * petiole_to_leaf ~ track(u"g/m^2/hr")
-end
+end

src/physiology/photosynthesis.jl

@@ -11,16 +11,16 @@ Photosynthesis
 
 #     "Gas exchange model for sunlit leaves"
 
-#     sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width, drought_factor =drought_factor) ~ ::GasExchange
+#     sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width, water_stress =water_stress) ~ ::GasExchange
 
 #     "Gas exchange model for shaded leaves"
-#     shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, drought_factor =drought_factor) ~ ::GasExchange
+#     shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, water_stress =water_stress) ~ ::GasExchange
 
     "Gas exchange model for sunlit leaves"
-    sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width, s=s, drought_factor=drought_factor) ~ ::GasExchange
+    sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width, s=s, water_stress=water_stress) ~ ::GasExchange
 
     "Gas exchange model for shaded leaves"
-    shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, s=s, drought_factor=drought_factor) ~ ::GasExchange
+    shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, s=s, water_stress=water_stress) ~ ::GasExchange
 
 
     #=================
@@ -41,6 +41,11 @@ Photosynthesis
         a + b
     end ~ track(u"mmol/m^2/s" #= H2O =#)
 
+    "Net solar radiation absorbed (sunlit + shaded)"
+    R_sw(a=sunlit_gasexchange.R_sw_total, b=shaded_gasexchange.R_sw_total): radiation_absorbed_J_per_s => begin
+        a + b
+    end ~ track(u"W/m^2")
+
     CO2_weight => 44.0098 ~ preserve(u"g/mol")
     C_weight => 12.0107 ~ preserve(u"g/mol")
     CH2O_weight => 30.031 ~ preserve(u"g/mol")

src/rhizosphere/rhizosphere.jl

@@ -1,7 +1,8 @@
 include("soil.jl")
 include("waterbalance.jl")
+include("soilenergy.jl")
 
 """
 Rhizosphere encompasses all rhizosphere-related systems as mixins.
 """
-@system Rhizosphere(Soil, WaterBalance)
+@system Rhizosphere(Soil, SoilEnergy, WaterBalance)

src/rhizosphere/soilenergy.jl

@@ -0,0 +1,52 @@
+"""
+`SoilEnergy` keeps track of soil energy balance.
+potential soil surface evaporation is determined from a modified Priestly-Taylor model
+(Yan & Shugart, 2010) - doi:10.1029/2009JD013598
+"""
+
+@system SoilEnergy begin
+
+    #=========
+    Parameters
+    =========#
+
+    λ: latent_heat_of_vaporization_at_25 => 44 ~ preserve(u"kJ/mol", parameter) # should be a function of temp?
+    Cp: specific_heat_of_air => 29.3 ~ preserve(u"J/mol/K", parameter)
+
+    pt: priestly_taylor_parameter => 1.35 ~ preserve(parameter) # (Yan and Shugart, 2010)
+
+    ghr: ground_heat_ratio => 0 ~ preserve(parameter) # ASSUMPTION: ground heat flux is negligible
+
+    H2O_weight => 18.01528 ~ preserve(u"g/mol")
+    H2O_density => 997 ~ preserve(u"kg/m^3")
+
+    #=================
+    =================#
+
+    R_soil(solrad, R_sw): radiation_input_to_soil_surface => begin
+        solrad - R_sw
+    end ~ track(u"W/m^2")
+
+    G(R_soil,ghr): ground_heat_flux => begin
+        ghr * R_soil
+    end ~ track(u"W/m^2")
+
+    psychrometric_constant(Cp,P_air,λ) => begin
+        Cp * P_air / λ
+    end ~ track(u"kPa/K")
+
+    λE_soil(pt,RH,R_soil,G,d=VPD_Δ,g=psychrometric_constant): soil_latent_heat_flux => begin
+        pt * RH * d / (d + g) * (R_soil - G)
+    end ~ track(u"W/m^2")
+
+    E0_soil(λE_soil,λ): potential_soil_evaporation_mmol_per_m2_s => begin
+        λE_soil / λ
+    end ~ track(u"mmol/m^2/s")
+
+    # Conversion to mm/hr to match water balance.
+    "potential surface evaporation in mm/hr"
+    potential_surface_evaporation(E0_soil, w=H2O_weight, d=H2O_density) => begin
+        E0_soil * w / d
+    end ~ track(u"mm/hr")
+
+end

src/rhizosphere/waterbalance.jl

@@ -71,15 +71,34 @@ Transpiration
 
     "Intercepted rain"
     rainInterception(interception, rain) => interception * rain ~ track(u"mm/hr")
+
+    "Exess rain and irrigation after transpiration"
+    excessInput(rain, rainInterception, irrigation, transpiration) => begin
+        rain - rainInterception + irrigation - transpiration
+    end ~ track(u"mm/hr", min=0)
+    
+    "Soil surface evaporation modifier"
+    beta(SW, WP, field_capacity) => begin
+        (SW - WP) / (field_capacity - WP)
+    end ~ track
+
+    "Soil surface evaporation"
+    surface_evaporation(potential_surface_evaporation, beta, excessInput) => begin
+        if excessInput > potential_surface_evaporation
+            potential_surface_evaporation
+        else
+            excessInput + beta * (potential_surface_evaporation - excessInput)
+        end
+    end ~ track(u"mm/hr", max=potential_surface_evaporation)
 
-    "Canopy evapotranspiration" # Looks like missing soil surface evaporation 
-    evapotranspiration(transpiration, rainInterception) => begin
-        transpiration + rainInterception
+    "Canopy evapotranspiration"
+    evapotranspiration(transpiration, rainInterception, surface_evaporation) => begin
+        transpiration + rainInterception + surface_evaporation
     end ~ track(u"mm/hr", max=SWhour)
 
     "Potential canopy evapotranspiration"
-    potential_evapotranspiration(transpiration, rainInterception) => begin
-        transpiration + rainInterception
+    potential_evapotranspiration(transpiration, rainInterception, surface_evaporation) => begin
+        transpiration + rainInterception + surface_evaporation
     end ~ track(u"mm/hr")
 
     "Hourly excess soil water"
@@ -166,4 +185,15 @@ Transpiration
             irrigation_rate
         end
     end ~ track(u"mm/hr")
+
+    # Relative drought factor from CROPGRO. Used for N_uptake_conversion_factor.
+    # Captures water stress due to both drought and water logging through reduction in stomatal conductance
+    "Relative water stress factor"
+    water_stress(SW, minSW, field_capacity, soil_saturation, WP) => begin
+        if SW > field_capacity
+            1.0 - (SW - field_capacity) / (soil_saturation - field_capacity)
+        else
+            1 * ((SW - WP) / (field_capacity - WP))
+        end
+    end ~ track(min=0.1, max=1) 
 end