AgroForestry: treecover on cropland and betr

core/macros.gms

@@ -97,12 +97,12 @@ $macro m_linear_cell_data_interpol(output,x,input_x1,input_x2,input_y1,input_y2)
 * macro for simple carbon stocks
 $macro m_carbon_stock(land,carbon_density,item) \
             (land(j2,item) * sum(ct,carbon_density(ct,j2,item,ag_pools)))$(sameas(stockType,"actual")) + \
-            (land(j2,item) * sum(ct,carbon_density(ct,j2,item,ag_pools)))$(sameas(stockType,"actualNoAcEst"));
+            (land(j2,item) * sum(ct,carbon_density(ct,j2,item,ag_pools)))$(sameas(stockType,"actualNoAcEst"))
 
 * macro for carbon stocks with age classes
 $macro m_carbon_stock_ac(land,carbon_density,sets,sets_sub) \
             sum((&&sets), land(j2,&&sets) * sum(ct, carbon_density(ct,j2,&&sets,ag_pools)))$(sameas(stockType,"actual")) + \
-            sum((&&sets_sub), land(j2,&&sets_sub) * sum(ct, carbon_density(ct,j2,&&sets_sub,ag_pools)))$(sameas(stockType,"actualNoAcEst"));
+            sum((&&sets_sub), land(j2,&&sets_sub) * sum(ct, carbon_density(ct,j2,&&sets_sub,ag_pools)))$(sameas(stockType,"actualNoAcEst"))
 
 * macros for peatland module
 $macro m58_LandMerge(land,landForestry,set) \

modules/30_crop/endo_apr21/declarations.gms

@@ -5,6 +5,10 @@
 *** |  MAgPIE License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: [email protected]
 
+scalars
+ s30_shift number of 5-year age-classes corresponding to current time step length (1)
+;
+
 parameters
  p30_avl_cropland(t,j)              Total available land for crop cultivation (mio. ha)
  p30_country_snv_weight(i)          SNV policy country weight per region (1)
@@ -16,13 +20,22 @@ parameters
  i30_avl_cropland_iso(iso)          Available land area for cropland at ISO level (mio. ha)
  p30_snv_scenario_fader(t_all)      SNV scenario fader (1)
  p30_rotation_scenario_fader(t_all) Crop rotation scenario fader (1)
+ p30_treecover(t,j,ac)              Cropland tree cover per age class (mio. ha)
+ pc30_treecover(j,ac)               Cropland tree cover per age class in current time step (mio. ha)
+ p30_carbon_density_ac(t,j,ac,ag_pools)  Carbon density cropland tree cover (tC per ha)
+ p30_treecover_bii_coeff(bii_class_secd,potnatveg)  BII coefficient for cropland treecover (1)
+ p30_treecover_min_shr(t,j)         Minimum share of treecover on cropland (1)
+ p30_treecover_scenario_fader(t_all)  Cropland treecover scenario fader (1)
+ p30_betr_min_shr(t,j)             Minimum share of bioenergy trees on cropland (1)
+ p30_betr_scenario_fader(t_all)    Bioenergy trees scenario fader (1)
 ;
 
 positive variables
 * Fallow land is cropland which is temporarily fallow. Croparea+fallow=cropland
  vm_fallow(j)                       Fallow land (mio. ha)
  vm_area(j,kcr,w)                   Agricultural production area (mio. ha)
  vm_rotation_penalty(i)             Penalty for violating rotational constraints (USD05MER)
+ v30_treecover(j,ac)                Cropland tree cover per age class (mio. ha)
 ;
 
 equations
@@ -36,13 +49,18 @@ equations
  q30_bv_per(j,potnatveg)            Biodiversity value of perennial cropland (mio. ha)
  q30_land_snv(j)                    Land constraint for the SNV policy in cropland areas (mio. ha)
  q30_land_snv_trans(j)              Land transition constraint for SNV policy in cropland areas (mio. ha)
+ q30_treecover(j)                   Land constraint cropland treecover (mio. ha)
+ q30_treecover_est(j,ac)            Land constraint cropland treecover (mio. ha)
+ q30_betr(j)                        Land constraint betr (mio. ha)
+ q30_bv_treecover(j,potnatveg)      Biodiversity value for cropland treecover (Mha)
 ;
 
 *#################### R SECTION START (OUTPUT DECLARATIONS) ####################
 parameters
  ov_fallow(t,j,type)                      Fallow land (mio. ha)
  ov_area(t,j,kcr,w,type)                  Agricultural production area (mio. ha)
  ov_rotation_penalty(t,i,type)            Penalty for violating rotational constraints (USD05MER)
+ ov30_treecover(t,j,ac,type)              Cropland tree cover per age class (mio. ha)
  oq30_cropland(t,j,type)                  Total cropland calculation (mio. ha)
  oq30_avl_cropland(t,j,type)              Available cropland constraint (mio. ha)
  oq30_rotation_max(t,j,crp30,w,type)      Local maximum rotational constraints (mio. ha)
@@ -53,6 +71,10 @@ parameters
  oq30_bv_per(t,j,potnatveg,type)          Biodiversity value of perennial cropland (mio. ha)
  oq30_land_snv(t,j,type)                  Land constraint for the SNV policy in cropland areas (mio. ha)
  oq30_land_snv_trans(t,j,type)            Land transition constraint for SNV policy in cropland areas (mio. ha)
+ oq30_treecover(t,j,type)                 Land constraint cropland treecover (mio. ha)
+ oq30_treecover_est(t,j,ac,type)          Land constraint cropland treecover (mio. ha)
+ oq30_betr(t,j,type)                      Land constraint betr (mio. ha)
+ oq30_bv_treecover(t,j,potnatveg,type)    Biodiversity value for cropland treecover (Mha)
 ;
 *##################### R SECTION END (OUTPUT DECLARATIONS) #####################
 

modules/30_crop/endo_apr21/equations.gms

@@ -10,7 +10,7 @@
 *' the sum of crop and water supply type specific land requirements:
 
  q30_cropland(j2)  ..
-   sum((kcr,w), vm_area(j2,kcr,w)) =e= vm_land(j2,"crop");
+   sum((kcr,w), vm_area(j2,kcr,w)) + vm_fallow(j2) + sum(ac, v30_treecover(j2,ac)) =e= vm_land(j2,"crop");
 
 *' We assume that crop production can only take place on suitable cropland area.
 *' We use a suitability index (SI) map from @zabel_global_2014 to exclude areas
@@ -40,6 +40,19 @@
  q30_land_snv_trans(j2) ..
          sum(land_snv, vm_lu_transitions(j2,"crop",land_snv)) =g= sum(ct, p30_snv_relocation(ct,j2));
 
+ q30_treecover(j2) ..
+ sum(ac, v30_treecover(j2,ac)) =g= sum(ct, p30_treecover_min_shr(ct,j2)) * vm_land(j2,"crop");
+
+ q30_betr(j2) ..
+    sum(w, vm_area(j2,"betr",w)) =g=
+     vm_land(j2,"crop") * sum(ct, p30_betr_min_shr(ct,j2));
+
+*' This constraint distributes additions to forestry land over ac_est,
+*' which depends on the time step length (e.g. ac0 and ac5 for a 10 year time step).
+
+ q30_treecover_est(j2,ac_est) ..
+ v30_treecover(j2,ac_est) =e= sum(ac_est2, v30_treecover(j2,ac_est2))/card(ac_est2);
+
 *' As additional constraints minimum and maximum rotational constraints limit
 *' the placing of crops. On the one hand, these rotational constraints reflect
 *' crop rotations limiting the share a specific crop can cover of the total area
@@ -68,7 +81,8 @@
 
  q30_carbon(j2,ag_pools,stockType) ..
  vm_carbon_stock(j2,"crop",ag_pools,stockType) =e=
-   m_carbon_stock(vm_land,fm_carbon_density,"crop");
+   m_carbon_stock(vm_land,fm_carbon_density,"crop") + 
+   m_carbon_stock_ac(v30_treecover,p30_carbon_density_ac,"ac","ac_sub");
 
 *' The biodiversity value for cropland is calculated separately for annual and perennial crops:
  q30_bv_ann(j2,potnatveg) .. vm_bv(j2,"crop_ann",potnatveg)
@@ -79,3 +93,8 @@
           =e=
           (vm_land(j2,"crop") - sum((crop_ann30,w), vm_area(j2,crop_ann30,w)))
           * fm_bii_coeff("crop_per",potnatveg) * fm_luh2_side_layers(j2,potnatveg);
+
+ q30_bv_treecover(j2,potnatveg) .. vm_bv(j2,"crop_tree",potnatveg)
+          =e=
+          sum(bii_class_secd, sum(ac_to_bii_class_secd(ac,bii_class_secd), v30_treecover(j2,ac)) * 
+          p30_treecover_bii_coeff(bii_class_secd,potnatveg)) * fm_luh2_side_layers(j2,potnatveg);

modules/30_crop/endo_apr21/input.gms

@@ -29,6 +29,17 @@ s30_snv_relocation_data_x1      First reference value in SNV target cropland dat
 s30_snv_relocation_data_x2      Second reference value in SNV target cropland data (1) / 0.5 /
 s30_rotation_scenario_start     Rotation scenario start year      / 2020 /
 s30_rotation_scenario_target    Rotation scenario target year     / 2050 /
+s30_treecover_plantation        Growth curve switch for tree cover on cropland 0=natveg 1=plantations (1) / 0 /
+s30_treecover_bii_coeff         BII coefficent to be used for tree cover on cropland 0=secondary vegetation 1=timber plantations (1) / 0 /
+s30_treecover_min_shr           Share of treecover on cropland (1) / 0 /
+s30_treecover_min_shr_noselect  Share of treecover on cropland (1) / 0 /
+s30_treecover_scenario_start    Cropland treecover scenario start year       / 2020 /
+s30_treecover_scenario_target   Cropland treecover scenario target year      / 2030 /
+s30_treecover_decrease          Cropland treecover can decrease (1) or not (0) / 0 /
+s30_betr_min_shr                Minimum share of bioenergy trees on cropland (1) / 0 /
+s30_betr_min_shr_noselect       Minimum share of bioenergy trees on cropland (1) / 0 /
+s30_betr_scenario_start         Bioenergy trees scenario start year       / 2020 /
+s30_betr_scenario_target        Bioenergy trees scenario target year      / 2050 /
 ;
 
 * Set-switch for countries affected by regional SNV policy
@@ -122,3 +133,13 @@ $ondelim
 $include "./modules/30_crop/endo_apr21/input/SNVTargetCropland.cs3"
 $offdelim
 ;
+
+********* Cropland tree cover *******************************************
+
+parameter f30_treecover(j) Tree cover on cropland in 2019 (mio. ha)
+/
+$ondelim
+$include "./modules/30_crop/endo_apr21/input/CroplandTreecover.cs2"
+$offdelim
+/
+;

modules/30_crop/endo_apr21/input/files

@@ -7,3 +7,5 @@ f30_croparea_w_initialisation.cs3
 avl_cropland_iso.cs3
 SNVTargetCropland.cs3
 SNVTargetCropland_0.5.mz
+CroplandTreecover.cs2
+CroplandTreecover_0.5.mz

modules/30_crop/endo_apr21/postsolve.gms

@@ -5,11 +5,13 @@
 *** |  MAgPIE License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: [email protected]
 
+pc30_treecover(j,ac) = v30_treecover.l(j,ac);
 
 *#################### R SECTION START (OUTPUT DEFINITIONS) #####################
  ov_fallow(t,j,"marginal")                      = vm_fallow.m(j);
  ov_area(t,j,kcr,w,"marginal")                  = vm_area.m(j,kcr,w);
  ov_rotation_penalty(t,i,"marginal")            = vm_rotation_penalty.m(i);
+ ov30_treecover(t,j,ac,"marginal")              = v30_treecover.m(j,ac);
  oq30_cropland(t,j,"marginal")                  = q30_cropland.m(j);
  oq30_avl_cropland(t,j,"marginal")              = q30_avl_cropland.m(j);
  oq30_rotation_max(t,j,crp30,w,"marginal")      = q30_rotation_max.m(j,crp30,w);
@@ -20,9 +22,14 @@
  oq30_bv_per(t,j,potnatveg,"marginal")          = q30_bv_per.m(j,potnatveg);
  oq30_land_snv(t,j,"marginal")                  = q30_land_snv.m(j);
  oq30_land_snv_trans(t,j,"marginal")            = q30_land_snv_trans.m(j);
+ oq30_treecover(t,j,"marginal")                 = q30_treecover.m(j);
+ oq30_treecover_est(t,j,ac,"marginal")          = q30_treecover_est.m(j,ac);
+ oq30_betr(t,j,"marginal")                      = q30_betr.m(j);
+ oq30_bv_treecover(t,j,potnatveg,"marginal")    = q30_bv_treecover.m(j,potnatveg);
  ov_fallow(t,j,"level")                         = vm_fallow.l(j);
  ov_area(t,j,kcr,w,"level")                     = vm_area.l(j,kcr,w);
  ov_rotation_penalty(t,i,"level")               = vm_rotation_penalty.l(i);
+ ov30_treecover(t,j,ac,"level")                 = v30_treecover.l(j,ac);
  oq30_cropland(t,j,"level")                     = q30_cropland.l(j);
  oq30_avl_cropland(t,j,"level")                 = q30_avl_cropland.l(j);
  oq30_rotation_max(t,j,crp30,w,"level")         = q30_rotation_max.l(j,crp30,w);
@@ -33,9 +40,14 @@
  oq30_bv_per(t,j,potnatveg,"level")             = q30_bv_per.l(j,potnatveg);
  oq30_land_snv(t,j,"level")                     = q30_land_snv.l(j);
  oq30_land_snv_trans(t,j,"level")               = q30_land_snv_trans.l(j);
+ oq30_treecover(t,j,"level")                    = q30_treecover.l(j);
+ oq30_treecover_est(t,j,ac,"level")             = q30_treecover_est.l(j,ac);
+ oq30_betr(t,j,"level")                         = q30_betr.l(j);
+ oq30_bv_treecover(t,j,potnatveg,"level")       = q30_bv_treecover.l(j,potnatveg);
  ov_fallow(t,j,"upper")                         = vm_fallow.up(j);
  ov_area(t,j,kcr,w,"upper")                     = vm_area.up(j,kcr,w);
  ov_rotation_penalty(t,i,"upper")               = vm_rotation_penalty.up(i);
+ ov30_treecover(t,j,ac,"upper")                 = v30_treecover.up(j,ac);
  oq30_cropland(t,j,"upper")                     = q30_cropland.up(j);
  oq30_avl_cropland(t,j,"upper")                 = q30_avl_cropland.up(j);
  oq30_rotation_max(t,j,crp30,w,"upper")         = q30_rotation_max.up(j,crp30,w);
@@ -46,9 +58,14 @@
  oq30_bv_per(t,j,potnatveg,"upper")             = q30_bv_per.up(j,potnatveg);
  oq30_land_snv(t,j,"upper")                     = q30_land_snv.up(j);
  oq30_land_snv_trans(t,j,"upper")               = q30_land_snv_trans.up(j);
+ oq30_treecover(t,j,"upper")                    = q30_treecover.up(j);
+ oq30_treecover_est(t,j,ac,"upper")             = q30_treecover_est.up(j,ac);
+ oq30_betr(t,j,"upper")                         = q30_betr.up(j);
+ oq30_bv_treecover(t,j,potnatveg,"upper")       = q30_bv_treecover.up(j,potnatveg);
  ov_fallow(t,j,"lower")                         = vm_fallow.lo(j);
  ov_area(t,j,kcr,w,"lower")                     = vm_area.lo(j,kcr,w);
  ov_rotation_penalty(t,i,"lower")               = vm_rotation_penalty.lo(i);
+ ov30_treecover(t,j,ac,"lower")                 = v30_treecover.lo(j,ac);
  oq30_cropland(t,j,"lower")                     = q30_cropland.lo(j);
  oq30_avl_cropland(t,j,"lower")                 = q30_avl_cropland.lo(j);
  oq30_rotation_max(t,j,crp30,w,"lower")         = q30_rotation_max.lo(j,crp30,w);
@@ -59,4 +76,8 @@
  oq30_bv_per(t,j,potnatveg,"lower")             = q30_bv_per.lo(j,potnatveg);
  oq30_land_snv(t,j,"lower")                     = q30_land_snv.lo(j);
  oq30_land_snv_trans(t,j,"lower")               = q30_land_snv_trans.lo(j);
+ oq30_treecover(t,j,"lower")                    = q30_treecover.lo(j);
+ oq30_treecover_est(t,j,ac,"lower")             = q30_treecover_est.lo(j,ac);
+ oq30_betr(t,j,"lower")                         = q30_betr.lo(j);
+ oq30_bv_treecover(t,j,potnatveg,"lower")       = q30_bv_treecover.lo(j,potnatveg);
 *##################### R SECTION END (OUTPUT DEFINITIONS) ######################

modules/30_crop/endo_apr21/preloop.gms

@@ -9,6 +9,8 @@
 * sigmoidal interpolation between start year and target year
 m_sigmoid_time_interpol(p30_snv_scenario_fader,s30_snv_scenario_start,s30_snv_scenario_target,0,1);
 m_sigmoid_time_interpol(p30_rotation_scenario_fader,s30_rotation_scenario_start,s30_rotation_scenario_target,0,1);
+m_sigmoid_time_interpol(p30_treecover_scenario_fader,s30_treecover_scenario_start,s30_treecover_scenario_target,0,1);
+m_sigmoid_time_interpol(p30_betr_scenario_fader,s30_betr_scenario_start,s30_betr_scenario_target,0,1);
 
 * linear interpolation to estimate the cropland that
 * requires relocation due to SNV policy
@@ -37,3 +39,23 @@ p30_country_dummy(policy_countries30) = 1;
 * Countries are weighted by available cropland area.
 i30_avl_cropland_iso(iso) = f30_avl_cropland_iso(iso,"%c30_marginal_land%");
 p30_country_snv_weight(i) = sum(i_to_iso(i,iso), p30_country_dummy(iso) * i30_avl_cropland_iso(iso)) / sum(i_to_iso(i,iso), i30_avl_cropland_iso(iso));
+
+* Initial tree cover on cropland is assumed to be equally distributed among all age-classes
+pc30_treecover(j,ac) = f30_treecover(j)/card(ac);
+
+*' Switch for tree cover on cropland:
+*' 0 = Use natveg growth curve towards LPJmL natural vegetation
+*' 1 = Use plantation growth curve (faster than natveg) towards LPJmL natural vegetation
+if(s30_treecover_plantation = 0,
+ p30_carbon_density_ac(t,j,ac,ag_pools) = pm_carbon_density_ac(t,j,ac,ag_pools);
+elseif s30_treecover_plantation = 1,
+ p30_carbon_density_ac(t,j,ac,ag_pools) = pm_carbon_density_ac_forestry(t,j,ac,ag_pools);
+);
+
+** set bii coefficients
+p30_treecover_bii_coeff(bii_class_secd,potnatveg) = 0;
+if(s30_treecover_bii_coeff = 0,
+ p30_treecover_bii_coeff(bii_class_secd,potnatveg) = fm_bii_coeff(bii_class_secd,potnatveg)
+elseif s30_treecover_bii_coeff = 1,
+ p30_treecover_bii_coeff(bii_class_secd,potnatveg) = fm_bii_coeff("timber",potnatveg)
+);

modules/30_crop/endo_apr21/presolve.gms

@@ -50,3 +50,43 @@ p30_snv_relocation(t,j)$(p30_snv_relocation(t, j) > p30_max_snv_relocation(t,j))
 *' Area potentially available for cropping
 p30_avl_cropland(t,j) = f30_avl_cropland(j,"%c30_marginal_land%") * (1 - p30_snv_shr(t,j));
 *' @stop
+
+* Growth of trees on cropland is modelled by shifting age-classes according to time step length.
+s30_shift = m_timestep_length_forestry/5;
+* example: ac10 in t = ac5 (ac10-1) in t-1 for a 5 yr time step (s30_shift = 1)
+    p30_treecover(t,j,ac)$(ord(ac) > s30_shift) = pc30_treecover(j,ac-s30_shift);
+* account for cases at the end of the age class set (s30_shift > 1) which are not shifted by the above calculation
+    p30_treecover(t,j,"acx") = p30_treecover(t,j,"acx")
+                  + sum(ac$(ord(ac) > card(ac)-s30_shift), pc30_treecover(j,ac));
+
+pc30_treecover(j,ac) = p30_treecover(t,j,ac);
+v30_treecover.l(j,ac) = p30_treecover(t,j,ac);
+
+p30_treecover_min_shr(t,j) = p30_treecover_scenario_fader(t) *
+  (s30_treecover_min_shr * sum(cell(i,j), p30_country_snv_weight(i))
+  + s30_treecover_min_shr_noselect * sum(cell(i,j), 1-p30_country_snv_weight(i)));
+
+*define ac_est and ac_sub
+ac_est(ac) = no;
+ac_est(ac) = yes$(ord(ac) <= (m_yeardiff_forestry(t)/5));
+
+ac_sub(ac) = no;
+ac_sub(ac) = yes$(ord(ac) > (m_yeardiff_forestry(t)/5));
+
+* fix before s30_treecover_scenario_start to pc30_treecover(j,ac)?
+if(m_year(t) <= s30_treecover_scenario_start,
+ v30_treecover.fx(j,ac) = pc30_treecover(j,ac);
+else
+ v30_treecover.lo(j,ac_est) = 0;
+ v30_treecover.up(j,ac_est) = Inf;
+ if(s30_treecover_decrease = 1,
+  v30_treecover.lo(j,ac_sub) = 0;
+  v30_treecover.up(j,ac_sub) = pc30_treecover(j,ac_sub);
+ else
+  v30_treecover.fx(j,ac_sub) = pc30_treecover(j,ac_sub);  
+ );
+);
+
+p30_betr_min_shr(t,j) = p30_betr_scenario_fader(t) *
+  (s30_betr_min_shr * sum(cell(i,j), p30_country_snv_weight(i))
+  + s30_betr_min_shr_noselect * sum(cell(i,j), 1-p30_country_snv_weight(i)));

modules/32_forestry/dynamic_feb21/presolve.gms

@@ -5,13 +5,6 @@
 *** |  MAgPIE License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: [email protected]
 
-*define ac_est and ac_sub
-ac_est(ac) = no;
-ac_est(ac) = yes$(ord(ac) <= (m_yeardiff_forestry(t)/5));
-
-ac_sub(ac) = no;
-ac_sub(ac) = yes$(ord(ac) > (m_yeardiff_forestry(t)/5));
-
 *Reduction of ac_est is not possible.
 v32_hvarea_forestry.fx(j,ac_est) = 0;
 v32_land_reduction.fx(j,type32,ac_est) = 0;

modules/44_biodiversity/bii_target/sets.gms

@@ -8,10 +8,10 @@
 sets
 
   landcover44 land cover classes used in bii calculation
-    / crop_ann, crop_per, manpast, rangeland, urban, aff_ndc, aff_co2p, primforest, secdforest, other, plant /
+    / crop_ann, crop_per, crop_tree, manpast, rangeland, urban, aff_ndc, aff_co2p, primforest, secdforest, other, plant /
 
   bii_class44 bii coefficent land cover classes
-    / crop_ann, crop_per, manpast, rangeland, urban, primary, secd_mature, secd_young, timber /
+    / crop_ann, crop_per, crop_tree, 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 /

modules/44_biodiversity/bv_btc_mar21/sets.gms

@@ -8,10 +8,10 @@
 sets
 
   landcover44 land cover classes used in bii calculation
-    / crop_ann, crop_per, manpast, rangeland, urban, aff_ndc, aff_co2p, primforest, secdforest, other, plant /
+    / crop_ann, crop_per, crop_tree, manpast, rangeland, urban, aff_ndc, aff_co2p, primforest, secdforest, other, plant /
 
   bii_class44 bii coefficent land cover classes
-    / crop_ann, crop_per, manpast, rangeland, urban, primary, secd_mature, secd_young, timber /
+    / crop_ann, crop_per, crop_tree, 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 /