Source code for landlab.components.vegetation_dynamics.vegetation_dynamics

import numpy as np

from landlab import Component

_VALID_METHODS = {"Grid"}


[docs]def assert_method_is_valid(method): if method not in _VALID_METHODS: raise ValueError("%s: Invalid method name" % method)
[docs]class Vegetation(Component): """Landlab component that simulates net primary productivity, biomass and leaf area index at each cell based on inputs of root-zone average soil moisture. Ref: Zhou, X., Istanbulluoglu, E., & Vivoni, E. R. (2013). Modeling the ecohydrological role of aspect controlled radiation on tree grass shrub coexistence in a semiarid climate. Water Resources Research, 49(5), 2872-2895. .. codeauthor:: Sai Nudurupati and Erkan Istanbulluoglu Examples -------- >>> from landlab import RasterModelGrid >>> from landlab.components import Vegetation Create a grid on which to simulate vegetation dynamics. >>> grid = RasterModelGrid((5,4), xy_spacing=(0.2, 0.2)) The grid will need some input data. To check the names of the fields that provide the input to this component, use the *input_var_names* class property. >>> sorted(Vegetation.input_var_names) # doctest: +NORMALIZE_WHITESPACE ['surface__evapotranspiration', 'surface__potential_evapotranspiration_30day_mean', 'surface__potential_evapotranspiration_rate', 'vegetation__plant_functional_type', 'vegetation__water_stress'] >>> sorted(Vegetation.units) # doctest: +NORMALIZE_WHITESPACE [('surface__evapotranspiration', 'mm'), ('surface__potential_evapotranspiration_30day_mean', 'mm'), ('surface__potential_evapotranspiration_rate', 'mm'), ('vegetation__cover_fraction', 'None'), ('vegetation__dead_biomass', 'g m^-2 d^-1'), ('vegetation__dead_leaf_area_index', 'None'), ('vegetation__live_biomass', 'g m^-2 d^-1'), ('vegetation__live_leaf_area_index', 'None'), ('vegetation__plant_functional_type', 'None'), ('vegetation__water_stress', 'None')] Provide all the input fields. >>> grid['cell']['vegetation__plant_functional_type']= ( ... np.zeros(grid.number_of_cells, dtype=int)) >>> grid['cell']['surface__evapotranspiration'] = ( ... 0.2 * np.ones(grid.number_of_cells)) >>> grid['cell']['surface__potential_evapotranspiration_rate']= np.array([ ... 0.25547770, 0.25547770, 0.22110221, ... 0.22110221, 0.24813062, 0.24813062]) >>> grid['cell']['surface__potential_evapotranspiration_30day_mean']= ( ... np.array([0.25547770, 0.25547770, 0.22110221, ... 0.22110221, 0.24813062, 0.24813062])) >>> grid['cell']['vegetation__water_stress'] = ( ... 0.01 * np.ones(grid.number_of_cells)) Instantiate the 'Vegetation' component. >>> Veg = Vegetation(grid) >>> Veg.grid.number_of_cell_rows 3 >>> Veg.grid.number_of_cell_columns 2 >>> Veg.grid is grid True >>> import numpy as np >>> sorted(Vegetation.output_var_names) # doctest: +NORMALIZE_WHITESPACE ['vegetation__cover_fraction', 'vegetation__dead_biomass', 'vegetation__dead_leaf_area_index', 'vegetation__live_biomass', 'vegetation__live_leaf_area_index'] >>> np.all(grid.at_cell['vegetation__live_leaf_area_index'] == 0.) True >>> Veg.update() >>> np.all(grid.at_cell['vegetation__live_leaf_area_index'] == 0.) False References ---------- **Required Software Citation(s) Specific to this Component** None Listed **Additional References** Zhou, X., Istanbulluoglu, E., and Vivoni, E. R.: Modeling the ecohydrological role of aspect-controlled radiation on tree-grass-shrub coexistence in a semiarid climate, Water Resour. Res., 49, 2872– 2895, doi:10.1002/wrcr.20259, 2013. """ _name = "Vegetation" _unit_agnostic = False _info = { "surface__evapotranspiration": { "dtype": float, "intent": "in", "optional": False, "units": "mm", "mapping": "cell", "doc": "actual sum of evaporation and plant transpiration", }, "surface__potential_evapotranspiration_30day_mean": { "dtype": float, "intent": "in", "optional": False, "units": "mm", "mapping": "cell", "doc": "30 day mean of surface__potential_evapotranspiration", }, "surface__potential_evapotranspiration_rate": { "dtype": float, "intent": "in", "optional": False, "units": "mm", "mapping": "cell", "doc": "potential sum of evaporation and potential transpiration", }, "vegetation__cover_fraction": { "dtype": float, "intent": "out", "optional": False, "units": "None", "mapping": "cell", "doc": "fraction of land covered by vegetation", }, "vegetation__dead_biomass": { "dtype": float, "intent": "out", "optional": False, "units": "g m^-2 d^-1", "mapping": "cell", "doc": ( "weight of dead organic mass per unit area - measured in terms " "of dry matter" ), }, "vegetation__dead_leaf_area_index": { "dtype": float, "intent": "out", "optional": False, "units": "None", "mapping": "cell", "doc": "one-sided dead leaf area per unit ground surface area", }, "vegetation__live_biomass": { "dtype": float, "intent": "out", "optional": False, "units": "g m^-2 d^-1", "mapping": "cell", "doc": ( "weight of green organic mass per unit area - measured in terms " "of dry matter" ), }, "vegetation__live_leaf_area_index": { "dtype": float, "intent": "out", "optional": False, "units": "None", "mapping": "cell", "doc": "one-sided green leaf area per unit ground surface area", }, "vegetation__plant_functional_type": { "dtype": int, "intent": "in", "optional": False, "units": "None", "mapping": "cell", "doc": ( "classification of plants (int), grass=0, shrub=1, tree=2, bare=3, " "shrub_seedling=4, tree_seedling=5" ), }, "vegetation__water_stress": { "dtype": float, "intent": "in", "optional": False, "units": "None", "mapping": "cell", "doc": "parameter that represents nonlinear effects of water deficit on plants", }, }
[docs] def __init__( self, grid, Blive_init=102.0, Bdead_init=450.0, ETthreshold_up=3.8, ETthreshold_down=6.8, Tdmax=10.0, w=0.55, WUE_grass=0.01, LAI_max_grass=2.0, cb_grass=0.0047, cd_grass=0.009, ksg_grass=0.012, kdd_grass=0.013, kws_grass=0.02, WUE_shrub=0.0025, LAI_max_shrub=2.0, cb_shrub=0.004, cd_shrub=0.01, ksg_shrub=0.002, kdd_shrub=0.013, kws_shrub=0.02, WUE_tree=0.0045, LAI_max_tree=4.0, cb_tree=0.004, cd_tree=0.01, ksg_tree=0.002, kdd_tree=0.013, kws_tree=0.01, WUE_bare=0.01, LAI_max_bare=0.01, cb_bare=0.0047, cd_bare=0.009, ksg_bare=0.012, kdd_bare=0.013, kws_bare=0.02, method="Grid", PETthreshold_switch=0, Tb=24.0, Tr=0.01, ): """ Parameters ---------- grid: RasterModelGrid A grid. Blive_init: float, optional Initial value for vegetation__live_biomass. Converted to field. Bdead_init: float, optional Initial value for vegetation__dead_biomass. Coverted to field. ETthreshold_up: float, optional Potential Evapotranspiration (PET) threshold for growing season (mm/d). ETthreshold_down: float, optional PET threshold for dormant season (mm/d). Tdmax: float, optional Constant for dead biomass loss adjustment (mm/d). w: float, optional Conversion factor of CO2 to dry biomass (Kg DM/Kg CO2). WUE: float, optional Water Use Efficiency - ratio of water used in plant water lost by the plant through transpiration (KgCO2Kg-1H2O). LAI_max: float, optional Maximum leaf area index (m2/m2). cb: float, optional Specific leaf area for green/live biomass (m2 leaf g-1 DM). cd: float, optional Specific leaf area for dead biomass (m2 leaf g-1 DM). ksg: float, optional Senescence coefficient of green/live biomass (d-1). kdd: float, optional Decay coefficient of aboveground dead biomass (d-1). kws: float, optional Maximum drought induced foliage loss rate (d-1). method: str Method name. Tr: float, optional Storm duration (hours). Tb: float, optional Inter-storm duration (hours). PETthreshold_switch: int, optional Flag to indiate the PET threshold. This controls whether the threshold is for growth (1) or dormancy (any other value). """ super().__init__(grid) self.Tb = Tb self.Tr = Tr self.PETthreshold_switch = PETthreshold_switch self._method = method assert_method_is_valid(self._method) self.initialize( Blive_init=Blive_init, Bdead_init=Bdead_init, ETthreshold_up=ETthreshold_up, ETthreshold_down=ETthreshold_down, Tdmax=Tdmax, w=w, WUE_grass=WUE_grass, LAI_max_grass=LAI_max_grass, cb_grass=cb_grass, cd_grass=cd_grass, ksg_grass=ksg_grass, kdd_grass=kdd_grass, kws_grass=kws_grass, WUE_shrub=WUE_shrub, LAI_max_shrub=LAI_max_shrub, cb_shrub=cb_shrub, cd_shrub=cd_shrub, ksg_shrub=ksg_shrub, kdd_shrub=kdd_shrub, kws_shrub=kws_shrub, WUE_tree=WUE_tree, LAI_max_tree=LAI_max_tree, cb_tree=cb_tree, cd_tree=cd_tree, ksg_tree=ksg_tree, kdd_tree=kdd_tree, kws_tree=kws_tree, WUE_bare=WUE_bare, LAI_max_bare=LAI_max_bare, cb_bare=cb_bare, cd_bare=cd_bare, ksg_bare=ksg_bare, kdd_bare=kdd_bare, kws_bare=kws_bare, ) self.initialize_output_fields() self._cell_values = self._grid["cell"] self._Blive_ini = self._Blive_init * np.ones(self._grid.number_of_cells) self._Bdead_ini = self._Bdead_init * np.ones(self._grid.number_of_cells)
@property def Tb(self): """Storm duration (hours).""" return self._Tb @Tb.setter def Tb(self, Tb): if Tb < 0.0: raise ValueError("Tb must be non-negative") self._Tb = Tb @property def Tr(self): """Inter-storm duration (hours).""" return self._Tr @Tr.setter def Tr(self, Tr): assert Tr >= 0 self._Tr = Tr @property def PETthreshold_switch(self): """Flag to indiate the PET threshold. This controls whether the threshold is for growth (1) or dormancy (any other value). """ return self._PETthreshold_switch @PETthreshold_switch.setter def PETthreshold_switch(self, PETthreshold_switch): self._PETthreshold_switch = PETthreshold_switch
[docs] def initialize( self, Blive_init=102.0, Bdead_init=450.0, ETthreshold_up=3.8, ETthreshold_down=6.8, Tdmax=10.0, w=0.55, WUE_grass=0.01, LAI_max_grass=2.0, cb_grass=0.0047, cd_grass=0.009, ksg_grass=0.012, kdd_grass=0.013, kws_grass=0.02, WUE_shrub=0.0025, LAI_max_shrub=2.0, cb_shrub=0.004, cd_shrub=0.01, ksg_shrub=0.002, kdd_shrub=0.013, kws_shrub=0.02, WUE_tree=0.0045, LAI_max_tree=4.0, cb_tree=0.004, cd_tree=0.01, ksg_tree=0.002, kdd_tree=0.013, kws_tree=0.01, WUE_bare=0.01, LAI_max_bare=0.01, cb_bare=0.0047, cd_bare=0.009, ksg_bare=0.012, kdd_bare=0.013, kws_bare=0.02, ): # GRASS = 0; SHRUB = 1; TREE = 2; BARE = 3; # SHRUBSEEDLING = 4; TREESEEDLING = 5 """ Parameters ---------- grid: RasterModelGrid A grid. Blive_init: float, optional Initial value for vegetation__live_biomass. Converted to field. Bdead_init: float, optional Initial value for vegetation__dead_biomass. Coverted to field. ETthreshold_up: float, optional Potential Evapotranspiration (PET) threshold for growing season (mm/d). ETthreshold_down: float, optional PET threshold for dormant season (mm/d). Tdmax: float, optional Constant for dead biomass loss adjustment (mm/d). w: float, optional Conversion factor of CO2 to dry biomass (Kg DM/Kg CO2). WUE: float, optional Water Use Efficiency - ratio of water used in plant water lost by the plant through transpiration (KgCO2Kg-1H2O). LAI_max: float, optional Maximum leaf area index (m2/m2). cb: float, optional Specific leaf area for green/live biomass (m2 leaf g-1 DM). cd: float, optional Specific leaf area for dead biomass (m2 leaf g-1 DM). ksg: float, optional Senescence coefficient of green/live biomass (d-1). kdd: float, optional Decay coefficient of aboveground dead biomass (d-1). kws: float, optional Maximum drought induced foliage loss rate (d-1). """ self._vegtype = self._grid["cell"]["vegetation__plant_functional_type"] self._WUE = np.choose( self._vegtype, [WUE_grass, WUE_shrub, WUE_tree, WUE_bare, WUE_shrub, WUE_tree], ) # Water Use Efficiency KgCO2kg-1H2O self._LAI_max = np.choose( self._vegtype, [ LAI_max_grass, LAI_max_shrub, LAI_max_tree, LAI_max_bare, LAI_max_shrub, LAI_max_tree, ], ) # Maximum leaf area index (m2/m2) self._cb = np.choose( self._vegtype, [cb_grass, cb_shrub, cb_tree, cb_bare, cb_shrub, cb_tree] ) # Specific leaf area for green/live biomass (m2 leaf g-1 DM) self._cd = np.choose( self._vegtype, [cd_grass, cd_shrub, cd_tree, cd_bare, cd_shrub, cd_tree] ) # Specific leaf area for dead biomass (m2 leaf g-1 DM) self._ksg = np.choose( self._vegtype, [ksg_grass, ksg_shrub, ksg_tree, ksg_bare, ksg_shrub, ksg_tree], ) # Senescence coefficient of green/live biomass (d-1) self._kdd = np.choose( self._vegtype, [kdd_grass, kdd_shrub, kdd_tree, kdd_bare, kdd_shrub, kdd_tree], ) # Decay coefficient of aboveground dead biomass (d-1) self._kws = np.choose( self._vegtype, [kws_grass, kws_shrub, kws_tree, kws_bare, kws_shrub, kws_tree], ) # Maximum drought induced foliage loss rates (d-1) self._Blive_init = Blive_init self._Bdead_init = Bdead_init self._ETthresholdup = ETthreshold_up # Growth threshold (mm/d) self._ETthresholddown = ETthreshold_down # Dormancy threshold (mm/d) self._Tdmax = Tdmax # Constant for dead biomass loss adjustment self._w = w # Conversion factor of CO2 to dry biomass self._Blive_ini = self._Blive_init * np.ones(self._grid.number_of_cells) self._Bdead_ini = self._Bdead_init * np.ones(self._grid.number_of_cells)
[docs] def update(self): """Update fields with current loading conditions. This method looks to the properties ``PETthreshold_switch``, ``Tb``, and ``Tr`` and uses their values to calculate the new field values. """ PETthreshold_ = self._PETthreshold_switch Tb = self._Tb Tr = self._Tr PET = self._cell_values["surface__potential_evapotranspiration_rate"] PET30_ = self._cell_values["surface__potential_evapotranspiration_30day_mean"] ActualET = self._cell_values["surface__evapotranspiration"] Water_stress = self._cell_values["vegetation__water_stress"] self._LAIlive = self._cell_values["vegetation__live_leaf_area_index"] self._LAIdead = self._cell_values["vegetation__dead_leaf_area_index"] self._Blive = self._cell_values["vegetation__live_biomass"] self._Bdead = self._cell_values["vegetation__dead_biomass"] self._VegCov = self._cell_values["vegetation__cover_fraction"] if PETthreshold_ == 1: PETthreshold = self._ETthresholdup else: PETthreshold = self._ETthresholddown for cell in range(0, self._grid.number_of_cells): WUE = self._WUE[cell] LAImax = self._LAI_max[cell] cb = self._cb[cell] cd = self._cd[cell] ksg = self._ksg[cell] kdd = self._kdd[cell] kws = self._kws[cell] # ETdmax = self._ETdmax[cell] LAIlive = min(cb * self._Blive_ini[cell], LAImax) LAIdead = min(cd * self._Bdead_ini[cell], (LAImax - LAIlive)) NPP = max((ActualET[cell] / (Tb + Tr)) * WUE * 24.0 * self._w * 1000, 0.001) if self._vegtype[cell] == 0: if PET30_[cell] > PETthreshold: # Growing Season Bmax = (LAImax - LAIdead) / cb Yconst = 1 / ( (1 / Bmax) + (((kws * Water_stress[cell]) + ksg) / NPP) ) Blive = (self._Blive_ini[cell] - Yconst) * np.exp( -(NPP / Yconst) * ((Tb + Tr) / 24.0) ) + Yconst Bdead = ( self._Bdead_ini[cell] + (Blive - max(Blive * np.exp(-1 * ksg * Tb / 24.0), 0.00001)) ) * np.exp(-1 * kdd * min(PET[cell] / self._Tdmax, 1.0) * Tb / 24.0) else: # Senescense Blive = max( self._Blive_ini[cell] * np.exp((-2) * ksg * Tb / 24.0), 1 ) Bdead = max( ( self._Bdead_ini[cell] + ( self._Blive_ini[cell] - ( max( self._Blive_ini[cell] * np.exp((-2) * ksg * Tb / 24.0), 0.000001, ) ) ) * np.exp( (-1) * kdd * min(PET[cell] / self._Tdmax, 1.0) * Tb / 24.0 ), 0.0, ) ) elif self._vegtype[cell] == 3: Blive = 0.0 Bdead = 0.0 else: Bmax = LAImax / cb Yconst = 1.0 / ( (1.0 / Bmax) + (((kws * Water_stress[cell]) + ksg) / NPP) ) Blive = (self._Blive_ini[cell] - Yconst) * np.exp( -(NPP / Yconst) * ((Tb + Tr) / 24.0) ) + Yconst Bdead = ( self._Bdead_ini[cell] + (Blive - max(Blive * np.exp(-ksg * Tb / 24.0), 0.00001)) ) * np.exp(-kdd * min(PET[cell] / self._Tdmax, 1.0) * Tb / 24.0) LAIlive = min(cb * (Blive + self._Blive_ini[cell]) / 2.0, LAImax) LAIdead = min( cd * (Bdead + self._Bdead_ini[cell]) / 2.0, (LAImax - LAIlive) ) if self._vegtype[cell] == 0: Vt = 1.0 - np.exp(-0.75 * (LAIlive + LAIdead)) else: # Vt = 1 - np.exp(-0.75 * LAIlive) Vt = 1.0 self._LAIlive[cell] = LAIlive self._LAIdead[cell] = LAIdead self._VegCov[cell] = Vt self._Blive[cell] = Blive self._Bdead[cell] = Bdead self._Blive_ini = self._Blive self._Bdead_ini = self._Bdead