Not able to implement clear cut from a restart simulation

[sdeherto]

Dear ED developers
I'm trying to implement some regrowth simulations after a clear cut disturbance event. I'm doing this through the land use change files as this would be most comparable to the cases I try to compare to. I setup a new simulation (restart from previous runs history files) and provide a luc file similar to:

WEST.LONGITUDE = 26
EAST.LONGITUDE = 26.5
SOUTH.LATITUDE = 0.5
NORTH.LATITUDE = 1
BLOCK.AREA = 3066351986.63614
FIRST.LUYEAR = 1400
LAST.LUYEAR = 2300
N.PFT.HARVEST = 3
HARVEST.PFT = 2 3 4
MINDBH.1ARY = 0. 0. 0.
HARVPROB.1ARY = 1. 1. 1.
MINDBH.2ARY = 0. 0. 0.
HARVPROB.2ARY = 1. 1. 1.

year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
1993 0 0 0 10000000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

So What I assume should happen here is that during the first year of my simulation (which is 1993) all my vegetation should become pasture, however when I explore the ABG of my simulation this doesn't happen in this single year but is spread out over a longer period (see below).

Any idea why this is happening and how I could approach this to have an actual clear cut of my forest over the span of 1 year?

Many thanks in advance for any input!
Best wishes
Steven

[mpaiao]

Hi @sdeherto, this is unexpected indeed, with the transition rate you imposed in 1993 you should have almost all the area assigned to pastures. I see you are using some of the old R scripts I developed, I think they should produce a plot with the relative area by disturbance type, how does that plot look like? If that plot is not available, let me know, I may still have a version that does this somewhere.

A couple of additional thoughts considering what you are planning to do.

1. Your land use file has 1994 multiple times, I think it should be one line per year.
2. If your goal is to simulate clear-cut followed by forest regrowth, I would set the disturbance rate in column vs to a very high number, set columns sbh and vbh to -1.0 (a special flag telling to use selective logging parameters). Another alternative would be to have a very high disturbance like you have followed by a very large ps "disturbance" (abandonment) in 1994.
3. I don't know which version of ED2 you are using, but note that the land use file and header changed in one of the most recent commits. I noted that the documentation was not up to date, I just edited the Wiki page to describe the additional columns.

[sdeherto]

Hi @mpaiao, thanks for the quick response! I don't think I have that plot in my figures so I would definitely be interested to check out this other version of the plots you mentioned. I tried to implement your recommendations, however it is still not behaving as I would expect. Here are some point by point responses:

1. Thanks for spotting this mistake I fixed it!
2. I implemented both options (as shown below), however it didn't change the outcome
   year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
   1993 0 0 0 0 0 0 0 0 0 0 10000000000 -1 0 -1 0 0 0 0 0
   1994 0 0 0 0 0 0 0 0 0 10000000000 0 0 0 0 0 0 0 0 0
   1995 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   1996 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   1997 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3. Thanks for pointing this out! I updated the files following these descriptions!

[mpaiao]

If you are using the (near) latest version of ED2, the script ED2/ED/Template/Template/plot_monthly.r has the code for generating the plots by disturbance type. I think the script you are using is based on that one, so you may be able to adapt it.

ED2/ED/Template/Template/plot_monthly.r

Lines 1925 to 2195 in 396cd57

	#---------------------------------------------------------------------------------------#
	# Time series by LU. #
	#---------------------------------------------------------------------------------------#
	for (v in sequence(ntslu)){
	thistslu = tslu[[v]]
	vnam = thistslu$vnam
	description = thistslu$desc
	unit = thistslu$unit
	plog = thistslu$plog
	plotit = thistslu$plt
	#----- Check whether the user wants to have this variable plotted. ------------------#
	if (plotit && any(sellu)){
	#---------------------------------------------------------------------------------#
	# Check whether the time series directory exists. If not, create it. #
	#---------------------------------------------------------------------------------#
	outdir = file.path(outpref,"tslu")
	if (! file.exists(outdir)) dir.create(outdir)
	cat0(" + ",description," time series for all LUs.")
	#----- Load variable -------------------------------------------------------------#
	if (vnam %in% names(lu)){
	thisvar = lu[[vnam]]
	if (plog){
	xylog = "y"
	badlog = ! (thisvar %gt% 0)
	thisvar[badlog] = NA
	}else{
	xylog = ""
	}#end if
	}else{
	thisvar = array(NA,dim=c(ntimes,nlu+1))
	}#end if
	#---------------------------------------------------------------------------------#
	#----- Loop over output formats. -------------------------------------------------#
	for (o in sequence(nout)){
	fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
	dummy = open.plot( fichier = fichier
	, outform = outform[o]
	, size = eysize
	, ptsz = ptsz
	, depth = depth
	)#end open.plot
	#------------------------------------------------------------------------------#
	# Find the limit, make some room for the legend, and in case the field is #
	# a constant, nudge the limits so the plot command will not complain. #
	#------------------------------------------------------------------------------#
	xlimit = pretty.xylim(u=as.numeric(datum$tomonth),fracexp=0.0,is.log=FALSE)
	ylimit = pretty.xylim(u=thisvar[,sellu] ,fracexp=0.0,is.log=plog )
	if (plog){
	xylog = "y"
	ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
	, exp(sqrt(ylimit[2]^3/ylimit[1]))
	)#end c
	}else{
	xylog = ""
	ydrought = c(ylimit[1] - 0.5 * diff(ylimit), ylimit[2] + 0.5 * diff(ylimit))
	}#end if
	#------------------------------------------------------------------------------#
	#----- Plot settings. ---------------------------------------------------------#
	letitre = paste(description,lieu,sep=" - ")
	ley = desc.unit(desc=description,unit=unit)
	cols = lucols[sellu]
	legs = lunames[sellu]
	#------------------------------------------------------------------------------#
	#------------------------------------------------------------------------------#
	# Split the plot into two windows. #
	#------------------------------------------------------------------------------#
	par(par.user)
	layout(mat=rbind(2,1),heights=c(1.-f.leg,f.leg))
	#------------------------------------------------------------------------------#
	#------------------------------------------------------------------------------#
	# First plot: legend. #
	#------------------------------------------------------------------------------#
	par(mar=c(0.1,4.6,0.1,2.1))
	plot.new()
	plot.window(xlim=c(0,1),ylim=c(0,1))
	legend( x = "bottom"
	, inset = 0.0
	, legend = legs
	, col = cols
	, lwd = lwidth
	, ncol = min(3,pretty.box(n.sellu)$ncol)
	, title = expression(bold("Land use type"))
	, xpd = TRUE
	, bty = "n"
	)#end legend
	#------------------------------------------------------------------------------#
	#------------------------------------------------------------------------------#
	# Main plot. #
	#------------------------------------------------------------------------------#
	par(mar=c(4.1,4.6,4.1,2.1))
	plot.new()
	plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
	axis(side=1,at=whenplot8$levels,labels=whenplot8$labels,padj=whenplot8$padj)
	axis(side=2,las=1)
	box()
	title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7,log=xylog)
	if (drought.mark){
	for (n in sequence(ndrought)){
	rect(xleft = drought[[n]][1],ybottom = ydrought[1]
	,xright = drought[[n]][2],ytop = ydrought[2]
	,col = grid.colour,border=NA)
	}#end for
	}#end if
	#----- Plot grid. -------------------------------------------------------------#
	if (plotgrid){
	abline(v=whenplot8$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
	}#end if
	#----- Plot lines. ------------------------------------------------------------#
	for (n in sequence(nlu+1)){
	if (sellu[n]){
	lines(datum$tomonth,thisvar[,n],type="l",col=lucols[n],lwd=lwidth)
	}#end if
	}#end for
	#------------------------------------------------------------------------------#
	#----- Close the device. ------------------------------------------------------#
	dummy = close.plot(outform=outform[o])
	#------------------------------------------------------------------------------#
	}#end for outform
	#---------------------------------------------------------------------------------#
	}#end if (tseragbpft)
	#------------------------------------------------------------------------------------#
	}#end for tseries
	#---------------------------------------------------------------------------------------#
	#---------------------------------------------------------------------------------------#
	# Plot disturbance rate by disturbance transition. #
	#---------------------------------------------------------------------------------------#
	if (tserdist && any(seldist)){
	cat0(" + Disturbance rate time series for all disturbances.")
	for (o in sequence(nout)){
	fichier = file.path(outpref,paste0("disturb-",suffix,".",outform[o]))
	dummy = open.plot( fichier = fichier
	, outform = outform[o]
	, size = eysize
	, ptsz = ptsz
	, depth = depth
	)#end open.plot
	#---------------------------------------------------------------------------------#
	# Find the limit, make some room for the legend, and in case the field is a #
	# constant, nudge the limits so the plot command will not complain. #
	#---------------------------------------------------------------------------------#
	xlimit = pretty.xylim(u=as.numeric(datum$tomonth),fracexp=0.0,is.log=FALSE)
	ylimit = NULL
	mylucols = NULL
	mylulegs = NULL
	n = 0
	for (jlu in sequence(nlu)){
	for (ilu in sequence(nlu)){
	n = n + 1
	if (seldist[ilu,jlu]){
	ylimit = c(ylimit,lu$dist[,ilu,jlu])
	mylucols = c(mylucols,distcols [n])
	mylulegs = c(mylulegs,distnames[n])
	}#end if
	}#end for
	}#end for
	ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=FALSE)
	ydrought = c(ylimit[1] - 0.5 * diff(ylimit), ylimit[2] + 0.5 * diff(ylimit))
	#---------------------------------------------------------------------------------#
	#----- Plot settings. ------------------------------------------------------------#
	letitre = paste("Disturbance rates",lieu,sep=" - ")
	#---------------------------------------------------------------------------------#
	#---------------------------------------------------------------------------------#
	# Split the plot into two windows. #
	#---------------------------------------------------------------------------------#
	par(par.user)
	layout(mat=rbind(2,1),heights=c(1.-f.leg,f.leg))
	#---------------------------------------------------------------------------------#
	#---------------------------------------------------------------------------------#
	# First plot: legend. #
	#---------------------------------------------------------------------------------#
	par(mar=c(0.1,4.6,0.1,2.1))
	plot.new()
	plot.window(xlim=c(0,1),ylim=c(0,1))
	legend( x = "bottom"
	, inset = 0.0
	, bg = background
	, legend = mylulegs
	, col = mylucols
	, lwd = lwidth
	, ncol = min(3,pretty.box(n)$ncol)
	, title = expression(bold("Transition"))
	, xpd = TRUE
	, bty = "n"
	)#end legend
	#---------------------------------------------------------------------------------#
	#---------------------------------------------------------------------------------#
	# Main plot. #
	#---------------------------------------------------------------------------------#
	par(mar=c(4.1,4.6,4.1,2.1))
	plot.new()
	plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
	axis(side=1,at=whenplot8$levels,labels=whenplot8$labels,padj=whenplot8$padj)
	axis(side=2,las=1)
	box()
	title( main = letitre
	, xlab = "Year"
	, ylab = desc.unit(desc="Disturbance rate",unit=untab$oneoyr)
	, cex.main = 0.7
	)#end title
	if (drought.mark){
	for (n in sequence(ndrought)){
	rect(xleft = drought[[n]][1],ybottom = ydrought[1]
	,xright = drought[[n]][2],ytop = ydrought[2]
	,col = grid.colour,border=NA)
	}#end for
	}#end if
	#----- Plot grid. ----------------------------------------------------------------#
	if (plotgrid){
	abline(v=whenplot8$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
	}#end if
	#----- Plot lines. ---------------------------------------------------------------#
	n = 0
	for (jlu in sequence(nlu)){
	for (ilu in sequence(nlu)){
	n = n + 1
	if (seldist[ilu,jlu]){
	lines(datum$tomonth,lu$dist[,ilu,jlu],type="l"
	,col=distcols[n],lwd=lwidth)
	}#end if
	}#end for
	}#end for
	#---------------------------------------------------------------------------------#
	#----- Close the device. ---------------------------------------------------------#
	dummy = close.plot(outform=outform[o])
	#---------------------------------------------------------------------------------#
	}#end for outform
	#------------------------------------------------------------------------------------#
	}#end if
	#---------------------------------------------------------------------------------------#

[femeunier]

Hi @sdeherto, this is a restart run, right? Does it start in 1993? What if you restart the simulation one year before the land use change event? Does it also happen in near-bare ground simulations? Probably good to add an empty line for the first year in the land use file as well:

1992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1993 0 0 0 0 0 0 0 0 0 0 10000000000 -1 0 -1 0 0 0 0 0
1994 0 0 0 0 0 0 0 0 0 10000000000 0 0 0 0 0 0 0 0 0
1995 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1996 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[sdeherto]

Hi @mpaiao and @femeunier
Many thanks for the help so far, I've found the root of the problem, I noticed that I had forgotten to turn the IANTH_DISTURB flag to 1 in the ED2IN file, this implies that I actually don't think the .lu file was being read in during these previous simulations.
I have adapted this now and also included a grass PFT (number 16), as this seems to be required when using land use, as indicated by the PASTURE_STOCK and AGRI_STOCK flags in ED2IN.

So to recap, now I setup a restart simulation with two PFTs (3 and 16) and anthropogenic disturbance set on. However, when I try to run the model it crashes when reading in the .lu file:

Loading anthropogenic disturbance. Node: 01
[-] filelist_f: Checking prefix: /data/gent/vo/000/gvo00074/land_use/regrowth/dist_IC_Yaya_2002.lu
Prefix length: 65

• Showing all 1 files:
  [-] File #: 1 /data/gent/vo/000/gvo00074/land_use/regrowth/dist_IC_Yaya_2002.lu
  forrtl: severe (24): end-of-file during read, unit -5, file Internal List-Directed Read
  Image PC Routine Line Source
  ed_2.2-opt-master 000000000107D7BD Unknown Unknown Unknown
  ed_2.2-opt-master 000000000107CD4C Unknown Unknown Unknown
  ed_2.2-opt-master 0000000001011766 ed1_fileinfo_ 242 ed_filelist.F90
  ed_2.2-opt-master 0000000000D44274 landuse_init_mp_r 144 landuse_init.f90
  ed_2.2-opt-master 00000000006B538F ed_init_history_m 622 ed_init_history.f90
  ed_2.2-opt-master 0000000000411161 ed_driver_ 196 ed_driver.F90
  ed_2.2-opt-master 0000000000410414 MAIN__ 289 edmain.F90
  ed_2.2-opt-master 000000000040F6CD Unknown Unknown Unknown
  libc.so.6 00001517B9629590 Unknown Unknown Unknown
  libc.so.6 00001517B9629640 __libc_start_main Unknown Unknown
  ed_2.2-opt-master 000000000040F5E5 Unknown Unknown Unknown

I have tried starting my simulation earlier as @femeunier suggested or making an entry for each year in the lu file but I keep encountering this same issue. My lu file currently reads as such (example for a run starting in 2002):

WEST.LONGITUDE = -4
EAST.LONGITUDE = -3.5
SOUTH.LATITUDE = 5.5
NORTH.LATITUDE = 6
BLOCK.AREA = 3066351986.63614
FIRST.LUYEAR = 1400
LAST.LUYEAR = 2300
SKID.AREA = 0.6
N.PFT.HARVEST = 3
HARVEST.PFT = 2 3 4
MINDBH.SLOG = 0. 0. 0.
HARVPROB.SLOG = 1. 1. 1.
MINDBH.FPLT = 0. 0. 0.
HARVPROB.FPLT = 1. 1. 1.
SKID_DBH_THRESH = 1. 1. 1.
SKID_S_GTHARV = 0. 0. 0.
SKID_S_LTHARV = 0. 0. 0.
FELLING_S_LTHARV = 0. 0. 0.
year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
2002 0 0 0 0 0 0 0 0 0 0 10000000000 -1 0 -1 0 0 0 0 0
2003 0 0 0 0 0 0 0 0 0 10000000000 0 0 0 0 0 0 0 0 0
2004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2005 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2006 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

I cannot find a clear mistake in the structure of my.lu file, but I'm probably missing something again?
Thanks for all the help in nay case, I feel like we are getting close to implementing this properly now!

Cheers
Steven

[mpaiao]

@sdeherto Glad to hear that the issue was something as simple as turning on IANTH_DISTURB, although this raises the question of why the biomass was collapsing to zero in your previous runs...

In any case, the current problem seems to be an error reading in the lu files. My suggestion is to take a look at the landuse_init.f90 file of your version of ED2 and make sure it is consistent with the current file format. The header changed in a recent pull request, so it may be worth double checking that the new lines added are expected in your version.

[sdeherto]

Hi @mpaiao and @femeunier
I've solved the issue, it was a multitude of things:

• the error described in the last posts was just due to not abiding to the naming convention for the landuse file, these should always end with lat and lon in the name (e.g. prefix-lat0.5lon25.lu)
• The weird delayed collapse of my vegetation I showed in the initial plot was just caused by restarting from a history run with one PFT defined (nr 4) while having a single different PFT active in this simulation (PFT nr 3). This seems to cause ED to first let the history PFT die off before the establishment of the new PFT occurs. So in the end this had nothing to do with disturbance at all!
• By correctly implementing the land use file and doing a restart with consistent PFTs it now implements the clear cut as expected!

[sdeherto]

One more thing related to this: I was testing different ways to implement the disturbance (i.e. option 1: transition from primary to secondary land, vs or option 2: transition from primary to pasture, ps) and I noticed that while the logging was complete for the second case it didn't do a complete clear cut for the case where we convert to secondary land. In fact for option 1, only about 50% of the area was logged.

option 1
year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
2017 0 0 0 0 0 0 0 0 0 0 100000000000 -1 0 -1 0 0 0 0 0

option 2
year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
2017 0 0 0 10000000000000 0 0 0 0 0 0 0 -1 0 -1 0 0 0 0 0

After looking into this somewhat deeper with the help of @femeunier, we found that by using the -1 flag for biomass that should be harvested it was actually only removing patches that are indicated as mature by the model, thus small and young patches remained intact. We were able to force the model to remove these by also adding high values for sbh and sbh2 to the transition matrix (see below).

solution
year cp pc pv vp vc cv sc cs sp ps vs sbh f_sbh vbh f_vbh sbh2 f_sbh2 vbh2 f_vbh2
2017 0 0 0 0 0 0 0 0 0 0 100000000000 1000000000000 0 100000000000000 0 1000000000000000 0 100000000000 0

Thanks for all the help, this works now!

[mpaiao]

Great to hear that it worked, @sdeherto! Regarding your last comment on the logging only impacting mature forests, there is a parameter in ED2 that controls the minimum age of a patch to be considered mature (mature_harvest_age). This can be defined through the xml parameter file. Maybe this parameter, along with some other logging and forestry parameters, should be eventually transferred to the lu file, because logging and forest plantation settings are site- and scenario-specific.

[femeunier]

Hi @mpaiao,

I wasn't aware of that parameter and looking at the code now, settings mature_harvest_age = 0 in the config file should lead to the same results, indeed. Agreed it could be moved to the lu files ..

The only thing that is still bugging me with those simulations is that we mimic a slah and burn agriculture event with clear-cut logging, which is not the same. Am I correct that prescribing fire is currently not possible with ED2? I see that this idea was left open by @mdietze with the event file but not implemented, right?

[mdietze]

For the events file (which I haven't touched in many years) I think the question is what specific pools do you want to change, and in what ways, that would be different in a slash and burn scenario from a clear cut? Yes, the fire event is still not implemented, but the harvest event does let you separately control the fractions of above ground wood, below ground biomass, foliage, and storage removed. I think the main thing missing would be the ability to manipulate the soil pools. Also worth noting is that the events code wasn't really writen for complex age and cohort structures -- it was added at a time when my lab was focused on agriculture and coppice forestry (biofuels). But I do think the underlying data structures are flexible and fairly easily extensible