Merge pull request #28 from FrankErrickson/oxidised_methane

Update CH4 to CO2 oxidation equation

Author: davidanthoff

src/MimiFAIR.jl

@@ -92,7 +92,7 @@ function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year
     set_param!(m, :ch4_cycle, :fossil_frac, gas_fractions.ch4_fossil)
     set_param!(m, :ch4_cycle, :oxidation_frac, 0.61)
     set_param!(m, :ch4_cycle, :mol_weight_CH₄, gas_data[gas_data.gas .== "CH4", :mol_weight][1])
-    set_param!(m, :ch4_cycle, :mol_weight_CO₂, gas_data[gas_data.gas .== "CO2", :mol_weight][1])
+    set_param!(m, :ch4_cycle, :mol_weight_C, gas_data[gas_data.gas .== "C", :mol_weight][1])
     set_param!(m, :ch4_cycle, :emiss2conc_ch4, conversions[conversions.gases .== "CH4", :emiss2conc][1])
 
     # ---- Nitrous Oxide Cycle ---- #
@@ -110,7 +110,6 @@ function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year
     set_param!(m, :co2_cycle, :a, [0.2173, 0.2240, 0.2824, 0.2763])
     set_param!(m, :co2_cycle, :τ_CO₂, [10.0^6, 394.4, 36.54, 4.304])
     set_param!(m, :co2_cycle, :E, rcp_emissions.FossilCO2 .+ rcp_emissions.OtherCO2)
-    set_param!(m, :co2_cycle, :gtc2ppm, conversions[conversions.gases .== "CO2", :emiss2conc][1])
 
     # ---- Other Well-Mixed Greenhouse Gas Cycles ---- #
     set_param!(m, :other_ghg_cycles, :τ_other_ghg, gas_data[findall((in)(other_ghg_names), gas_data.gas), :lifetimes])    
@@ -223,6 +222,7 @@ function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year
     set_param!(m, :NOx_emiss, rcp_emissions.NOx)
     set_param!(m, :fix_pre1850_RCP, true)
     set_param!(m, :mol_weight_N, gas_data[gas_data.gas .== "N", :mol_weight][1])
+    set_param!(m, :gtc2ppm, conversions[conversions.gases .== "CO2", :emiss2conc][1])
 
     # ---------------------------------------------
     # Create connections between Mimi components.

src/components/ch4_cycle.jl

@@ -6,32 +6,37 @@
 
     emiss2conc_ch4    = Parameter()             # Conversion between ppb/ppt concentrations and Mt/kt emissions.
     CH₄_0             = Parameter()             # Initial (pre-industrial) atmospheric methane concentration (ppb).
-    τ_CH₄                 = Parameter()             # Atmospheric (e-folding) lifetime of methane.
+    τ_CH₄             = Parameter()             # Atmospheric (e-folding) lifetime of methane.
+    gtc2ppm           = Parameter()             # Conversion factor between GtC and ppm.
     oxidation_frac    = Parameter()             # Fraction of methane lost through reaction with hydroxyl radical that is converted to carbon dioxide.
     mol_weight_CH₄    = Parameter()             # Molecular mass of methane.
-    mol_weight_CO₂    = Parameter()             # Molecular mass of carbon dioxide.
+    mol_weight_C      = Parameter()             # Molecular mass of carbon.
     fossil_emiss_CH₄  = Parameter(index=[time]) # Fossil-fuel methane emissions (Mt CH₄ yr⁻¹).
     natural_emiss_CH₄ = Parameter(index=[time]) # Natural methane emissions (Mt CH₄ yr⁻¹).
     fossil_frac       = Parameter(index=[time]) # Fraciton of anthropogenic methane attributable to fossil sources.
 
-    CH₄            = Variable(index=[time])  # Atmospheric methane concentration (ppb).
-    oxidised_CH₄   = Variable(index=[time])  # Methane that has been oxidized to carbon dioxide.
-
+    CH₄              = Variable(index=[time])  # Atmospheric methane concentration (ppb).
+    oxidised_CH₄     = Variable(index=[time])  # Methane that has been oxidized to carbon dioxide (ppb).
+    oxidised_CH₄_GtC = Variable(index=[time])  # Methane that has been oxidized to carbon dioxide (GtC)
 
     function run_timestep(p, v, d, t)
 
         # Set initial methane concentration values.
         if is_first(t)
             v.CH₄[t] = p.CH₄_0
             v.oxidised_CH₄[t] = 0.0
+            v.oxidised_CH₄_GtC[t] = 0.0
         else
             # Calculate atmospheric methane concentration.
             emiss_prev = p.fossil_emiss_CH₄[t-1] + p.natural_emiss_CH₄[t]
             emiss_curr = p.fossil_emiss_CH₄[t] + p.natural_emiss_CH₄[t]
             v.CH₄[t] = v.CH₄[t-1] - v.CH₄[t-1] * (1.0 - exp(-1/p.τ_CH₄)) + 0.5 * (emiss_prev + emiss_curr) * (1.0/p.emiss2conc_ch4)
 
             # Calculate carbon dioxide from oxidized methane (bounded below at 0.0).
-            v.oxidised_CH₄[t] = max(((v.CH₄[t-1]-p.CH₄_0) * (1.0 - exp(-1.0/p.τ_CH₄)) * (p.mol_weight_CO₂/p.mol_weight_CH₄ * 0.001 * p.oxidation_frac * p.fossil_frac[t])), 0.0)
+            v.oxidised_CH₄[t] = max(((v.CH₄[t-1]-p.CH₄_0) * (1.0 - exp(-1.0/p.τ_CH₄)) * (p.mol_weight_C/p.mol_weight_CH₄ * 0.001 * p.oxidation_frac * p.fossil_frac[t])), 0.0)
+
+            # Also provide oxidised CH₄ in units of GtC.
+            v.oxidised_CH₄_GtC[t] = v.oxidised_CH₄[t] * p.gtc2ppm
         end
     end
 end