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.
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]:
In [3]:
docdict[V_m]
Out[3]:
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)
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]:
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]:
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]:
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]:
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]:
In [10]:
eq_gswmol = g_swmol == 1/(r_end + r_sp + r_vs)
units_check(eq_gswmol)
Out[10]:
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]:
In [12]:
# Eq. 5
eq_rbwmol = r_bwmol == B_l/k_dv
units_check(eq_rbwmol)
Out[12]:
In [13]:
docdict[k_dv]
Out[13]:
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]:
In [15]:
# Pore area for circular pores
eq_Ap = A_p == pi*r_p^2
units_check(eq_Ap)
Out[15]:
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]:
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]:
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]:
In [20]:
# Saving session so that it can be loaded using load_session().
save_session('temp/stomatal_cond_eqs.sobj')
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]: