# -*- coding: latin-1 -*-
from __future__ import division # use '//' to do integer division
from math import sqrt, log, exp
from openalea.cnwgrass.gasexchange import parameters
"""
gasexchange.model
~~~~~~~~~~~~~~~~~~~
Model of photosynthesis based on Farquhar's approach.
The model includes the dependence of photosynthesis to organ temperature and nitrogen content.
Internal CO2 and organ temperature are found numerically.
"""
# TODO: extract all parameters and put them in farqhuar.parameters
[docs]
def organ_temperature(w, z, Zh, Ur, PAR, gsw, Ta, Ts, RH, organ_name):
"""
Energy balance for the estimation of organ temperature
:param float w: organ characteristic dimension (m) to be considered for heat transfer through forced convection (by wind).
For a leaf: its width (more related to wind direction than length), for cylindric stem elements: diameter.
:param float z: organ height from soil (m)
:param float Zh: canopy height (m)
:param float Ur: wind speed (m s-1) at the reference height (zr), e.g. top of the canopy + 2m (in the case of wheat, Ur can be approximated as the wind speed at 2m from soil)
:param float PAR: absorbed PAR (µmol m-2 s-1)
:param float gsw: stomatal conductance to water vapour (mol m-2 s-1)
:param float Ta: air temperature (degree C)
:param float Ts: organ temperature (degree C). Ts = Ta at the first iteration of the numeric resolution
:param float RH: Relative humidity (decimal fraction)
:param str organ_name: name of the organ to which belongs the element (used to distinguish lamina from cylindric organs)
:return: Ts (organ temperature, degree C), Tr (organ transpiration rate, mm s-1)
:rtype: (float, float)
"""
d = parameters.Zh_d * Zh #: Zero plane displacement height (m) #TODO: should we adapt this calculation before heading ?
Zo = parameters.Zh_Zo * Zh #: Roughness length (m)
# TODO: Temporary patch to avoid div 0 error
Ur = max(Ur, parameters.Ur_min)
#: Wind speed
u_star = (Ur * parameters.K) / log((parameters.ZR - d) / Zo) #: Friction velocity (m s-1)
Uh = (u_star / parameters.K) * log((Zh - d) / Zo) #: Wind speed at the top of canopy (m s-1)
u = Uh * exp(parameters.A * (z / Zh - 1)) #: Wind speed at organ height (m s-1), from Campbell and Norman (1998), second edition.
#: Boundary layer resistance to heat (s m-1). See Finnigan J, Raupach M. 1987 and Monteith JL. 1973 for basic equations.
if organ_name == 'blade':
rbh = parameters.rhb_blade_A * sqrt(w / u) #: Case of horizontal planes submitted to forced convection
else:
rbh = w / (parameters.rhb_other_A * ((u * w) / parameters.rhb_other_B) ** parameters.rhb_other_C) #: Case of vertical cylinders submitted to forced convection
#: Turbulence resistance to heat (s m-1)
ra = 1 / (parameters.K ** parameters.ra_expo * Ur) * (log((parameters.ZR - d) / Zo)) ** parameters.ra_expo #: Aerodynamic resistance integrated from zr to z0 + d
#: Net absorbed radiation Rn (PAR and NIR, J m-2 s-1)
RGa = (PAR * parameters.PARa_to_RGa) / parameters.Watt_to_PPFD #: Global absorbed radiation by organ (J m-2 s-1).
es_Ta = parameters.s_C * exp((parameters.s_B * Ta) / (parameters.s_A + Ta)) #: Saturated vapour pressure of the air (kPa), Ta in degree Celsius
V = RH * es_Ta #: Vapour pressure of the air (kPa)
# fvap = 0.56 - 0.079*sqrt(10*V) #: Fraction of vapour pressure
# tau = RGa/parameters.I0 #: Atmospheric transmissivity (dimensionless)
# fclear = 0.1 + 0.9*max(0, min(1, (tau-0.2)/0.5)) #: Fraction sky clearness
Rn = RGa # # NB: this only accounts for the visible radiations. More general equation would be RGa + 0.95 * parameters.SIGMA * ((Ta + parameters.KELVIN_DEGREE)**4 + 0.97 * (Tsoil + parameters.KELVIN_DEGREE)**4 - 2*(Ts + parameters.KELVIN_DEGREE)**4)
# General equation is Rn = RGa + epsilon*Ra + epsilon*sigma*(Ts_feuilles_voisines + parameters.KELVIN_DEGREE)**4 - epsilon*sigma*(Ts + parameters.KELVIN_DEGREE)**4
# if Ra unavailable, use Ra = sigma*(Tair + parameters.KELVIN_DEGREE)**4*fvap*fclear
#: Transpiration (mm s-1), Penman-Monteith
if Ts == Ta:
Ta_K = Ta + parameters.KELVIN_DEGREE
s = ((parameters.s_B * parameters.s_A) / (Ta_K + parameters.s_A) ** parameters.s_expo) * es_Ta #: Slope of the curve relating saturation vapour pressure to temperature (kPa K-1)
else:
es_Tl = parameters.s_C * exp((parameters.s_B * Ts) / (parameters.s_A + Ts)) #: Saturated vapour pressure at organ level (kPa), Ts in degree Celsius
Ts_K, Ta_K = Ts + parameters.KELVIN_DEGREE, Ta + parameters.KELVIN_DEGREE
s = (es_Tl - es_Ta) / (Ts_K - Ta_K) #: Slope of the curve relating saturation vapour pressure to temperature (kPa K-1)
VPDa = es_Ta - V
rbw = parameters.rbh_rbw * rbh #: Boundary layer resistance for water (s m-1)
gsw_physic = (gsw * parameters.R * (Ts + parameters.KELVIN_DEGREE)) / parameters.PATM #: Stomatal conductance to water in physical units (m s-1). Relation given by A. Tuzet (2003)
rswp = 1 / gsw_physic #: Stomatal resistance for water (s m-1)
Tr = max(0., (s * Rn + (parameters.RHOCP * VPDa) / (rbh + ra)) / (
parameters.LAMBDA * (s + parameters.GAMMA * ((rbw + ra + rswp) / (rbh + ra))))) #: mm s-1
#: Organ temperature
Ts = Ta + ((rbh + ra) * (Rn - parameters.LAMBDA * Tr)) / parameters.RHOCP
return Ts, Tr
[docs]
def stomatal_conductance_BWB(Ag, An, surfacic_nitrogen, ambient_CO2, RH):
"""
Ball, Woodrow, and Berry model of stomatal conductance (1987)
:param float Ag: gross assimilation rate (µmol m-2 s-1)
:param float An: net assimilation rate (µmol m-2 s-1)
:param float surfacic_nitrogen: surfacic nitrogen content(g m-2) including or not structural nitrogen depending on parameter.MODEL_VERSION
:param float ambient_CO2: Air CO2 (µmol mol-1)
:param float RH: Relative humidity (decimal fraction)
:return: gsw (mol m-2 s-1)
:rtype: float
"""
Cs = ambient_CO2 - An * (parameters.K_Cs / parameters.GB) #: CO2 concentration at organ surface (µmol mol-1 or Pa). From Prieto et al. (2012). GB in mol m-2 s-1
m = parameters.PARAM_N['delta1'] * surfacic_nitrogen ** parameters.PARAM_N['delta2'] #: Scaling factor dependence to surfacic_nitrogen (dimensionless). This function is maintained
# although I'm not convinced that it should be taken into account
gsw = (parameters.GSMIN + m * ((Ag * RH) / Cs)) #: Stomatal conductance to water vapour (mol m-2 s-1), from Braune et al. (2009), Muller et al. (2005): using Ag rather than An.
return gsw
[docs]
def stomatal_conductance_Leuning(Ag, An, Ta, ambient_CO2, RH):
"""
Leuning model of stomatal conductance to water (1995)
:param float Ag: gross assimilation rate (µmol m-2 s-1)
:param float An: net assimilation rate (µmol m-2 s-1)
:param float Ta: air temperature (degree C)
:param float ambient_CO2: Air CO2 (µmol mol-1)
:param float RH: Relative humidity (decimal fraction)
:return: gsw (mol m-2 s-1)
:rtype: float
"""
es_Ta = parameters.s_C * exp((parameters.s_B * Ta) / (parameters.s_A + Ta)) #: Saturated vapour pressure of the air (kPa), Ta in degree Celsius
V = RH * es_Ta #: Vapour pressure of the air (kPa)
VPDa = es_Ta - V
fw = 1 / (1 + VPDa / parameters.D0)
Cs = ambient_CO2 - An * (parameters.K_Cs / parameters.GB) #: CO2 concentration at organ surface (µmol mol-1 or Pa). From Prieto et al. (2012). GB in mol m-2 s-1
gamma = parameters.GAMMA0 * (1 + parameters.GAMMA1 * (Ta - parameters.T_ref) + parameters.GAMMA2 * ((Ta - parameters.T_ref) * (Ta - parameters.T_ref)))
gsw = parameters.GSMIN + (1.6 * parameters.m * Ag * fw) / (Cs - gamma)
return gsw
[docs]
def stomatal_conductance_Tuzet(Ag, An, Ta, ambient_CO2, water_potential):
"""
Tuzet model of stomatal conductance to water (2003)
:param float Ag: gross assimilation rate (µmol m-2 s-1)
:param float An: net assimilation rate (µmol m-2 s-1)
:param float Ta: air temperature (degree C)
:param float ambient_CO2: Air CO2 (µmol mol-1)
:param float water_potential: water potential of the organ (Mpa)
:return: gsw (mol m-2 s-1)
:rtype: float
"""
fw = 1 / (1 + (water_potential / parameters.water_potential_ref) ** parameters.n)
Cs = ambient_CO2 - An * (parameters.K_Cs / parameters.GB) #: CO2 concentration at organ surface (µmol mol-1 or Pa). From Prieto et al. (2012). GB in mol m-2 s-1
gamma = parameters.GAMMA0 * (1 + parameters.GAMMA1 * (Ta - parameters.T_ref) + parameters.GAMMA2 * ((Ta - parameters.T_ref) * (Ta - parameters.T_ref)))
gsw = parameters.GSMIN + (1.6 * parameters.m * Ag * fw) / (Cs - gamma)
return gsw
[docs]
def stomatal_conductance_hydraulics(Ag, An, Ta, ambient_CO2, RH, water_potential):
"""
Model of stomatal conductance to water coupling Tuzet and Leuning
:param float Ag: gross assimilation rate (µmol m-2 s-1)
:param float An: net assimilation rate (µmol m-2 s-1)
:param float Ta: air temperature (degree C)
:param float ambient_CO2: Air CO2 (µmol mol-1)
:param float RH: Relative humidity (decimal fraction)
:param float water_potential: water potential of the organ (Mpa)
:return: gsw (mol m-2 s-1)
:rtype: float
"""
fpsi = 1 / (1 + (water_potential / parameters.water_potential_ref) ** parameters.n)
es_Ta = parameters.s_C * exp((parameters.s_B * Ta) / (parameters.s_A + Ta)) #: Saturated vapour pressure of the air (kPa), Ta in degree Celsius
V = RH * es_Ta #: Vapour pressure of the air (kPa)
VPDa = es_Ta - V
fvpd = 1 / (1 + VPDa / parameters.D0)
Cs = ambient_CO2 - An * (parameters.K_Cs / parameters.GB) #: CO2 concentration at organ surface (µmol mol-1 or Pa). From Prieto et al. (2012). GB in mol m-2 s-1
gamma = parameters.GAMMA0 * (1 + parameters.GAMMA1 * (Ta - parameters.T_ref) + parameters.GAMMA2 * ((Ta - parameters.T_ref) * (Ta - parameters.T_ref)))
gsw = parameters.GSMIN + (1.6 * parameters.m * Ag * fvpd * fpsi) / (Cs - gamma)
return gsw
STOMATAL_MODELS_MAPPING = {
'BWB': stomatal_conductance_BWB,
'Leuning': stomatal_conductance_Leuning,
'Tuzet': stomatal_conductance_Tuzet,
'hydraulics': stomatal_conductance_hydraulics
}
[docs]
def calculate_Ci(ambient_CO2, An, gsw):
"""
Calculates the internal CO2 concentration (Ci)
:param float ambient_CO2: air CO2 (µmol mol-1)
:param float An: net assimilation rate of CO2 (µmol m-2 s-1)
:param float gsw: stomatal conductance to water vapour (mol m-2 s-1)
:return: Ci (µmol mol-1)
:rtype: float
"""
Ci = ambient_CO2 - An * ((parameters.gsw_gs_CO2 / gsw) + (parameters.Ci_A / parameters.GB)) #: Intercellular concentration of CO2 (µmol mol-1)
# gsw and GB in mol m-2 s-1 so that (An * ((1.6/gs) + (1.37/parameters.GB)) is thus in µmol mol-1 as ambient_CO2
return Ci
def _f_temperature(pname, p25, T):
"""
Photosynthetic parameters relation to temperature
:param str pname: name of parameter
:param float p25: parameter value at 25 degree C
:param float T: organ temperature (degree C)
:return: p (parameter value at organ temperature)
:rtype: float
"""
Tk = T + parameters.KELVIN_DEGREE
deltaHa = parameters.PARAM_TEMP['deltaHa'][pname] #: Enthalpy of activation of parameter pname (kJ mol-1)
Tref = parameters.PARAM_TEMP['Tref']
f_activation = exp((deltaHa * (Tk - Tref)) / (parameters.R * 1E-3 * Tref * Tk)) #: Energy of activation (normalized to unity)
if pname in ('Vc_max', 'Jmax', 'TPU'):
deltaS = parameters.PARAM_TEMP['deltaS'][pname] #: entropy term of parameter pname (kJ mol-1 K-1)
deltaHd = parameters.PARAM_TEMP['deltaHd'][pname] #: Enthalpy of deactivation of parameter pname (kJ mol-1)
f_deactivation = (1 + exp((Tref * deltaS - deltaHd) / (Tref * parameters.R * 1E-3))) / (
1 + exp((Tk * deltaS - deltaHd) / (Tk * parameters.R * 1E-3))) #: Energy of deactivation (normalized to unity)
else:
f_deactivation = 1
p = p25 * f_activation * f_deactivation
return p
def _inhibition_by_NSC(NSC):
"""
Calculates the relative diminution of Ag due to inhibition by NSC. Adapted from Azcon-Bieto 1983
:param float NSC: Surfacic content of water-soluble carbohydrates (µmol C m-2)
:return: Relative diminution (dimensionless)
:rtype: float
"""
if NSC <= parameters.WSC_min:
return 0
else:
return min(parameters.Inhibition_max * (NSC - parameters.WSC_min) / (
parameters.K_Inhibition + NSC - parameters.WSC_min), 1)
[docs]
def calculate_photosynthesis(PAR, surfacic_nitrogen, NSC_Retroinhibition, surfacic_NSC, Ts, Ci):
"""
Computes photosynthesis rate following Farquhar's model with regulation by organ temperature and nitrogen content.
In this version, most of the parameters are derived from Braune et al. (2009) on barley and Evers et al. (2010) for N dependencies.
:param float PAR: PAR absorbed (µmol m-2 s-1)
:param float surfacic_nitrogen: surfacic nitrogen content(g m-2) including or not structural nitrogen depending on parameter.MODEL_VERSION
:param bool NSC_Retroinhibition: if True, Ag is inhibited by surfacic NSC (Non-Structural Carbohydrates).
:param float surfacic_NSC: surfacic content of NSC (Non-Structural Carbohydrates) (µmol C m-2).
:param float Ts: organ temperature (degree C)
:param float Ci: internal CO2 (µmol mol-1), Ci = 0.7*CO2air for the first iteration
:return: Ag (µmol m-2 s-1), An (µmol m-2 s-1), Rd (µmol m-2 s-1)
:rtype: (float, float, float)
"""
#: RuBisCO parameters dependence to temperature
Kc = _f_temperature('Kc', parameters.KC25, Ts)
Ko = _f_temperature('Ko', parameters.KO25, Ts)
Gamma = _f_temperature('Gamma', parameters.GAMMA25, Ts)
#: RuBisCO-limited carboxylation rate
Sna_Vcmax25 = parameters.PARAM_N['S_surfacic_nitrogen']['Vc_max25']
surfacic_nitrogen_min_Vcmax25 = parameters.PARAM_N['surfacic_nitrogen_min']['Vc_max25']
Vc_max25 = Sna_Vcmax25 * (surfacic_nitrogen - surfacic_nitrogen_min_Vcmax25) #: Relation between Vc_max25 and surfacic_nonstructural_nitrogen (µmol m-2 s-1)
Vc_max = _f_temperature('Vc_max', Vc_max25, Ts) #: Relation between Vc_max and temperature (µmol m-2 s-1)
Ac = (Vc_max * (Ci - Gamma)) / (Ci + Kc * (1 + parameters.O2 / Ko)) #: Rate of assimilation under Vc_max limitation (µmol m-2 s-1)
#: RuBP regeneration-limited carboxylation rate via electron transport
ALPHA = parameters.PARAM_N['S_surfacic_nitrogen']['alpha'] * surfacic_nitrogen + parameters.PARAM_N['beta'] #: Relation between ALPHA and surfacic_nitrogen (mol e- mol-1 photon)
Sna_Jmax25 = parameters.PARAM_N['S_surfacic_nitrogen']['Jmax25']
surfacic_nitrogen_min_Jmax25 = parameters.PARAM_N['surfacic_nitrogen_min']['Jmax25']
Jmax25 = Sna_Jmax25 * (surfacic_nitrogen - surfacic_nitrogen_min_Jmax25) #: Relation between Jmax25 and surfacic_nitrogen (µmol m-2 s-1)
Jmax = _f_temperature('Jmax', Jmax25, Ts) #: Relation between Jmax and temperature (µmol m-2 s-1)
J = ((Jmax + ALPHA * PAR) - sqrt((Jmax + ALPHA * PAR) ** parameters.J_expo - parameters.J_A * parameters.THETA * ALPHA * PAR * Jmax)) / (
parameters.J_B * parameters.THETA) #: Electron transport rate (Muller et al. (2005), Evers et al. (2010)) (µmol m-2 s-1)
Aj = (J * (Ci - Gamma)) / (parameters.Aj_A * Ci + parameters.Aj_B * Gamma) #: Rate of assimilation under RuBP regeneration limitation (µmol m-2 s-1)
#: Triose phosphate utilisation-limited carboxylation rate --- NOT USED, calculated just for information
Sna_TPU25 = parameters.PARAM_N['S_surfacic_nitrogen']['TPU25']
surfacic_nitrogen_min_TPU25 = parameters.PARAM_N['surfacic_nitrogen_min']['TPU25']
TPU25 = Sna_TPU25 * (surfacic_nitrogen - surfacic_nitrogen_min_TPU25) #: Relation between TPU25 and surfacic_nitrogen (µmol m-2 s-1)
TPU = _f_temperature('TPU', TPU25, Ts) #: Relation between TPU and temperature (µmol m-2 s-1)
Vomax = (Vc_max * Ko * Gamma) / (parameters.Vomax_A * Kc * parameters.O2) #: Maximum rate of Vo (µmol m-2 s-1) (µmol m-2 s-1)
Vo = (Vomax * parameters.O2) / (parameters.O2 + Ko * (1 + Ci / Kc)) #: Rate of oxygenation of RuBP (µmol m-2 s-1)
Ap = (1 - Gamma / Ci) * (parameters.Ap_A * TPU + Vo) #: Rate of assimilation under TPU limitation (µmol m-2 s-1).
# I think there was a mistake in the paper of Braune et al. (2009) where they wrote Ap = (1-Gamma/Ci)*(3*TPU) + Vo
# A more recent expression of Ap was given by S. v Caemmerer in her book (2000): AP = (3TPU * (Ci-Gamma))/(Ci-(1+3alpha)*Gamma),
# where 0 < alpha > 1 is the fraction of glycolate carbon not returned to the chloroplast, but I couldn't find any estimation of alpha for wheat
#: Gross assimilation rate (µmol m-2 s-1)
if NSC_Retroinhibition:
Ag = min(Ac, Aj) * (1 - _inhibition_by_NSC(surfacic_NSC))
else:
Ag = min(Ac, Aj, Ap)
#: Mitochondrial respiration rate of organ in light Rd (processes other than photorespiration)
Rdark25 = parameters.PARAM_N['S_surfacic_nitrogen']['Rdark25'] * (surfacic_nitrogen - parameters.PARAM_N['surfacic_nitrogen_min'][
'Rdark25']) #: Relation between Rdark25 (respiration in obscurity at 25 degree C) and surfacic_nitrogen (µmol m-2 s-1)
Rdark = _f_temperature('Rdark', Rdark25, Ts) #: Relation between Rdark and temperature (µmol m-2 s-1)
Rd = Rdark * (parameters.Rd_A + (1 - parameters.Rd_A) * parameters.Rd_B ** (PAR / parameters.Rd_C)) # Found in Muller et al. (2005), eq. 19 (µmol m-2 s-1)
#: Net C assimilation (µmol m-2 s-1)
if Ag <= 0: # Occurs when Ci is lower than Gamma or when (surfacic_nitrogen - surfacic_nitrogen_min)<0, in these cases there is no net assimilation (Farquhar, 1980; Caemmerer, 2000)
Ag, An = 0, 0
else:
An = Ag - Rd
return Ag, An, Rd
[docs]
def calculate_surfacic_nitrogen(nitrates, amino_acids, proteins, Nstruct, green_area):
"""Surfacic content of nitrogen
:param float nitrates: amount of nitrates (µmol N)
:param float amino_acids: amount of amino_acids (µmol N)
:param float proteins: amount of proteins (µmol N)
:param float Nstruct: structural N (g)
:param float green_area: green area (m-2)
:return: Surfacic nitrogen (g m-2)
:rtype: float
"""
mass_N_tot = (nitrates + amino_acids + proteins) * 1E-6 * parameters.N_MOLAR_MASS + Nstruct
return mass_N_tot / green_area
[docs]
def calculate_surfacic_nonstructural_nitrogen(nitrates, amino_acids, proteins, green_area):
"""Surfacic content of non-structural nitrogen
:param float nitrates: amount of nitrates (µmol N)
:param float amino_acids: amount of amino_acids (µmol N)
:param float proteins: amount of proteins (µmol N)
:param float green_area: green area (m-2)
:return: Surfacic non-structural nitrogen (g m-2)
:rtype: float
"""
mass_N_tot = (nitrates + amino_acids + proteins) * 1E-6 * parameters.N_MOLAR_MASS
return mass_N_tot / green_area
[docs]
def calculate_surfacic_photosynthetic_proteins(proteins, green_area):
"""Surfacic content of photosynthetic proteins
:param float proteins: amount of proteins (µmol N)
:param float green_area: green area (m-2)
:return: Surfacic non-structural nitrogen (g m-2)
:rtype: float
"""
mass_N_prot = proteins * 1E-6 * parameters.N_MOLAR_MASS
return mass_N_prot / green_area
[docs]
def calculate_surfacic_nonstructural_nitrogen_Farquhar(surfacic_photosynthetic_proteins):
"""Estimate of non structural SLN used in Farquhar
:param float surfacic_photosynthetic_proteins: surfacic proteins content (µmol N m-2)
:return: Surfacic non-structural nitrogen (g m-2)
:rtype: float
"""
return surfacic_photosynthetic_proteins * parameters.Psurf_to_SLNnonstruct
[docs]
def calculate_surfacic_WSC(sucrose, starch, fructan, green_area):
"""Surfacic content of water soluble carbohydrates # TODO: rename because starch is not water-soluble.
:param float sucrose: amount of sucrose (µmol C)
:param float starch: amount of starch (µmol C)
:param float fructan: amount of fructan (µmol C)
:param float green_area: green area (m-2)
:return: Surfacic content of water soluble carbohydrates (µmol C m-2)
:rtype: float
"""
return (sucrose + starch + fructan) / green_area
