Code Monkey home page Code Monkey logo

gday's Introduction

GDAY model

GDAY (Generic Decomposition And Yield) is a simple ecosystem model that simulates carbon, nitrogen, and water dynamics at the stand scale (Comins and McMurtrie, 1993; Medlyn et al. 2000; Corbeels et al. 2005a,b).

The model can be run at either a daily time step, or a 30-minute time step. When the model is run at the sub-daily timescale, photosynthesis is calculated using a two-leaf (sunlit/shade) approximation (de Pury and Farquhar, 1997; Wang and Leuning, 1998), otherwise photosynthesis is calculated following Sands (1995;1996). The sub-daily approach (photosynthesis & leaf energy balance) mirrors MAESPA, without the complexity of the radiation treatment. In the standard model the water balance is represented simply, with two (fixed) soil water "buckets", which represent a top soil (e.g. 5 cm) and a larger root-zone. If you are using the sub-daily version, there is now the option to use a SPA-style representation of hydraulics. GDAY-SPA resolves multiple soil layers, soil and leaf water potential and limits gas exchange following the Emax approach (see MAESPA for more details).

GDAY uses a modified version of the CENTURY model to simulate soil carbon and nutrient dynamics (Parton et al. 1987; 1993).

Installation

To get the code your best route is probably to fork the repository, there is a nice explanation on github.

The model is coded entirely in C without any dependancies. The wrapper files for the example scripts and the script used to change parameter options, are written in python. The old python version is still online.

There is a Makefile in the src directory...

$ make clean ; make

The Makefile will need to be edited by hand to set the $ARCH flag, which sets the installation path. Currently it is hardwired to my computer.

Running the model

A simple model usage can be displayed by calling GDAY as follows:

$ gday -u

Presently, there are only two options which the user can set via the command line. All model options, including all of the model parameters are customisable via the parameter file.

To spin-up GDAY:

$ gday -s -p param_file.cfg

To run GDAY:

$ gday -p param_file.cfg

When the model is run it expects to find its "model state" (i.e. from a previous spin-up) in the parameter file. This state is automatically written the parameter file after the initial spin-up when the "print_options" flag has been set to "end", rather than "daily".

Parameter file

GDAY expects a parameter file to be supplied as an argument (-p filename) on the command line. Parameter files follow the standard .INI format, although only a simple INI parser has been coded into GDAY.

Parameter files are broken down into 6 section, namely [git], [files], [params], [control], [state] and [print]. The order of these sections shouldn't make any difference. The basic element contained in the parameter file is the key or property. Every key has a name and a value, delimited by an equals sign (=). The name appears to the left of the equals sign.

eac = 79430.0

As mentioned above, keys are grouped into sections:

[control]
alloc_model = allometric
deciduous_model = false

or

[files]
cfg_fname = params/NCEAS_DUKE_model_youngforest_amb.cfg
met_fname = met_data/DUKE_met_data_amb_co2.csv
out_fname = outputs/D1GDAYDUKEAMB.csv
out_param_fname = params/NCEAS_DUKE_model_simulation_amb.cfg

As all the model parameters are accessible via this file, these files can be quite long. Clearly it isn't necessary to list every parameter. The recommended approach is to use the base file and then customise whichever parameters are required via a shell script, e.g. see the python wrapper script. This file just lists the parameters which needs to be changed and calls adjust_gday_param_file.py to swap the parameters (listed as a python dictionary) with the default parameter. I have also written an equivalent version in R adjust_gday_param_file.R. I should highlight that I wouldn't necessarily trust the default values :).

Finally, the options to print different the state and flux variables on the fly is a nice hangover from the python implementation. Sadly, this functionality doesn't actually exist in the C code, instead all the state and flux variables used in the FACE intercomparisons are dumped as standard.

When I have time I will write something more extensive (ha), but information about what different variable names refer to are listed in the header file, which documents the different structures (i.e. control, state, params).

The git hash allows you to connect which version of the model code produced which version of the model output. I'd argue for maintaining this functionality, but if you don't use git or wish to ignore me, filling this line with gibberish and disabling the shell command in the Makefile should allow you to do this.

Potential gotchas

  • The parameter alpha_j which represents the quantum yield of electron transport (mol mol-1) is the intrinsic quantum yield (i.e. per unit APAR). For the two-leaf version of the model, alpha_j should be divided by (1.0 - omega), where omega is the leaf scattering coefficient of PAR (leaf reflectance and transmittance combined). Currently we are assuming omega = 0.15 (radiation.c), this is currently hardwired.

  • The deciduous phenology scheme does not currently work with the two-leaf version of the model (can be fixed).

  • Wind speed must be > 0.0, some flux files have timesteps where this isn't the case and it will crash the model. Similarly, with the Medlyn gs model, VPD must be greater than 0.05 kPa.

Meteorological driving file

30-minute file:

Variable Description Units
year
doy day of year [1-365/6]
hod hour of day [0-47]
rain rainfall mm 30 min-1
par photosynthetically active radiation umol m-2 s-1
tair air temperature deg C
tsoil soil temperature deg C
vpd vapour pressure deficit kPa
co2 CO2 concentration ppm
ndep nitrogen deposition t ha-1 30 min-1
nfix biological nitrogen fixation t ha-1 30 min-1
wind wind speed m s-1
press atmospheric pressure kPa

Day file:

Variable Description Units
year
doy day of year [0-365/6]
tair (daylight) air temperature deg C
rain rainfall mm day-1
tsoil soil temperature deg C
tam morning air temperature deg C
tpm afternoon air temperature deg C
tmin minimum (day) air temperature deg C
tmax minimum (day) air temperature deg C
tday day average air temperature (24 hrs) deg C
vpd_am morning vapour pressure deficit kPa
vpd_pm afternoon vapour pressure deficit kPa
co2 CO2 concentration ppm
ndep nitrogen deposition t ha-1 day-1
nfix biological nitrogen fixation t ha-1 day-1
wind wind speed m s-1
press atmospheric pressure kPa
wind_am morning wind speed m s-1
wind_pm afternoon wind speed m s-1
par_am morning photosynthetically active radiation MJ m-2 d-1
par_pm afternoon photosynthetically active radiation MJ m-2 d-1

Nitrogen inputs

Nitrogen (N) entering the system via biological N fixation (BNF; tonnes ha-1 yr-1) and N deposition (tonnes ha-1 yr-1) are prescribed and passed via the met file. If information isn't available from the experiment GDAY is being applied to, BNF can be calculated as a function of evapotranspiration (ET) based on Cleveland et al. 1999.

Following Smith et al. (2014), Biogeosciences and Wieder et al. (2015), ERL, we also suggest the conservation BNF equation (Fig. 1). For estimates of ET you can either use values by PFT based on Table 1 in Cleveland or use the sum of canopy evaporation and transpiration. Wieder et al (pg 3) argued that using total ET leads to a high bias in BNF estimates in arid regions.

BNF (kg N ha-1 yr-1) is then calculated as a function of ET:

BNF = 0.102 * (ET * mm_2_cm) + 0.524

Hydraulics

From SPA we borrow the multi-layer soil scheme, which considers infiltration and drainage between layers. We also implement the soil-to-leaf hydraulics from SPA, which includes weighting soil water potential. We limit gas exchange following the Emax approach (Duursma et al. 2008). This approach therefore assumes isohydric behaviour, i.e. the plant maintains a leaf water potential above a critical minimum value.

We do not currently implement the thermal calculations which would allow you to estimate soil temperature.

Model spinup

There are now two options to spin the model up: (i) brute force, i.e. continusly recycling the met data and re-running GDAY (Flag - spinup_method = brute); and (ii) an implementation of the Semi-Analytical Solution (SAS) for steady-state approximation (Flag - spinup_method = sas). The SAS approach is about 60% faster than the slow spinup so I'd suggest using it.

The SAS method is after:

  • Xia, J., Y. Luo, Y.-P. Wang, and O. Hararuk. 2013. Traceable components of terrestrial carbon storage capacity in biogeochemical models. Global Change Biology 19:2104-2116

Example run

The example directory has two python scripts which provide an example of how one might set about running the model. example.py simulates the DUKE FACE experiment and run_experiment.py is just nice a wrapper script around this which produces a plot at the end comparing the data to the observations.

cd example/
python run_experiment.py

This should pop a plot of NPP, LAI and transpiration onto your screen. This example tends to break from time to time when I change various options, so please let me know if it isn't working!

NB to use this wrapper script you will need to have an installation of the Pandas and Matplotlib libraries installed. If you are a python user this is fairly standard.

Key References

  1. Comins, H. N. and McMurtrie, R. E. (1993) Long-Term Response of Nutrient-Limited Forests to CO2 Enrichment; Equilibrium Behavior of Plant-Soil Models. Ecological Applications, 3, 666-681.
  2. Medlyn, B. E., McMurtrie, R. E., Dewar, R. C. and Jeffreys, M. P. (2000), Soil processes dominate the long-term response of forest net primary productivity to increased temperature and atmospheric CO2 concentration, Canadian Journal of Forest Research, 30, 873–888.
  3. Corbeels M, McMurtrie RE, Pepper DA, O’Connell AM (2005a) A process-based model for nitrogen cycling in forest planta- tions Part I. Structure, calibration and analysis of decomposi- tion model. Ecological Modelling, 187, 426–448.
  4. Corbeels M, McMurtrie RE, Pepper DA, O’Connell AM (2005b) A process-based model for nitrogen cycling in forest planta- tions Part II. Simulating growth and nitrogen mineralisation of Eucalyptus globulus plantations in south-western Australia. Ecological Modelling, 187, 449–474.
  5. Sands PJ (1995) Modelling canopy production. II. From single-leaf photosynthetic parameters to daily canopy photosynthesis. Australian Journal of Plant Physiology, 22, 603-614.
  6. Sands PJ (1996) Modelling canopy production. III. Canopy light-utilisation efficiency and its sensitivity to physiological environmental variables. Australian Journal of Plant Physiology, 23, 103-114.
  7. de Pury, D.G.G., Farquhar, G.D. (1997) Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models. Plant, Cell and Environment, 20, 537-557.
  8. Wang, Y-P. and Leuning R. (1998) A two-leaf model for canopy conductance, photosynthesis and portioning of available energy I: Model description and comparison with a multi-layered model. Agricultural and Forest Meteorology, 91, 89–111.
  9. Parton, W.J., Schimel, D.S., Cole, C.V., and Ojima, D.S. (1987) Analysis of factors controlling soil organic matter levels in Great Plains grasslands. Soil Sc. Soc. Am. J. 51: 1173–1179.
  10. Parton, W.J., Scurlock, J.M.O., Ojima, D.S., Gilmanov, T.G., Scholes, R.J., Schimel, D.S., Kirchner, T., Menaut, J.-C., Seastedt, T., Garcia Moya, E. Kamnalrut, A., and Kinyamario, J.I. (1993) Observations and modeling of biomass and soil or- ganic matter dynamics for the grassland biome worldwide. Global Biogeochem. Cycles, 7: 785–809.
  11. Duursma, RA and Kolari, P and Perämäki, M and Nikinmaa, E and Hari, P and Delzon, S and Loustau, D and Ilvesniemi, H and Pumpanen, J and Mäkelä, A. (2008) Predicting the decline in daily maximum transpiration rate of two pine stands during drought based on constant minimum leaf water potential and plant hydraulic conductance. Tree physiology, 28, 265-276.
  12. Duusma, R. A. and Medlyn, B. E. (2012) MAESPA: a model to study interactions between water limitation, environmental drivers and vegetation function at tree and stand levels, with an example application to [CO2] x drought interactions. Geoscientific Model Development, 5, 919-940.

Contacts

gday's People

Contributors

mdekauwe avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

gday's Issues

Respiration

Someone needs to add maintenance and growth respiration functions

S.Hem pheno

Currently doesn't work

  • Added logic to calculate leaf on/off dates
  • Added logic to calculate transfer fluxes respecting s.hem growth

Outstanding Issues:

  • in gday.c things are being zero'd at the end of the year which is of course the middle of the year. However I quickly tried moving this to DOY==182 and that didn't fix it.
  • the stress calculations are also assuming year end == 365
  • the biggy...I'm assuming the pheno can be calculated from the current year when really it should be taking in the last half ot the previous year and ditto the final growth period at year end. Obviously matters, but less than why it currently doesn't work

SPA spin-up

If you're using SPA then the time involved to spin-up the model is pretty costly.

Test case: model spin-up & then a 12 year simulation at a 30 minute timestep for Tumbarumba

First, the impact of SAS on simulation time = ~57% improvement using the SAS spin-up.
SPA off, brute force spin-up (i.e. our classic GDAY just at 30 mins): 2m 20
SPA off, SAS spin-up: 0 m 59

But turning SPA on, leads to a ridiculously long simulation time.
SPA on, SAS spin-up: 17 m 46

Causes?

  1. Drainage through the soil water profile, this involves all the costly calls to numerical recipes for the integration funcs (odein, rkck, rkqs).
  2. Calculating a new rooting distribution at the end of the year, this involves a call to zbrent to do a minimisation.
  3. Using 20 soil layers, obviously that also impacts on (1) as well.
  4. Potentially too strict a stopping criteria for the passive pool.

Solution 1

Assume a fixed drainage rate, I took an average across 10,000 SPA iterations for the Tumbarumba soil, finding a rate ~7%. So I tested the assumption that 10% of water in a layer can drain during a step for the spin-up (obviously one would test for a more robust value!). Then for the actual simulations I go back to using the SPA calculations for drainage. This gets you a total run time of:

SPA on, SAS spin-up, fast drainage: 3 m 43.6

Solution 2

Either (1) assume a fixed rooting distribution as in CABLE, or (2) fix the slope term in the rooting distribution function from SPA, to remove the call to zbrent. I tested 2 and it makes things a few seconds faster, so this is a small issue.

Solution 3

A bit arbitrary...We could play with a few different soil layers numbers and see at what point we lose important resolution.

I tested using 6 layers instead of 20:

SPA on, SAS spin-up, 6 layers: 5 m 56

Solution 4

Relax the stopping criteria. Currently we stop once the change in the passive pool is < 0.05 t/ha, if I make that 0.5 t/ha

SPA on, SAS spin up, delta passive < 0.5: 8 m 3

Moving to the impact on simulated fluxes ...

Solution 1: Drainage

Here EXP is the fixed rate drainage and CTRL is standard GDAY-SPA. NB. the fixed drainage rate is only used during spin-up, during the actual simulation we switch back.

Dry down 1

tumbarumba_spa_comparison_20020801_20030801

Dry down 2

tumbarumba_spa_comparison_20040101_20040101

Makes no differences, so I think that is a viable spin-up assumption?

Solution 2: 6 layers vs. 20

As you'd expect this has more of an impact, not sure what the right balance is here. For this test I just made 6 equal sized layers, clearly a few smaller ones and a single larger one might be preferable.

Dry down 1

tumbarumba_spa_comparison_20020801_20030801

Dry down 2

tumbarumba_spa_comparison_20040101_20040101

Solution 4: relaxed stopping criteria

Dry down 1

tumbarumba_spa_comparison_20020801_20030801

Dry down 2

tumbarumba_spa_comparison_20040101_20040101

How fast could we go?

Fixed drainage, relaxed stopping criteria: 1 m 42s
Fixed drainage, relaxed stopping criteria, 6 layers: 0 m 48s

LUE component

Target: Develop a model that has grass phenology and species competition

Issues: I can't get the LUE component to produce a reasonable value.

Not sure what I did wrong but here's something I found:

  1. The met file from the example folder did not have the columns in the same order as the model assumes to read in. This is because the model reads each column in a fix order (hardwired) without checking the column names. I change the columns in the met file for now but will try to change the code to check the column names;
  2. The f->apar in L542 in photosynthesis is timed by 60.0 * 60.0 * daylen; which assmues the apar is mumol m-2 s-1 but the apar should have the same unit asm->par which is mj m-2 d-1;
  3. I'm still don't have a solution yet but will update as it moves on.
  4. Realted to 2. PAR is defined as the sum of par_am and par_pm from met data: m->par = ma->par_am[c->day_idx] + ma->par_pm[c->day_idx]. par_ in the met data should be in mumol m-2 s-1. This makes the sum of the two not informative?
  5. I think based on the code, that the readme file listed wrong units. par_am and pm should be mj m-2 halfday-1.

Small unrealted things:

  1. Small typo in water_balace.c: "top_soil_loss" should be "topsoil_loss"?
  2. In read_metdata.c, "current_yr" did not has a initial value;

Windoes freindly changes:

  1. I silenced all the git related bits since it doesn't provide any use for now;
  2. The parameter and met data files are not fixed to be par.cfg and met.csv which should be in the same folder as the exe.

Interception

The interception code I've added borrows from CABLE and is consistent with the simple nature of the day model.

  • One could look at replacing this with the Rutter model
  • The current scheme ignores any drainage from the canopy store during each time step. There is "spill" but no drainage.
  • Need to consider whether we go down the route of using some of the wet fraction of the canopy to meet atmos. demand. Currently we are not as before.

Hydraulics

  • Calculate psi_soil and use psi_pd to set drought limitation, i.e. g1 = f(psi_pd)
  • Calculate psi_lead during stability loop, but for the moment just output it.
  • Also need to calculate psi_stem.
  • Loop at what point we hit PLC 88.

Sub-daily dynamics

  1. Add new flag "sub-daily", if called we look for 30min data (or if not found we call a weather generator?)
  2. read 30 minute data.
  3. add routine to determine direct/diffuse component of PAR
  4. add internal day loop over (48 time steps), if sun is up call new photosynthesis routine.
  5. add new photosynthesis routine, leaf energy balance, resolve Ci, VPD, etc. For the moment keep the SW as a bucket.
  6. Once the sun has gone down, skip out of sub-daily loop.
  7. Sum daily C, water.
  8. Join standard part of code and move C, N around the plant.

N cut back logic

Explore whether this works with new sunlit/shaded code:

options:

  • up-regulate rateuptake parameter temporarily, see if this works...
  • or stop leaf turnover.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.