In [1]:
%%capture storage
# The above redirects all output of the below commands to the variable 'storage' instead of displaying it.
# It can be viewed by typing: 'storage()'
# Setting up worksheet and importing equations from other worksheets
load('temp/Worksheet_setup.sage')
load_session('temp/leaf_enbalance_eqs.sobj')
fun_loadvars(dict_vars) # any variables defined using var2() are saved with their attributes in dict_vars. Here we re-define them based on dict_vars from the other worksheet.

Equations to compute stomatal conductance based on sizes and densities of stomata

Equation numbers in comments refer to Lehmann & Or, 2013

Definitions of additional variables


In [2]:
# Following Lehmann_2015_Effects
var2('B_l', 'Boundary layer thickness', meter, domain1='positive')
var2('g_sp', 'Diffusive conductance of a stomatal pore', mole/second/meter^2, domain1='positive')
var2('r_sp', 'Diffusive resistance of a stomatal pore', meter^2*second/mole, domain1='positive')
var2('r_vs', 'Diffusive resistance of a stomatal vapour shell', meter^2*second/mole, domain1='positive')
var2('r_bwmol', 'Leaf BL resistance in molar units',  meter^2*second/mole, domain1='positive', latexname = 'r_{bw,mol}')
var2('r_end', 'End correction, representing resistance between evaporating sites and pores', meter^2*second/mole, domain1='positive')
var2('A_p', 'Cross-sectional pore area', meter^2, domain1='positive')
var2('d_p', 'Pore depth', meter, domain1='positive')
var2('V_m', 'Molar volume of air', meter^3/mole, domain1='positive')
var2('n_p', 'Pore density', 1/meter^2, domain1='positive')
var2('r_p', 'Pore radius (for ellipsoidal pores, half the pore width)', meter, domain1='positive')
var2('l_p', 'Pore length', meter, domain1='positive')
var2('k_dv', 'Ratio $D_{va}/V_m$', mole/meter/second, domain1='positive', latexname='k_{dv}')
var2('s_p', 'Spacing between stomata', meter, domain1='positive')
var2('F_p', 'Fractional pore area (pore area per unit leaf area)', meter^2/meter^2, domain1='positive')


Out[2]:
F_p

In [3]:
docdict[V_m]


Out[3]:
'Molar volume of air'

Equations from Lehmann & Or, 2013


In [4]:
# Eqn1
eq_kdv = k_dv == D_va/V_m
print units_check(eq_kdv)
eq_gsp_A = g_sp == A_p/d_p*k_dv*n_p
print units_check(eq_gsp_A)
eq_Ap_rp = A_p == pi*r_p^2
print units_check(eq_Ap_rp)
eq_rsp_A = r_sp == 1/eq_gsp_A.rhs()
print units_check(eq_rsp_A)
eq_rsp_rp = r_sp == 1/eq_gsp_A.rhs().subs(eq_Ap_rp)
print units_check(eq_rsp_rp)


mole/(meter*second) == mole/(meter*second)
mole/(meter^2*second) == mole/(meter^2*second)
meter^2 == meter^2
meter^2*second/mole == meter^2*second/mole
meter^2*second/mole == meter^2*second/mole

In [5]:
# Eq. 2(c)
eq_rvs_B = r_vs == (1/(4*r_p) - 1/(pi*s_p))*1/(k_dv*n_p)
units_check(eq_rvs_B)


Out[5]:
meter^2*second/mole == meter^2*second/mole

In [6]:
# Eq. 2(d)
eq_rvs_S = r_vs == (1/(2*pi*r_p) - 1/(pi*s_p))*1/(k_dv*n_p)
units_check(eq_rvs_B)


Out[6]:
meter^2*second/mole == meter^2*second/mole

In [7]:
# Below Eq. 3(b)
assume(n_p>0)
eq_sp_np = s_p == 1/sqrt(n_p)
units_check(eq_sp_np)


Out[7]:
meter == sqrt(n_p)/sqrt(n_p/meter^2)

In [8]:
# Eq. 2(a)
eq_rend_rp = r_end == 1/(4*r_p*k_dv*n_p)
units_check(eq_rend_rp)


Out[8]:
meter^2*second/mole == meter^2*second/mole

In [9]:
# Eq. 2(b)
eq_rend_PW = r_end == log(4*(l_p/2)/r_p)/(2*pi*l_p/2*k_dv*n_p)
units_check(eq_rend_PW)


Out[9]:
meter^2*second/mole == meter^2*second/mole

In [10]:
eq_gswmol = g_swmol == 1/(r_end + r_sp + r_vs)
units_check(eq_gswmol)


Out[10]:
mole/(meter^2*second) == mole/(meter^2*second)

It actually does not seem to make much sense to add r_end to the resistances, as it does not reflect the geometry of the inter-cellular air space! Also eq_rvs_B is the correct one to use for stomata, as eq_rvs_S is valid for droplets, not holes!


In [11]:
assume(s_p>0)
eq_np_sp = solve(eq_sp_np, n_p)[0]
units_check(eq_np_sp)


Out[11]:
meter^(-2) == meter^(-2)

In [12]:
# Eq. 5
eq_rbwmol = r_bwmol == B_l/k_dv
units_check(eq_rbwmol)


Out[12]:
meter^2*second/mole == meter^2*second/mole

In [13]:
docdict[k_dv]


Out[13]:
'Ratio $D_{va}/V_m$'

Below, we compare our system of equations to Eq. 7a in Lehmann & Or (2015), i.e. we check if $\frac{r_{bwmol}}{r_{tot}} = 1/(1 + s_p^2/B_l*(1/(2*r_p) + d_p/(\pi*r_p^2)))$:


In [14]:
# Verification of Eq. 7a
eq1 = (r_bwmol/(2*r_end + r_sp + r_bwmol)).subs(eq_rend_rp, eq_rsp_rp, eq_rvs_B, eq_rbwmol).subs(eq_np_sp)
eq1.simplify_full().show()
eq7a = 1/(1 + s_p^2/B_l*(1/(2*r_p) + d_p/(pi*r_p^2)))
eq7a.show()
(eq1 - eq7a).simplify_full()


Out[14]:
0

Additional equations for calculating conductances of laser perforated foils


In [15]:
# Pore area for circular pores
eq_Ap = A_p == pi*r_p^2
units_check(eq_Ap)


Out[15]:
meter^2 == meter^2

In [16]:
# To deduce pore radius from pore area, assuming circular pore:
eq_rp_Ap = solve(eq_Ap_rp, r_p)[0]
units_check(eq_rp_Ap)


Out[16]:
meter == sqrt(A_p*meter^2/pi)/sqrt(A_p/pi)

In [17]:
# For regular pore distribution we can derive the pore density from porosity and pore area
eq_np = n_p == F_p/A_p

In [18]:
# Conversion from mol/m2/s to m/s, i.e. from g_swmol to g_sw
eq_gsw_gswmol = solve(eq_gtwmol_gtw(g_tw = g_sw, g_twmol = g_swmol), g_sw)[0]
units_check(eq_gsw_gswmol)


Out[18]:
meter/second == meter/second

In [19]:
# Simplification, neglecting effect of leaf-air temperature difference
eq_gsw_gswmol_iso = eq_gsw_gswmol(T_l = T_a).simplify_full()
units_check(eq_gsw_gswmol_iso)


Out[19]:
meter/second == meter/second

In [20]:
# Saving session so that it can be loaded using load_session().
save_session('temp/stomatal_cond_eqs.sobj')

Table of symbols


In [21]:
# Creating dictionary to substitute names of units with shorter forms
var('m s J Pa K kg mol')
subsdict = {meter: m, second: s, joule: J, pascal: Pa, kelvin: K, kilogram: kg, mole: mol}
var('N_Re_L N_Re_c N_Le N_Nu_L N_Gr_L N_Sh_L')
dict_varnew = {Re: N_Re_L, Re_c: N_Re_c, Le: N_Le, Nu: N_Nu_L, Gr: N_Gr_L, Sh: N_Sh_L}
dict_varold = {v: k for k, v in dict_varnew.iteritems()}
variables = sorted([str(variable.subs(dict_varnew)) for variable in udict.keys()],key=str.lower)
tableheader = [('Variable', 'Description (value)', 'Units')]
tabledata = [('Variable', 'Description (value)', 'Units')]
for variable1 in variables:
    variable2 = eval(variable1).subs(dict_varold)
    variable = str(variable2)
    tabledata.append((eval(variable),docdict[eval(variable)],fun_units_formatted(variable)))

table(tabledata, header_row=True)


Out[21]:
Variable Description (value) Units
Cross-sectional pore area m
Fraction of one-sided leaf area covered by stomata (1 if stomata are on one side only, 2 if they are on both sides) 1
Fraction of projected area exchanging sensible heat with the air (2) 1
Thermal diffusivity of dry air m s
Boundary layer thickness m
Specific heat of dry air (1010) J K kg
Concentration of water in the free air mol m
Concentration of water in the leaf air space mol m
Pore depth m
Binary diffusion coefficient of water vapour in air m s
Latent heat flux from leaf J m s
Transpiration rate in molar units mol m s
Longwave emmissivity of the leaf surface (1.0) 1
Fractional pore area (pore area per unit leaf area) 1
Gravitational acceleration (9.81) m s
Boundary layer conductance to water vapour m s
Boundary layer conductance to water vapour mol m s
Diffusive conductance of a stomatal pore mol m s
Stomatal conductance to water vapour m s
Stomatal conductance to water vapour mol m s
Total leaf conductance to water vapour m s
Total leaf layer conductance to water vapour mol m s
Average 1-sided convective transfer coefficient J K m s
Sensible heat flux from leaf J m s
Thermal conductivity of dry air J K m s
Ratio mol m s
Characteristic length scale for convection (size of leaf) m
Pore length m
Latent heat of evaporation (2.45e6) J kg
Molar mass of nitrogen (0.028) kg mol
Molar mass of oxygen (0.032) kg mol
Molar mass of water (0.018) kg mol
Grashof number 1
Lewis number 1
Nusselt number 1
Pore density m
Critical Reynolds number for the onset of turbulence 1
Reynolds number 1
Sherwood number 1
Kinematic viscosity of dry air m s
Air pressure Pa
Partial pressure of nitrogen in the atmosphere Pa
Partial pressure of oxygen in the atmosphere Pa
Vapour pressure in the atmosphere Pa
Saturation vapour pressure at air temperature Pa
Vapour pressure inside the leaf Pa
Prandtl number (0.71) 1
Boundary layer resistance to water vapour, inverse of s m
Leaf BL resistance in molar units s m mol
End correction, representing resistance between evaporating sites and pores s m mol
Longwave radiation away from leaf J m s
Molar gas constant (8.314472) J K mol
Pore radius (for ellipsoidal pores, half the pore width) m
Solar shortwave flux J m s
Diffusive resistance of a stomatal pore s m mol
Stomatal resistance to water vapour, inverse of s m
Total leaf resistance to water vapour, s m
Diffusive resistance of a stomatal vapour shell s m mol
Density of dry air kg m
Density of air at the leaf surface kg m
Spacing between stomata m
Stefan-Boltzmann constant (5.67e-8) J K m s
Air temperature K
Leaf temperature K
Radiative temperature of objects surrounding the leaf K
Molar volume of air m mol
Wind velocity m s