Ignore historic TAU patterns + land conversion costs forestry + Bugfix BII + Bugfix NPI/NDC afforestation

CHANGELOG.md

@@ -9,14 +9,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### changed
+- **13_tc_** added switch to ignore historic tau patterns in historic time steps (new default)
+- **39_landconversion** lower costs for expansion of forestry land
 
 ### added
 
 ### removed
+- **core** "removed sets ac_young and ac_mature (no longer needed due to changes in 44_biodiversity)
 
 ### fixed
-
-
+- **32_foresty** BII coefficients for CO2 price driven afforestation
+- **32_foresty** growth curve CO2 price driven afforestation
+- **32_foresty** NPI/NDC afforestation infeasibility
+- **44_biodiversity** ac0 included in pricing of biodiversity loss
 
 ## [4.3.4] - 2021-04-30
 

config/default.cfg

@@ -205,7 +205,10 @@ cfg$gms$s12_hist_interest_hic_noselect <- "0.04"    # def = 0.04
 cfg$gms$tc <- "endo_jun18"              # def = endo_jun18
 
 # * tc cost scenario: low, medium or high
-cfg$gms$c13_tccost <- "medium"
+cfg$gms$c13_tccost <- "medium"		# def = medium
+
+# * ignore historial tau (1) or use it as lower bound (0)
+cfg$gms$s13_ignore_tau_historical <- 1		# def = 1
 
 # ***---------------------    14_yield    --------------------------------------
 # * (biocorrect): cropland TC applied on pasture yields
@@ -443,9 +446,23 @@ cfg$gms$forestry  <- "dynamic_feb21"              # def = dynamic_feb21
 # * afforestation planing horizon (years)
 cfg$gms$s32_planing_horizon <- 50                # def = 50
 
-# * switch for using natveg (0) or plantation (1) growth curves for afforestation
+# * Settings for CO2 price driven afforestation (Growth curve and BII)
+
+# * Growth curve for CO2 price driven afforestation
+# * Switch for using natural vegetation (0) or plantation (1) growth curves towards LPJmL natural 
+# * vegetation carbon density. 
+# * Afforestation following plantation growth curves reflects managed or assistent regrowth,
+# * and might also include non-native species. 
 cfg$gms$s32_aff_plantation <- 0		# def = 0
 
+# * BII coefficient for CO2 price driven afforestation
+# * Switch for using secondary vegetation (0) or timber (1) BII coefficients for CO2 price driven afforestation
+# * The recommend setting is to map the BII coefficient to the choice of the growth curve.
+# * natural vegetation growth curve (0): secondary vegetation BII coefficient (0)
+# * plantation growth curve (1): timber BII coefficient (1)
+# * However, afforestation with plantations (1) could possibly be done in a biodiversity friendly way (0)
+cfg$gms$s32_aff_bii_coeff <- 0		# def = 0
+
 # Afforestation policy
 # * ("none"): no prescribed afforestation
 # * ("npi"): prescribed afforestation based on NPI policies

core/sets.gms

@@ -257,14 +257,6 @@ sets
                     ac105,ac110,ac115,ac120,ac125,ac130,ac135,ac140,ac145,
                     ac150,ac155,acx /
 
-  ac_young(ac) age classes young
-                  / ac5,ac10,ac15,ac20,ac25,ac30 /
-
-  ac_mature(ac) age classes mature
-                     / ac35,ac40,ac45,ac50,ac55,ac60,ac65,ac70,ac75,ac80,ac85,ac90,
-                       ac95,ac100,ac105,ac110,ac115,ac120,ac125,ac130,ac135,ac140,
-                       ac145,ac150,ac155,acx /
-
   ac_est(ac) Dynamic subset of age classes for establishment
 
   ac_sub(ac) Dynamic subset of age classes excluding establishment

modules/13_tc/endo_jun18/input.gms

@@ -6,6 +6,10 @@
 *** |  Contact: magpie@pik-potsdam.de
 
 
+scalars
+ s13_ignore_tau_historical	ignore historial tau (1) or use it as lower bound (0) (binary) / 1 /
+;
+
 parameter fm_tau1995(i) Agricultural land use intensity tau in 1995 (1)
 /
 $ondelim

modules/13_tc/endo_jun18/presolve.gms

@@ -8,7 +8,7 @@
 
 pc13_land(i) = sum(cell(i,j),pcm_land(j,"crop"));
 
-if (sum(sameas(t_past,t),1) = 1,
+if (sum(sameas(t_past,t),1) = 1 AND s13_ignore_tau_historical = 0,
 	vm_tau.lo(i) =    f13_tau_historical(t,i);
 else
 	vm_tau.lo(i) =    pc13_tau(i);

modules/32_forestry/dynamic_feb21/declarations.gms

@@ -54,7 +54,9 @@ parameters
  p32_plantation_contribution(t_ext,i)               Share of roundwood production coming from timber plantations (percent)
  p32_ac_dist_flag(j,ac)                             Distribution flag with inverse weights according to age-classes (1)
  p32_ac_dist(j,ac)                                  Actual share of age-class distribution (1)
- ;
+ p32_bii_coeff(type32,bii_class_secd,potnatveg)		bii coeff (unitless)
+ p32_c_density_ac_fast_forestry(t_all,j,ac)         Carbon densities in plantations based on Braakhekke et al (tC per ha)
+;
 
 positive variables
  vm_cost_fore(i)                                    Forestry costs (Mio USD)
@@ -68,11 +70,11 @@ positive variables
  v32_cost_establishment(i)                          Cost of establishment calculated at the current time step (mio. USD)
  v32_hvarea_forestry(j,ac)                          Harvested area from timber plantations (mio. ha)
  vm_prod_forestry(j,kforestry)                      Production of woody biomass from commercial plantations (mio. tDM per yr)
- ;
+;
 
 variables
  vm_cdr_aff(j,ac,aff_effect)                        Expected bgc (CDR) and local bph effects of afforestation depending on planning horizon (mio. tC)
- ;
+;
 
 equations
  q32_cost_total(i)                                  Total forestry costs constraint (mio. USD)

modules/32_forestry/dynamic_feb21/equations.gms

@@ -99,17 +99,19 @@ sum(ac_est, v32_land(j2,"aff",ac_est)) =l= sum(ac, v32_land(j2,"aff",ac)) - sum(
 
 q32_bv_aff(j2,potnatveg) .. vm_bv(j2,"aff_co2p",potnatveg)
  				  =e=
-          sum(ac_mature, v32_land(j2,"aff",ac_mature)) * fm_bii_coeff("secd_mature",potnatveg) * fm_luh2_side_layers(j2,potnatveg)
-          + sum(ac_young, v32_land(j2,"aff",ac_young)) * fm_bii_coeff("secd_young",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v32_land(j2,"aff",ac)) * 
+          p32_bii_coeff("aff",bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);
 
 q32_bv_ndc(j2,potnatveg) .. vm_bv(j2,"aff_ndc",potnatveg)
  					=e=
-          sum(ac_mature, v32_land(j2,"ndc",ac_mature)) * fm_bii_coeff("secd_mature",potnatveg) * fm_luh2_side_layers(j2,potnatveg)
-        + sum(ac_young, v32_land(j2,"ndc",ac_young)) * fm_bii_coeff("secd_young",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v32_land(j2,"ndc",ac)) * 
+          p32_bii_coeff("ndc",bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);
 
 q32_bv_plant(j2,potnatveg) .. vm_bv(j2,"plant",potnatveg)
  					=e=
-          sum(ac, v32_land(j2,"plant",ac)) * fm_bii_coeff("timber",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v32_land(j2,"plant",ac)) * 
+          p32_bii_coeff("plant",bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);
+
 
 ************************************************************
 **** Timber production related equations in plantations ****

modules/32_forestry/dynamic_feb21/input.gms

@@ -35,10 +35,10 @@ scalars
   s32_free_land_cost              Very high cost for using non existing land for plantation establishment (USD per ha) /1000000/
   s32_max_aff_area                Maximum total global afforestation (mio. ha)    / Inf /
   s32_aff_plantation              Switch for using growth curves for afforestation 0=natveg 1=plantations (1) / 0 /
-  s32_timber_plantation           Switch for using growth curves for timber plantations 0=natveg 1=plantations (1) / 1 /
   s32_tcre_local                  Switch for local (1) or global (0) TRCE factors (1) / 1 /
   s32_forestry_int_rate           Global interest rate for plantations (percent) / 0.05 /
   s32_max_self_suff               Upper ceiling for the self sufficiency used in calculation for establishment decision (1) / 0.8 /
+  s32_aff_bii_coeff               BII coefficent to be used for CO2 price driven afforestation 0=natural vegetation 1=plantation (1) / 0 /
 ;
 
 parameter f32_aff_mask(j) afforestation mask (binary)

modules/32_forestry/dynamic_feb21/preloop.gms

@@ -120,14 +120,14 @@ loop(t_all,
   );
 
 ** Afforestation policies NPI and NDCs
-p32_aff_pol(t,j) = f32_aff_pol(t,j,"%c32_aff_policy%");
+p32_aff_pol(t,j) = round(f32_aff_pol(t,j,"%c32_aff_policy%"),6);
 
 * Calculate the remaining exogenous afforestation with respect to the maximum exogenous target over time.
 * `p32_aff_togo` is used to adjust `s32_max_aff_area` in the constraint `q32_max_aff`.
-p32_aff_togo(t) = sum(j, smax(t2, p32_aff_pol(t2,j)) - p32_aff_pol(t,j));
+p32_aff_togo(t) = smax(t2, sum(j, p32_aff_pol(t2,j))) - sum(j, p32_aff_pol(t,j));
 
 * Adjust the afforestation limit `s32_max_aff_area` upwards, if it is below the exogenous policy target.
-s32_max_aff_area = max(s32_max_aff_area, sum(j, smax(t2, p32_aff_pol(t2,j))) );
+s32_max_aff_area = max(s32_max_aff_area, smax(t2, sum(j, p32_aff_pol(t2,j))) );
 
 p32_cdr_ac(t,j,ac) = 0;
 
@@ -243,5 +243,19 @@ p32_gs_scaling_reg(i)$(f32_gs_relativetarget(i)>0) = f32_gs_relativetarget(i) /
 ** Calibration factors lower than 1 are set to 1
 p32_gs_scaling_reg(i)$(p32_gs_scaling_reg(i) < 1) = 1;
 
-** Update c-densitiy based on calibration factor for growing stocks
+** Save pm_carbon_density_ac_forestry in a parameter before upscaling to FAO growing stocks.
+** This allows to use plantation growth curves for CO2 price driven afforestation.
+p32_c_density_ac_fast_forestry(t_all,j,ac) = pm_carbon_density_ac_forestry(t_all,j,ac,"vegc");
+
+** Update c-density for timber plantations based on calibration factor to match FAO growing stocks
 pm_carbon_density_ac_forestry(t_all,j,ac,"vegc") = pm_carbon_density_ac_forestry(t_all,j,ac,"vegc") * sum(cell(i,j),p32_gs_scaling_reg(i));
+
+** set bii coefficients
+p32_bii_coeff(type32,bii_class_secd,potnatveg) = 0;
+if(s32_aff_bii_coeff = 0,
+ p32_bii_coeff("aff",bii_class_secd,potnatveg) = fm_bii_coeff(bii_class_secd,potnatveg)
+elseif s32_aff_bii_coeff = 1,
+ p32_bii_coeff("aff",bii_class_secd,potnatveg) = fm_bii_coeff("timber",potnatveg)
+);
+p32_bii_coeff("ndc",bii_class_secd,potnatveg) = fm_bii_coeff(bii_class_secd,potnatveg);
+p32_bii_coeff("plant",bii_class_secd,potnatveg) = fm_bii_coeff("timber",potnatveg);

modules/32_forestry/dynamic_feb21/presolve.gms

@@ -16,42 +16,36 @@ ac_sub(ac) = yes$(ord(ac) > (m_yeardiff_forestry(t)/5));
 v32_hvarea_forestry.fx(j,ac_est) = 0;
 v32_land_reduction.fx(j,type32,ac_est) = 0;
 
-** Start ndc **
-* Limit demand for prescribed NPI/NDC afforestation in `p32_aff_pol` if not enough suitable area (`p32_aff_pot`) for afforestation is available.
+** START ndc **
+* calc NPI/NDC afforestation per time step based on forest stock change
+   p32_aff_pol_timestep("y1995",j) = 0;
+   p32_aff_pol_timestep(t,j)$(ord(t)>1) = p32_aff_pol(t,j) - p32_aff_pol(t-1,j);
+* Suitable area (`p32_aff_pot`) for NPI/NDC afforestation
    p32_aff_pot(t,j) = (vm_land.l(j,"crop") - vm_land.lo(j,"crop")) + (vm_land.l(j,"past") - vm_land.lo(j,"past"));
-*correct ndc forest stock based on p32_aff_pot
-  if((ord(t) > 1),
-      p32_aff_pol(t,j)$(p32_aff_pol(t,j) - p32_aff_pol(t-1,j) > p32_aff_pot(t,j)) = p32_aff_pol(t-1,j) + p32_aff_pot(t,j);
-    );
-*calc ndc afforestation per time step based on forest stock changes
-  p32_aff_pol_timestep("y1995",j) = 0;
-  p32_aff_pol_timestep(t,j)$(ord(t)>1) = p32_aff_pol(t,j) - p32_aff_pol(t-1,j);
+* Limit prescribed NPI/NDC afforestation in `p32_aff_pol_timestep` if not enough suitable area (`p32_aff_pot`) for afforestation is available
+   p32_aff_pol_timestep(t,j)$(p32_aff_pol_timestep(t,j) > p32_aff_pot(t,j)) = p32_aff_pot(t,j);
+** END ndc **
 
 if(m_year(t) <= sm_fix_SSP2,
  p32_max_aff_area = Inf;
 else
  p32_max_aff_area = s32_max_aff_area;
 );
-** END ndc **
 
 *' @code
 
 *' Afforestation switch:
-*' 0 = Use natveg carbon densities for afforestation,
-*' 1 = Use plantation carbon densities for afforestation.
+*' 0 = Use natveg growth curve towards LPJmL natural vegetation
+*' 1 = Use plantation growth curve (faster than natveg) towards LPJmL natural vegetation
 if(s32_aff_plantation = 0,
  p32_carbon_density_ac(t,j,"aff",ac,ag_pools) = pm_carbon_density_ac(t,j,ac,ag_pools);
 elseif s32_aff_plantation = 1,
- p32_carbon_density_ac(t,j,"aff",ac,ag_pools) = pm_carbon_density_ac_forestry(t,j,ac,ag_pools);
-);
-*' Timber plantations switch:
-*' 0 = Use natveg carbon densities for timber plantations,
-*' 1 = Use plantation carbon densities for timber plantations.
-if(s32_timber_plantation = 0,
- p32_carbon_density_ac(t,j,"plant",ac,ag_pools) = pm_carbon_density_ac(t,j,ac,ag_pools);
-elseif s32_timber_plantation = 1,
- p32_carbon_density_ac(t,j,"plant",ac,ag_pools) = pm_carbon_density_ac_forestry(t,j,ac,ag_pools);
+ p32_carbon_density_ac(t,j,"aff",ac,ag_pools) = p32_c_density_ac_fast_forestry(t,j,ac);
 );
+
+*' Timber plantations carbon densities:
+p32_carbon_density_ac(t,j,"plant",ac,ag_pools) = pm_carbon_density_ac_forestry(t,j,ac,ag_pools);
+
 *' NDC carbon densities are natveg carbon densities.
 p32_carbon_density_ac(t,j,"ndc",ac,ag_pools) = pm_carbon_density_ac(t,j,ac,ag_pools);
 

modules/35_natveg/dynamic_feb21/equations.gms

@@ -39,13 +39,13 @@
 
  q35_bv_secdforest(j2,potnatveg) .. vm_bv(j2,"secdforest",potnatveg)
  					=e=
-          sum(ac_mature, v35_secdforest(j2,ac_mature)) * fm_bii_coeff("secd_mature",potnatveg) * fm_luh2_side_layers(j2,potnatveg)
-        + sum(ac_young, v35_secdforest(j2,ac_young)) * fm_bii_coeff("secd_young",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v35_secdforest(j2,ac)) * 
+          fm_bii_coeff(bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);
 
  q35_bv_other(j2,potnatveg) .. vm_bv(j2,"other",potnatveg)
  					=e=
-          sum(ac_mature, v35_other(j2,ac_mature)) * fm_bii_coeff("secd_mature",potnatveg) * fm_luh2_side_layers(j2,potnatveg)
-        + sum(ac_young, v35_other(j2,ac_young)) * fm_bii_coeff("secd_young",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v35_other(j2,ac)) * 
+          fm_bii_coeff(bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);
 
 *' NPI/NDC land protection policies are implemented as minium forest land and other land stock.
 

modules/39_landconversion/global_static_aug18/input.gms

@@ -6,6 +6,8 @@
 *** |  Contact: magpie@pik-potsdam.de
 
 scalars
-s39_cost_establish   Cost for establishing new land use (USD05MER per hectare)    /8000/
-s39_cost_clearing    Clearing costs linked to removed biomass (USD05MER per ton C)  /0/
+ s39_cost_establish_crop   Cost for establishing new land use (USD05MER per hectare)    /8000/
+ s39_cost_establish_past   Cost for establishing new land use (USD05MER per hectare)    /8000/
+ s39_cost_establish_forestry   Cost for establishing new land use (USD05MER per hectare)    /1000/
+ s39_cost_clearing    Clearing costs linked to removed biomass (USD05MER per ton C)  /0/
 ;

modules/39_landconversion/global_static_aug18/preloop.gms

@@ -6,7 +6,9 @@
 *** |  Contact: magpie@pik-potsdam.de
 
 i39_cost_establish(land) = 0;
-i39_cost_establish(land_establish39) = s39_cost_establish;
+i39_cost_establish("crop") = s39_cost_establish_crop;
+i39_cost_establish("past") = s39_cost_establish_past;
+i39_cost_establish("forestry") = s39_cost_establish_forestry;
 
 i39_cost_clearing(land) = 0;
 i39_cost_clearing(land_clearing39) = s39_cost_clearing;

modules/39_landconversion/global_static_aug18/realization.gms

@@ -9,8 +9,10 @@
 *' In the global_static_aug18 realization, per hectare land conversion costs are separated into 
 *' costs for expansion of cropland, pasture and forestry (establishment costs) and 
 *' costs for clearing of primary forest, secondary forest and other natural land (clearing costs).
-*' We assume a global cost factor of 8000 USD/ha (static over time) for establishment of managed land. 
-*' For clearing of natural vegetation we assume a global static cost factor of 5 USD/tC (based on @kreidenweis_pasture_2018). 
+*' By default, we assume a global cost factor of 8000 USD/ha (static over time) 
+*' for establishment of cropland and pasture, and 1000 USD/ha for forestry land.
+*' By default, clearing of natural vegetation is not priced. Plausible values range 
+*' between 0-5 USD/tC (based on @kreidenweis_pasture_2018). 
 *'
 *' @limitations Data availability for land conversion costs is very limited.
 

modules/44_biodiversity/bv_btc_mar21/sets.gms

@@ -9,9 +9,20 @@ sets
 	landcover44 land cover classes used in bii calculation
 		/ crop_ann, crop_per, manpast, rangeland, urban, aff_ndc, aff_co2p, primforest, secdforest, other, plant /
 
-	bii_class44 bii land cover classes
+	bii_class44 bii coefficent land cover classes
 		/ crop_ann, crop_per, manpast, rangeland, urban, primary, secd_mature, secd_young, timber /
 
+	bii_class_secd(bii_class44) bii coefficent land cover classes secondary vegetation
+		/ secd_mature, secd_young /
+
+	ac_to_bii_class_secd(ac,bii_class_secd) Mapping between forest ageclasses and bii coefficent land cover classes 
+	/ (ac0,ac5,ac10,ac15,ac20,ac25,ac30)    . (secd_young)
+	  (ac35,ac40,ac45,ac50,ac55,ac60,
+  	   ac65,ac70,ac75,ac80,ac85,ac90,
+	   ac95,ac100,ac105,ac110,ac115,
+	   ac120,ac125,ac130,ac135,ac140,
+	   ac145,ac150,ac155,acx)   			. (secd_mature) /
+
 	price_biodiv44	price paths biodiv loss
 		/ p0,p1,p1_p10,p10,p10_p100,p1_p1000,p10_p10000 /