IndEcol / RECC-ODYM

The RECC model
MIT License
21 stars 10 forks source link

Forestry model "GrowthCurve" #76

Open LolaRousseau opened 5 months ago

LolaRousseau commented 5 months ago

Hi @stefanpauliuk!

I would like to use the forestry model "GrowthCurve" and it gives me some strange results for RECC_System.StockDict['dS_1t'].Values[:,:,Carbon_loc]. There is some sort of a jump in the middle of the curve (see graph below) which after some investigation of the variables to get this one indicates that there is something happening when Forest_GrowthTable_timber is being diagonally filled (and there is some replication/periodicity happening because of the sizes time/cohort) and then np.diff(Forest_GrowthTable_timber.sum(axis=1))[SwitchTime-1::] is being used.

An idea we have been discussing with @CarrerF is to model the timber similarly to the fuel wood, but this would probably not capture that there is already forest growth happening from the timber in the buildings from past cohorts.

To note: I am not running the current version of the code but it seems like the part related to the "GrowthCurve" forestry model is the same. In addition, I have modified the code to be able to run it for any chosen year of start and year of end (in my case 2022 instead of 2015 and 2050 instead of 2060 - this is why the graph is only 28 years of timeframe).

Graph of RECC_System.StockDict['dS_1t'].Values[:,:,Carbon_loc]: image

LolaRousseau commented 5 months ago

Hi @stefanpauliuk! Here are some additional clarifications to my question after running the last version of the code on github with the last version of the database from Dropbox. At lines 2601-2602: Forest_GrowthTable_timber = np.zeros((Nc,Nc)) np.fill_diagonal(Forest_GrowthTable_timber, -1* SysVar_RoundwoodConstruc_c_1_2_r[:,mmr,mS,mR]) image

This, combined with the growth curve multiplier and then differenced (with np.diff) is used to fill in RECC_System.StockDict['dS_1t'].Values[1::,mmr,Carbon_loc] (line 2620) and gives this shape: image

This is causing the jump as shown in this figure: image

Why do we use SysVar_RoundwoodConstruc_c_1_2_r?

  1. It seems we are using future trends for wood demand to establish the past wood use.
  2. The end of the period causes spikes when using np.diff.

Modeling the timber similarly to the fire wood would give a more conservative, but continuous, forest stock dynamic. Another option could be to use the historic inflows of wood (Var_I for buildings) to determine the "real" historic wood consumption.

stefanpauliuk commented 5 months ago

Salut @LolaRousseau ! I will try to look into this tomorrow! Best regards, Stefan

stefanpauliuk commented 5 months ago

Salut @LolaRousseau to be honest, I have implemented but never worked with the ForestGrowthCurve model. I agree with you that the results for timer (material wood) are not plausible.

My intention was as follows: the GrowthCurve option is an at scale consequential approach, it accounts for the removal of the current harvest (F_1_2 timber and fuelwood) and for the subsequent regrowth after harvest. This removal and growth are recorded in Forest_GrowthTable_timber and Forest_GrowthTable_fuelwd.

These Growthtable are then set as the development of the stock S1 over time. Here, I agree with you (if I understood you correctly) that the assignment of harvest and regrowth for timber should be the same as for fuel wood, the old assignment was probably a mistake.

Also, the calculation of the stock change dS1 can be done directly from S1, which will remove the strange jumps that you saw in the plots above.

I have made the following modification on my version:

                # Assign growth table values to stock and stock changes in process 1:
                # only the forest carbon pool change relative to the baseline is quantified, not the total forest carbon stock.    
                RECC_System.StockDict['S_1t'].Values[:,:,mmr,Carbon_loc]              = Forest_GrowthTable_timber[0:Nt,:]
                RECC_System.StockDict['S_1f'].Values[:,SwitchTime-1::,mmr,Carbon_loc] = Forest_GrowthTable_fuelwd 

                RECC_System.StockDict['dS_1t'].Values[:,mmr,Carbon_loc]   = RECC_System.StockDict['S_1t'].Values[0,:,mmr,Carbon_loc].sum()
                RECC_System.StockDict['dS_1t'].Values[1::,mmr,Carbon_loc] = np.diff(RECC_System.StockDict['S_1t'].Values[:,:,mmr,Carbon_loc],axis=0).sum(axis=1)
                RECC_System.StockDict['dS_1f'].Values[:,mmr,Carbon_loc]   = RECC_System.StockDict['S_1f'].Values[0,:,mmr,Carbon_loc].sum()
                RECC_System.StockDict['dS_1f'].Values[1::,mmr,Carbon_loc] = np.diff(RECC_System.StockDict['S_1f'].Values[:,:,mmr,Carbon_loc],axis=0).sum(axis=1)

        if ScriptConfig['ForestryModel'] == 'CarbonNeutral': # Default option
...

With that, and a random trial run for a wood-intensive scneario for China, I get the following result: RECC_ForestGrowthCurveExample

The forest regrowth to replace the harvest takes time, which is why there is a delay in sequestration, as can be nicely seen in the plot. This is the consequential assessment of what happens when you harvest and regrow at the same place later.

In our current assessment, we look at the forest at the landscape level as a large industry and assume carbon neutrality at the landscape level. Actual forest carbon dynamics are then assessed with real forest growth models, like CRAFT: https://onlinelibrary.wiley.com/doi/full/10.1111/gcb.15004

happy to discuss further! Since the code above does what I intended to do with the GrowthCurve model, I suggest pushing this version to the main branch. what do you think Lola?

LolaRousseau commented 5 months ago

Hallo @stefanpauliuk! I was wondering if the past timber harvest (for buildings built before 2015) should have been accounted in some way but we are not looking at the embodied emissions of these buildings, so I do not know if the on-going sequestration due to past harvest is something you wanted to account for.

Modeling timber as fuel wood also makes sense, in this way it considers only harvest from 2015 and so on. Do you still have Forest_GrowthTable_timber as size Nc? I guess it could be of size Nt as it is for Forest_GrowthTable_fuelwd and the subsequent lines of code with Nc woulc be changed by Nt.

I'll incorporate these changes in my code.

stefanpauliuk commented 5 months ago

Hi @LolaRousseau very good point! The age-cohort structure of the existing stock allows us to calculate how much timber was harvested historically to build up the current stock. In 2015, there were around 8 Gt of Carbon in timber in the standing building stock. The current implementation though is for future harvest wood only.

Is this something you could add to the model? Take the sum over all building types in 2015 (or 2020) to get an estimate for final consumption (over c, starting in 1900) of timber into the standing stock. Then use 4_PY_TimberRoundWood_V1.1 to convert from structural timber to roundwood and then have the Forest growth table start in 1900?

I honestly speaking don't remember why the Timber Growth Table has size NcNc and the Fuelwood has NtNt. Maybe because I originally intended to do what you describe above, to add the regrowth due to historic harvest for timber.

For now, I'll commit and push the correction of the GrowthCurve impacts of future harvest, as quoted above, so that working part is in the official record. Would be cool if you we managed to add the regrowth due to historic harvest, too!

LolaRousseau commented 5 months ago

Hi @stefanpauliuk , I will look into that. I think it is possible to use Var_I (for buildings) which gives the inflows of materials. I keep you informed :)