Merge pull request #1929 from tabeado/fixes

adjust weathering rates naming and weathering calculation

modules/33_CDR/portfolio/datainput.gms

@@ -38,10 +38,14 @@ s33_step = 2.5;
 s33_costs_fix = 0.0422;
 s33_co2_rem_pot = 0.3 * 12/44;       !! default for basalt, for Olivine 1.1
 
-*' carbon removal rate: eqs 2+c1 in strefler, amann et al. (2018): wr = grain surface area based WR (10^-10.53 mol m^-2 s^-1) * molar weight of basalt/forsterite (140.7 g/mol) * 3.155^7 s/a * SSA(gs)
-s33_co2_rem_rate = 10**(-10.53) * 125 * 3.155*10**7 * 69.18*(cm_gs_ew**(-1.24));
-p33_co2_rem_rate("1") = -log(1-s33_co2_rem_rate * 0.94);
-p33_co2_rem_rate("2") = -log(1-s33_co2_rem_rate * 0.29);
+*' rock weathering rate (i.e. fraction of rock weathering per year) at ambient temperature (25 degree C), based on 
+*' eq 2 in strefler, amann et al. (2018): 
+*' wr = grain surface area based weathering rate (10^-10.53 mol m^-2 s^-1) * molar weight of basalt/forsterite (140.7 g/mol) 
+*'      * 3.155^7 s/a * specific surface area(depending on grain size cm_gs_ew)
+s33_rock_weath_rate_ambientT = 10**(-10.53) * 125 * 3.155*10**7 * 69.18*(cm_gs_ew**(-1.24));
+*' rock weathering rate for different climate grades: SI Tab F-1 of strefler, amann et al. (2018)
+p33_rock_weath_rate("1") = s33_rock_weath_rate_ambientT * 0.94;
+p33_rock_weath_rate("2") = s33_rock_weath_rate_ambientT * 0.29;
 
 *' JeS FE demand fit from Thorben: SI D in strefler, amann et al. (2018)
 p33_fedem("weathering", "feels") = 6.62 * cm_gs_ew**(-1.16);

modules/33_CDR/portfolio/declarations.gms

@@ -8,10 +8,10 @@
 scalars
 s33_capture_rate            "CO2 capture rate for capturing emissions, e.g., from burning natural gas" / 0.9 /
 
-s33_co2_rem_pot             "specific carbon removal potential [Gt C per Gt ground rock]"
-s33_co2_rem_rate            "carbon removal rate [fraction of annual reduction of total carbon removal potential]"
-s33_costs_fix               "fixed costs for mining, grinding, spreading [T$/Gt stone]"
-s33_step                    "size of bins in v33_weathering_onfield [Gt stone]"
+s33_co2_rem_pot                 "specific carbon removal potential [Gt C per Gt ground rock]"
+s33_rock_weath_rate_ambientT    "fraction of stone weathering per year at ambient temperature (25 degree C)"
+s33_costs_fix                   "fixed costs for mining, grinding, spreading [T$/Gt stone]"
+s33_step                        "size of bins in v33_weathering_onfield [Gt stone]"
 *JeS* GJ/tCO2 = EJ/Gt CO2 = 44/12 EJ/Gt C.
 
 s33_OAE_efficiency          "the amount of rock required to sequester 1GtC [Gt rock / GtC]"
@@ -22,7 +22,7 @@ s33_OAE_glo_limit           "global limit for OAE [tC / a]"
 parameters
 p33_fedem(all_te,all_enty)               "final energy demand of each technology [EJ/GtC] (for EW the unit is [EJ/Gt stone])"
 p33_LimRock(all_regi)                    "regional share of EW limit [fraction], calculated ex ante for a maximal annual amount of 8 Gt rock in D:\projects\CEMICS\paper_technical\supply_curve_transport_remind_regions.m"
-p33_co2_rem_rate(rlf)                    "carbon removal rate [fraction of annual reduction of total carbon removal potential], multiplied with grade factor"
+p33_rock_weath_rate(rlf)                 "fraction of stone weathering per year depending on climate grade (warm or temperate)"
 p33_EW_upScalingLimit(ttot)              "Annual growth rate limit on upscaling of mining & spreading rocks on fields"
 p33_EW_shortTermEW_Limit(all_regi)       "Limit on 2030 potential for enhanced weathering, defined in Gt rocks, based on % of land on which EW is applied"
 p33_EW_maxShareOfCropland(all_regi)      "Share of cropland that can be used for enhanced weathering. Limits maximum amount of rocks weathering."

modules/33_CDR/portfolio/equations.gms

@@ -122,20 +122,24 @@ q33_EW_capconst(t,regi)..
 *'  Calculation of the total amount of ground rock on the fields in timestep t.
 *'  The first part of the equation describes the decay of the rocks added until that time,
 *'  the rest describes the newly added rocks.
+*'  The amounts already on or newly spread on the fields are multiplied with the total fraction 
+*'  remaining for the next time step, i.e. the fraction not weathering in the time step years. 
+*'  This fraction is generally calculated as (1-p33_rock_weath_rate)**(time_step_years).
+*'  For better solver solution, it is rewritten according to a**b = exp(log(a)*b) as
+*'  exp(log(1-p33_rock_weath_rate) * time_step_years).
+
 ***---------------------------------------------------------------------------
 q33_EW_onfield_tot(ttot,regi,rlf_cz33,rlf)$(ttot.val ge max(2025, cm_startyear))..
     v33_EW_onfield_tot(ttot,regi,rlf_cz33,rlf)
     =e=
-    v33_EW_onfield_tot(ttot-1,regi,rlf_cz33,rlf) * exp(-p33_co2_rem_rate(rlf_cz33) * pm_ts(ttot))
+    v33_EW_onfield_tot(ttot-1,regi,rlf_cz33,rlf) * exp(log(1-p33_rock_weath_rate(rlf_cz33)) * pm_ts(ttot))
     + v33_EW_onfield(ttot-1,regi,rlf_cz33,rlf) * (
         sum(tall$(tall.val le (ttot.val - pm_ts(ttot)/2) and tall.val gt (ttot.val - pm_ts(ttot))),
-            exp(-p33_co2_rem_rate(rlf_cz33) * (ttot.val - tall.val))
+            exp(log(1-p33_rock_weath_rate(rlf_cz33)) * (ttot.val - tall.val)))
         )
-    )
     + v33_EW_onfield(ttot,regi,rlf_cz33,rlf) * (
         sum(tall$(tall.val le ttot.val and tall.val gt (ttot.val - pm_ts(ttot)/2)),
-            exp(-p33_co2_rem_rate(rlf_cz33) * (ttot.val-tall.val))
-        )
+            exp(log(1-p33_rock_weath_rate(rlf_cz33)) * (ttot.val-tall.val)))
     )
 ;
 
@@ -146,8 +150,8 @@ q33_EW_emi(t,regi)..
     vm_emiCdrTeDetail(t,regi, "weathering")
     =e=
     sum((rlf_cz33, rlf),
-        - v33_EW_onfield_tot(t,regi,rlf_cz33,rlf) * s33_co2_rem_pot * (1 - exp(-p33_co2_rem_rate(rlf_cz33)))
-    )
+        - v33_EW_onfield_tot(t,regi,rlf_cz33,rlf) * s33_co2_rem_pot * p33_rock_weath_rate(rlf_cz33)
+        )
     ;
 
 ***---------------------------------------------------------------------------