respiration updates

Author: AlexisRenchon

src/ClimaLand.jl

@@ -418,7 +418,7 @@ import .Soil:
     sublimation_source,
     compute_liquid_influx,
     compute_infiltration_energy_flux
-import .Soil.Biogeochemistry: soil_temperature, soil_moisture
+import .Soil.Biogeochemistry: soil_temperature, soil_moisture, soil_ice
 include("standalone/Snow/Snow.jl")
 using .Snow
 import .Snow: snow_boundary_fluxes!, maximum_snow_cover_fraction!

src/diagnostics/default_diagnostics.jl

@@ -480,7 +480,7 @@ function get_possible_diagnostics(model::EnergyHydrology)
 end
 
 function get_possible_diagnostics(model::SoilCO2Model)
-    return ["sco2", "hr", "scd", "scms", "so2", "soc", "soc_int"]
+    return ["sco2", "hr", "scd", "scms", "so2", "soc", "soc_int", "sco2_ppm"]
 end
 function get_possible_diagnostics(model::CanopyModel)
     diagnostics = [
@@ -617,7 +617,7 @@ function get_short_diagnostics(model::EnergyHydrology)
     return ["swc", "si", "sie", "tsoil", "et"]
 end
 function get_short_diagnostics(model::SoilCO2Model)
-    return ["sco2", "hr", "soc", "soc_int"]
+    return ["sco2", "hr", "soc", "soc_int", "sco2_ppm"]
 end
 function get_short_diagnostics(model::CanopyModel)
     diagnostics = ["gpp", "ct", "lai", "trans", "er", "sif"]

src/diagnostics/define_diagnostics.jl

@@ -979,8 +979,8 @@ function define_diagnostics!(land_model, possible_diags)
         short_name = "sco2",
         long_name = "Soil CO2",
         standard_name = "soil_co2",
-        units = "kg C m^3",
-        comments = "Concentration of CO2 in the porous air space of the soil. (depth resolved)",
+        units = "kg C m^-3",
+        comments = "Soil CO2 carbon mass per soil volume (air + dissolved) (depth resolved).",
         compute! = (out, Y, p, t) -> compute_soilco2!(out, Y, p, t, land_model),
     )
 
@@ -1002,7 +1002,7 @@ function define_diagnostics!(land_model, possible_diags)
         long_name = "Soil Organic Carbon",
         standard_name = "soil_organic_carbon",
         units = "kg C m^-3",
-        comments = "Mass concentration of soil organic carbon. (depth resolved)",
+        comments = "Soil organic carbon mass per soil volume (depth resolved).",
         compute! = (out, Y, p, t) -> compute_soc!(out, Y, p, t, land_model),
     )
     conditional_add_diagnostic_variable!(
@@ -1016,6 +1016,16 @@ function define_diagnostics!(land_model, possible_diags)
             compute_integrated_soc!(out, Y, p, t, land_model),
     )
 
+    # Soil CO2 in ppm (for NEON comparison)
+    add_diagnostic_variable!(
+        short_name = "sco2_ppm",
+        long_name = "Soil Pore Air CO2 Concentration",
+        standard_name = "soil_pore_air_co2_concentration",
+        units = "ppm",
+        comments = "CO2 concentration in soil pore air space, in parts per million by volume. Computed from air-equivalent CO2 concentration using ideal gas law. (depth resolved)",
+        compute! = (out, Y, p, t) -> compute_soilco2_ppm!(out, Y, p, t, land_model),
+    )
+
     # Soil water content
     conditional_add_diagnostic_variable!(
         possible_diags;

src/diagnostics/land_compute_methods.jl

@@ -621,6 +621,34 @@ end # Convert from kg C to mol CO2.
 # To convert from kg C to mol CO2, we need to multiply by:
 # [3.664 kg CO2/ kg C] x [10^3 g CO2/ kg CO2] x [1 mol CO2/44.009 g CO2] = 83.26 mol CO2/kg C
 
+# Soil CO2 in ppm (for comparison with NEON observations)
+# Converts air-equivalent CO2 concentration to ppm using ideal gas law:
+# ppm = (n_CO2 / n_air) × 10^6 = CO2_air_eq * R * T / (M_C * P) × 10^6
+function compute_soilco2_ppm!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::Union{SoilCanopyModel{FT}, LandModel{FT}},
+) where {FT}
+    params = land_model.soilco2.parameters
+    M_C = FT(params.M_C)
+    R = FT(LP.gas_constant(params.earth_param_set))
+
+    CO2_air_eq = p.soilco2.CO2_air_eq  # kg C / m³ air-equivalent
+    T_soil = p.soilco2.T               # K
+    P_sfc = p.drivers.P                # Pa
+
+    if isnothing(out)
+        out = zeros(axes(CO2_air_eq))
+        fill!(field_values(out), NaN)
+        @. out = CO2_air_eq * R * T_soil / (M_C * P_sfc) * FT(1e6)
+        return out
+    else
+        @. out = CO2_air_eq * R * T_soil / (M_C * P_sfc) * FT(1e6)
+    end
+end
+
 ## Other ##
 @diagnostic_compute "sw_albedo" Union{SoilCanopyModel, LandModel} p.α_sfc
 @diagnostic_compute "lw_up" Union{SoilCanopyModel, LandModel} p.LW_u
@@ -721,9 +749,18 @@ function compute_total_runoff!(
         out = zeros(soil.domain.space.surface) # Allocates
         fill!(field_values(out), NaN) # fill with NaNs, even over the ocean
     end
-    out .=
-        Runoff.get_surface_runoff(soil.boundary_conditions.top.runoff, Y, p) .+
-        Runoff.get_subsurface_runoff(soil.boundary_conditions.top.runoff, Y, p)
+    runoff = soil.boundary_conditions.top.runoff
+    surface = Runoff.get_surface_runoff(runoff, Y, p)
+    subsurface = Runoff.get_subsurface_runoff(runoff, Y, p)
+    if !isnothing(surface) && !isnothing(subsurface)
+        out .= surface .+ subsurface
+    elseif !isnothing(surface)
+        out .= surface
+    elseif !isnothing(subsurface)
+        out .= subsurface
+    else
+        fill!(out, 0)
+    end
     return out
 end
 

src/integrated/soil_energy_hydrology_biogeochemistry.jl

@@ -156,6 +156,15 @@ function soil_moisture(driver::PrognosticMet, p, Y, t, z)
     return p.soil.θ_l
 end
 
+"""
+    soil_ice(driver::PrognosticMet, p, Y, t, z)
+
+Returns the prognostic ice content from coupled soil model.
+"""
+function soil_ice(driver::PrognosticMet, p, Y, t, z)
+    return Y.soil.θ_i
+end
+
 
 function ClimaLand.get_drivers(model::LandSoilBiogeochemistry)
     bc = model.soil.boundary_conditions.top

src/simulations/initial_conditions.jl

@@ -97,6 +97,61 @@ function set_soil_initial_conditions!(
     return nothing
 end
 
+"""
+    set_soilco2_initial_conditions!(Y, p, land)
+
+Initialize soil CO2 state to be consistent with atmospheric CO2 and soil state.
+CO2 is stored as carbon mass per soil volume (kg C m⁻³ soil).
+"""
+function set_soilco2_initial_conditions!(Y, p, land)
+    FT = eltype(Y.soilco2.CO2)
+    soil = land.soil
+    soilco2 = land.soilco2
+    params = soilco2.parameters
+
+    θ_l = Y.soil.ϑ_l
+    θ_i = Y.soil.θ_i
+    θ_w = θ_l .+ θ_i
+    ν = soil.parameters.ν
+
+    ρc_s = ClimaLand.Soil.volumetric_heat_capacity.(
+        θ_l,
+        θ_i,
+        soil.parameters.ρc_ds,
+        soil.parameters.earth_param_set,
+    )
+    T_soil = ClimaLand.Soil.temperature_from_ρe_int.(
+        Y.soil.ρe_int,
+        θ_i,
+        ρc_s,
+        soil.parameters.earth_param_set,
+    )
+
+    θ_a = ClimaLand.Soil.Biogeochemistry.volumetric_air_content.(θ_w, ν)
+    R = FT(ClimaLand.Parameters.gas_constant(params.earth_param_set))
+    M_C = FT(params.M_C)
+
+    K_H = @. ClimaLand.Soil.Biogeochemistry.henry_constant(
+        params.K_H_co2_298,
+        params.dln_K_H_co2_dT,
+        T_soil,
+    )
+    β = @. ClimaLand.Soil.Biogeochemistry.beta_gas(K_H, R, T_soil)
+    θ_eff = @. ClimaLand.Soil.Biogeochemistry.effective_porosity(θ_a, θ_l, β)
+
+    @. Y.soilco2.CO2 = θ_eff * p.drivers.c_co2 * p.drivers.P * M_C / (R * T_soil)
+    Y.soilco2.O2_f .= FT(0.21)
+
+    # SOC: exponential decay with depth (same profile as DK-Sor)
+    # SOC(z) = SOC_bot + (SOC_top - SOC_bot) × exp(z / τ_soc)
+    SOC_top = FT(15.0)  # kgC/m³ at surface
+    SOC_bot = FT(0.5)   # kgC/m³ at depth
+    τ_soc = FT(1.0 / log(SOC_top / SOC_bot))  # ~0.294 m
+    z = ClimaCore.Fields.coordinate_field(axes(Y.soilco2.SOC)).z
+    @. Y.soilco2.SOC = SOC_bot + (SOC_top - SOC_bot) * exp(z / τ_soc)
+    return nothing
+end
+
 """
     clip_to_bounds(
     T::FT,
@@ -256,12 +311,6 @@ function make_set_initial_state_from_file(
         if atmos isa ClimaLand.PrescribedAtmosphere
             evaluate!(p.drivers.T, atmos.T, t0)
         end
-        # SoilCO2 IC
-        if !isnothing(land.soilco2)
-            Y.soilco2.CO2 .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
-            Y.soilco2.O2_f .= FT(0.21)    # atmospheric O2 volumetric fraction
-            Y.soilco2.SOC .= FT(5.0)      # default SOC concentration (kg C/m³)
-        end
         # Soil IC
         if enforce_constraints
             # get only the values over land
@@ -280,6 +329,11 @@ function make_set_initial_state_from_file(
             enforce_constraints,
         )
 
+        # SoilCO2 IC (requires soil state)
+        if !isnothing(land.soilco2)
+            set_soilco2_initial_conditions!(Y, p, land)
+        end
+
         # Snow IC
         # Use soil temperature at top to set IC
         soil = land.soil
@@ -396,10 +450,6 @@ function make_set_initial_state_from_file(
             evaluate!(p.drivers.T, atmos.T, t0)
         end
 
-        Y.soilco2.CO2 .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
-        Y.soilco2.O2_f .= FT(0.21)    # atmospheric O2 volumetric fraction
-        Y.soilco2.SOC .= FT(5.0)      # default SOC concentration (kg C/m³)
-
         if enforce_constraints
             # get only the values over land
             T_atmos = Array(parent(p.drivers.T))[:]
@@ -417,6 +467,9 @@ function make_set_initial_state_from_file(
             enforce_constraints,
         )
 
+        # SoilCO2 IC (requires soil state)
+        set_soilco2_initial_conditions!(Y, p, land)
+
         # Canopy IC
         # First determine if leaf water potential is in the file. If so, use
         # that to set the IC; otherwise choose steady state with the soil water.
@@ -551,7 +604,12 @@ function make_set_subseasonal_initial_conditions(
         # Set initial conditions that aren't read in from file
         Y.soilco2.CO2 .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
         Y.soilco2.O2_f .= FT(0.21)    # atmospheric O2 volumetric fraction
-        Y.soilco2.SOC .= FT(5.0)      # default SOC concentration (kg C/m³)
+        # SOC: exponential decay with depth (same profile as DK-Sor)
+        SOC_top = FT(15.0)  # kgC/m³ at surface
+        SOC_bot = FT(0.5)   # kgC/m³ at depth
+        τ_soc = FT(1.0 / log(SOC_top / SOC_bot))  # ~0.294 m
+        z = ClimaCore.Fields.coordinate_field(axes(Y.soilco2.SOC)).z
+        @. Y.soilco2.SOC = SOC_bot + (SOC_top - SOC_bot) * exp(z / τ_soc)
         Y.canopy.hydraulics.ϑ_l .= land.canopy.hydraulics.parameters.ν
 
         # Set snow T first to use in computing snow internal energy from IC file
@@ -781,12 +839,6 @@ function make_set_initial_state_from_atmos_and_parameters(
         if atmos isa ClimaLand.PrescribedAtmosphere
             evaluate!(p.drivers.T, atmos.T, t0)
         end
-        # SoilCO2 IC
-        if !isnothing(land.soilco2)
-            Y.soilco2.CO2 .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
-            Y.soilco2.O2_f .= FT(0.21)    # atmospheric O2 volumetric fraction
-            Y.soilco2.SOC .= FT(5.0)      # default SOC concentration (kg C/m³)
-        end
         (; θ_r, ν, ρc_ds) = land.soil.parameters
         @. Y.soil.ϑ_l = θ_r + (ν - θ_r) / 2
         Y.soil.θ_i .= FT(0.0)
@@ -805,6 +857,18 @@ function make_set_initial_state_from_atmos_and_parameters(
                 earth_param_set,
             )
 
+        # SoilCO2 IC (requires soil state)
+        if !isnothing(land.soilco2)
+            Y.soilco2.CO2 .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
+            Y.soilco2.O2_f .= FT(0.21)    # atmospheric O2 volumetric fraction
+            # SOC: exponential decay with depth (same profile as DK-Sor)
+            SOC_top = FT(15.0)  # kgC/m³ at surface
+            SOC_bot = FT(0.5)   # kgC/m³ at depth
+            τ_soc = FT(1.0 / log(SOC_top / SOC_bot))  # ~0.294 m
+            z = ClimaCore.Fields.coordinate_field(axes(Y.soilco2.SOC)).z
+            @. Y.soilco2.SOC = SOC_bot + (SOC_top - SOC_bot) * exp(z / τ_soc)
+        end
+
         Y.snow.S .= FT(0)
         Y.snow.S_l .= FT(0)
         Y.snow.U .= FT(0)