R/StableIsotopes.R

In griffithdan/r3PG: An R Implementation of the 3PG (Physiological Principles Predicting Growth) Model

#' Calculations for stable isotopes
#'
#' Internal functions for isotope calculations

calc_d13c <- function(T_av, CaMonthly, D13Catm, elev, GPPmolc, month, canopy_conductance, RGcGW, D13CTissueDif, aFracDiffu, bFracRubi){
    # calculating d13C by Liang Wei

    # Air pressure, kpa
    AirPressure = 101.3 * exp(-1 * elev / 8200)
    # Convert Unit of Atmospheric C, a ppm to part/part
    AtmCa = CaMonthly * 0.000001
    # Canopy conductance for water vapor in mol/m2s, unit conversion
    GwMol = canopy_conductance * 44.6 * (273.15 / (273.15 + T_av)) * (AirPressure / 101.3)
    # Canopy conductance for CO2 in mol/m2s
    GcMol = GwMol * RGcGW

    # GPP per second. Unit: mol/m2 s. GPPmolc divide by 24 hour/day and 3600 s/hr
    GPPmolsec = GPPmolc / (get_days_in_month(month) * 24 * 3600)

    # Calculating monthly average intercellular CO2 concentration. Ci = Ca - A/g
    InterCi = AtmCa - GPPmolsec / GcMol
    InterCiPPM = InterCi * 1000000

    # Calculating monthly d13C of new photosynthate, = d13Catm- a-(b-a) (ci/ca)
    D13CNewPS = D13Catm - aFracDiffu - (bFracRubi - aFracDiffu) * (InterCi / AtmCa)
    D13CTissue = D13CNewPS + D13CTissueDif
    d13c_list <- list(D13CTissue = D13CTissue, InterCiPPM = InterCiPPM)
    return(d13c_list)
}

d18O <- function(T_av, VPD, d18Osrc, canopy_transpiration_sec){
    #"""calculate d18O of leafwater and cellulose from tair, rh, and d18O of source water and vapor
    #"""
    #ratio between 18O and 16O (18O/16O) in VSMOW (Vienna-Standard Mean Ocean Water)
    r_vsmow = 0.0020052
    d18Osource = d18Osrc
    r_source = r_vsmow * (d18Osource / 1000. + 1) #;print d18Osource, r_source
    Tc = T_av
    Tk = Tc + 273.15 #[K]
    # now calculate the saturation and ambient vapor pressures based on rel. hum. and Tair
    esat = 6.112 * exp( 17.67 * Tc / (Tc + 243.5)) #(mbar) SVP versus Temp relationship from Jacobson 1999
    #this assumes that the VPD data being read in are in mbar, as values peak at ~15 (so could not be kPa)
    #also Wei et al. PC&E gives this unit in Table A1 and in the Appendix for VPD modifier function
    ea = esat - VPD
    epsilon_eq = exp(1137. / Tk ^ 2 - 0.4156 / Tk - 0.0020667) - 1.
    alpha_eq = exp(1137. / Tk ^ 2 - 0.4156 / Tk - 0.0020667); #print alpha_eq
    epsilon_eq = 1000. * epsilon_eq #[per mil] #;print(round(epsilon_eq,2))
    epsilon_eq2 = 2.644 - 3.206 * (1.e3 / Tk) + 1.534 * (1e6 / Tk ^ 2) #[per mil] ;print(round(epsilon_eq2,2))  #from Bottinga and Craig 1969
    alpha_eq2 = epsilon_eq2 / 1000. + 1. #print alpha_eq2
    d18Ovapor = d18Osource - epsilon_eq #[per mil] ;print d18Ovapor
    r_vapor = r_source / alpha_eq
    n = 1.0 #The exponent n depends on turbulence intensity - ranges from 0.5 for fully turbulent conditions to 1.0 for pure molecular diffusion [Mathieu and Bariac,1996].
    alpha_k = (1. / 0.97229) ^ n
    epsilon_k = 1000. * (alpha_k - 1) #[per mil]
    gs = 0.150 #[mol/m2/s] leaf stomatal conductance to water vapor
    gb = 1.42 #[mol/m2/s] leaf boundary conductance to water vapor
    epsilon_k2_weighted = (28. * 1 / gs + 22. * 1 / gb)/(1. / gs + 1. / gb) #[per mil] from Luz et al. 2009 (table 2), where gb is constant at 1.42 mol/m2/s
    alpha_k2 = epsilon_k2_weighted / 1000. + 1 # alpha = epsilon/1000 +1
    d18Oleaf = epsilon_eq + (1. - ea / esat) * (d18Osource + epsilon_k2_weighted) + (ea / esat) * d18Ovapor #[per mil]
    r_leaf = alpha_eq * (alpha_k2 * r_source * ((esat - ea) / esat) + r_vapor * (ea / esat))
    d18Oleaf_check = 1000. * (r_leaf / r_vsmow - 1.) #print round(d18Oleaf_check,1)
    #D18Oleaf = (d18Oleaf/1000. - d18Osource/1000.)/(1 + d18Osource/1000.) #expressed as enrichment above source water, i.e. 'big delta'
    D18Ovapor = (d18Ovapor - d18Osource)/(1. + d18Osource / 1000.) #expressed as enrichment above source water, i.e. 'big delta'
    D18Oleaf = epsilon_eq + epsilon_k2_weighted + (D18Ovapor - epsilon_k2_weighted) * (ea / esat)
    # L at steady state varies from ~.01m to ~0.2m (see Fig. 2 from Ferrio)
    L = 0.01 #[m] scaled effective path length (m) for water movement from the veins to the site of evaporation (Kleaf = 31.5 * L**(-0.43),
    C = 55.56e3 #[mol/m3] the molar concentration of water
    a_d = 1. / 1.026; a1 = 100.e-9; a2 = 577.; a3 = 145.
    D = a_d * a1 * exp(-a2 / (Tk - a3))
    ## D should be close to 2.66e-9 for H218O at 32.5 C (see Cuntz et al. 2007, p.897)
    #need to import ET from the model rather than the fixed rate here - add as canopy_transpiration_sec
    #0.05 #[mol/m2/s] leaf transpiration rate
    E = canopy_transpiration_sec
    peclet = (E * L) / (C * D)
    D18Oleaf_peclet = D18Oleaf * ((1 - exp(-peclet)) / peclet)
    epsilon_cell = 27. #[per mil] fractionation factor between carbonyl oxygen and water
    fO = 0.42 #[unitless] fraction of cellulose oxygen that exchanges atoms with xylem water during synthesis
    d18Ocell = fO * (d18Osource + epsilon_cell) + (1. - fO) * (d18Oleaf + epsilon_cell)
    d18Ocell_peclet = fO * (d18Osource + epsilon_cell) + (1. - fO) * (D18Oleaf_peclet + epsilon_cell)
    d18O_list <- list(d18Oleaf = d18Oleaf, d18Ocell = d18Ocell, d18Ocell_peclet = d18Ocell_peclet)
    return(d18O_list)
}