leaf_enbalance_eqs.ipynb

Kernel: SageMath 7.3

Leaf energy balance equations

Based on the following paper: Schymanski, S.J. and Or, D. (2016): Leaf-scale experiments reveal important omission in the Penman-Monteith equation. Hydrology and Earth System Sciences Discussions, p.1–33. doi: 10.5194/hess-2016-363.

Author: Stan Schymanski ([email protected])

Note: This worksheet is prepared for the open source software sage. It relies on definitions provided in Worksheet Worksheet_setup.

%%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()'

# Loading modules, settings and custom functions
load('temp/Worksheet_setup.sage')

Definition of variables

All variables are also listed in a Table at the end of this document.

var2('alpha_a', 'Thermal diffusivity of dry air',meter^2/second)
var2('a_s', 'Fraction of one-sided leaf area covered by stomata (1 if stomata are on one side only, 2 if they are on both sides)', meter^2/meter^2)
var2('a_sh', 'Fraction of projected area exchanging sensible heat with the air (2)', meter^2/meter^2, value = 2, latexname = 'a_{sh}')
var2('c_pa', 'Specific heat of dry air (1010) ', joule/kilogram/kelvin, latexname = 'c_{pa}', value = 1010)
#var2('c_pamol', 'Molar specific heat of dry air (29.19) ', joule/mole/kelvin, latexname = 'c_{pa,mol}', value = 29.19) # https://en.wikipedia.org/wiki/Heat_capacity#Specific_heat_capacity
var2('C_wa', 'Concentration of water in the free air ',mole/meter^3, latexname = 'C_{wa}')
var2('C_wl', 'Concentration of water in the leaf air space ',mole/meter^3, latexname = 'C_{wl}')
var2('D_va', 'Binary diffusion coefficient of water vapour in air',meter^2/second, latexname = 'D_{va}')
var2('E_lmol', 'Transpiration rate in molar units',mole/second/meter^2,latexname='E_{l,mol}')
var2('E_l', 'Latent heat flux from leaf',joule/second/meter^2)
var2('epsilon_l', 'Longwave emmissivity of the leaf surface (1.0)', value = 1, units = 1/1)
var2('g', 'Gravitational acceleration (9.81)', meter/second^2, value= 9.81)
var2('g_bw', 'Boundary layer conductance to water vapour ',meter/second,latexname='g_{bw}')
var2('g_bwmol', 'Boundary layer conductance to water vapour ',mole/meter^2/second,latexname='g_{bw,mol}')
var2('Gr', 'Grashof number', latexname = 'N_{Gr_L}')
var2('g_sw', 'Stomatal conductance to water vapour',meter/second,latexname='g_{sw}')
var2('g_swmol', 'Stomatal conductance to water vapour',mole/meter^2/second,latexname='g_{sw,mol}')
var2('g_tw', 'Total leaf conductance to water vapour',meter/second,latexname='g_{tw}')
var2('g_twmol', 'Total leaf layer conductance to water vapour',mole/meter^2/second,latexname='g_{tw,mol}')
var2('h_c', 'Average 1-sided convective transfer coefficient', joule/meter^2/second/kelvin)
var2('H_l', 'Sensible heat flux from leaf',joule/second/meter^2)
var2('k_a', 'Thermal conductivity of dry air',joule/second/meter/kelvin)
var2('lambda_E', 'Latent heat of evaporation (2.45e6)',joule/kilogram,value = 2.45e6)
var2('Le', 'Lewis number', latexname = 'N_{Le}')
var2('L_l', 'Characteristic length scale for convection (size of leaf)',meter)
var2('M_N2', 'Molar mass of nitrogen (0.028)',kilogram/mole,value = 0.028)
var2('M_O2', 'Molar mass of oxygen (0.032)',kilogram/mole,value = 0.032)
var2('M_w', 'Molar mass of water (0.018)',kilogram/mole,value = 0.018)
var2('nu_a', 'Kinematic viscosity of dry air',meter^2/second)
var2('Nu', 'Nusselt number',latexname = 'N_{Nu_L}')
var2('P_a', 'Air pressure',pascal)
var2('Pr', 'Prandtl number (0.71)',latexname = 'N_{Pr}',value = 0.71)
var2('P_N2', 'Partial pressure of nitrogen in the atmosphere', pascal, latexname = 'P_{N2}')
var2('P_O2', 'Partial pressure of oxygen in the atmosphere', pascal, latexname = 'P_{O2}')
var2('P_wa', 'Vapour pressure in the atmosphere', pascal, latexname = 'P_{wa}')
var2('P_was', 'Saturation vapour pressure at air temperature', pascal, latexname = 'P_{was}')
var2('P_wl', 'Vapour pressure inside the leaf', pascal, latexname = 'P_{wl}')
var2('r_bw', 'Boundary layer resistance to water vapour, inverse of $g_{bw}$', second/meter, latexname = 'r_{bw}')
var2('r_sw', 'Stomatal resistance to water vapour, inverse of $g_{sw}$', second/meter, latexname = 'r_{sw}')
var2('r_tw', 'Total leaf resistance to water vapour, $r_{bv} + r_{sv}$', second/meter, latexname = 'r_{tw}')
var2('Re_c', 'Critical Reynolds number for the onset of turbulence',latexname='N_{Re_c}')
var2('Re', 'Reynolds number',latexname='N_{Re_L}')
var2('rho_a', 'Density of dry air', kilogram/meter^3)
var2('rho_al', 'Density of air at the leaf surface', kilogram/meter^3)
var2('R_ll', 'Longwave radiation away from leaf',joule/second/meter^2, latexname = 'R_{ll}')
var2('R_mol', 'Molar gas constant (8.314472)',joule/mole/kelvin,latexname='R_{mol}', value = 8.314472)
var2('R_s', 'Solar shortwave flux',joule/second/meter^2)
var2('Sh', 'Sherwood number',latexname = 'N_{Sh_L}')
var2('sigm', 'Stefan-Boltzmann constant (5.67e-8)',joule/second/meter^2/kelvin^4,'real',latexname='\\sigma',value=5.67e-8)
var2('T_a', 'Air temperature',kelvin)
var2('T_l', 'Leaf temperature',kelvin)
var2('T_w', 'Radiative temperature of objects surrounding the leaf', kelvin)
var2('v_w', 'Wind velocity',meter/second)

v_w

Mathematical derivations

Leaf energy balance

The material below follows closely derivations published previously \citep{schymanski_stomatal_2013}, with re-organisation of equations for greater consistence with the present paper.

The leaf energy balance is determined by the dominant energy fluxes between the leaf and its surroundings, including radiative, sensible, and latent energy exchange (linked to mass exchange).

In this study we focus on steady-state conditions, in which the energy balance can be written as:

{eq_Rs_enbal}

$$R_s = R_{ll} + H_l + E_l$$

where $R_s$ absorbed short wave radiation, $R_{ll}$ is the net longwave balance, i.e. the emitted minus the absorbed, $H_l$ is the sensible heat flux away from the leaf and $E_l$ is the latent heat flux away from the leaf. In the above, extensive variables are defined per unit leaf area. Following our previous work \citep{schymanski_stomatal_2013}, this study considers spatially homogeneous planar leaves, i.e. homogenous illumination and a negligible temperature gradient between the two sides of the leaf. The net longwave emission is represented by the difference between blackbody radiation at leaf temperature ($T_l$) and that at the temperature of the surrounding objects ($T_w$) \citep{Monteith_principles_2007}:

{eq_Rll}

$$R_{ll} = a_{sH} \epsilon_l \sigma (T_{l}^{4} - T_{w}^{4})$$

where $a_{sH}$ is the fraction of projected leaf area exchanging radiative and sensible heat (2 for a planar leaf, 1 for a soil surface), $\epsilon_l$ is the leaf's longwave emmissivity ($\approx 1$) and $\sigma$ is the Stefan-Boltzmann constant. Total convective heat transport away from the leaf is represented as:

{eq_Hl}

$$H_l = a_{sH} h_c (T_l - T_a)$$

where $h_c$ is the average one-sided convective heat transfer coefficient, determined by properties of the leaf boundary layer.

Latent heat flux ($E_l$, W~m$^{-2}$) is directly related to the transpiration rate ($E_{l,mol}$) by:

{eq_El}

$$E_l = E_{l,mol} M_w \lambda_E$$

where $M_w$ is the molar mass of water and $\lambda_E$ the latent heat of vaporisation. $E_{l,mol}$ (mol~m$^{-2}$~s$^{-1}$) was computed in molar units as a function of the concentration of water vapour within the leaf ($C_{wl}$, mol m$^{-3}$) and in the free air ($C_{wa}$, mol m$^{-3}$) \citep[Eq. 6.8]{incropera_fundamentals_2006}:

{eq_Elmol}

$$E_{l,mol} = g_{tw}(C_{wl} - C_{wa})$$

where $g_{tw}$ (m~s$^{-1}$) is the total leaf conductance for water vapour, dependent on stomatal ($g_{sw}$) and boundary layer conductance ($g_{bw}$) in the following way:

{eq_gtw}

$$g_{tw} = \frac{1}{\frac{1}{g_{sw}} + \frac{1}{g_{bw}}}$$

eq_Rs_enbal = R_s == R_ll + H_l + E_l
units_check(eq_Rs_enbal).simplify_full()

kilogram/second^3 == kilogram/second^3

eq_Rll = R_ll == a_sh*sigm*epsilon_l*(T_l^4 - T_w^4)
units_check(eq_Rll).simplify_full()

kilogram/second^3 == kilogram/second^3

eq_Hl = H_l == a_sh*h_c*(T_l - T_a)
units_check(eq_Hl).simplify_full()

kilogram/second^3 == kilogram/second^3

eq_El = E_l == E_lmol*M_w*lambda_E
units_check(eq_El).simplify_full()

kilogram/second^3 == kilogram/second^3

ustr = str(units_check(M_w*lambda_E))
print str((M_w*lambda_E).subs(cdict)) + ustr

44100.0000000000kilogram*meter^2/(mole*second^2)

eq_Elmol = E_lmol == g_tw*(C_wl - C_wa)
units_check(eq_Elmol).simplify_full()

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

eq_gtw = g_tw == 1/(1/g_sw + 1/g_bw)
eq_gtw.simplify_full().show()
units_check(eq_gtw).simplify_full()

meter/second == meter/second

Boundary layer conductance to water vapour

The total leaf conductance to water vapour is determined by the boundary layer and stomatal conductances and equal to 1 over the sum of their respectife resistances ($g_{tw} = 1/(r_{sw} + r_{bw})$. The boundary layer conductance for water vapour is equivalent to the mass transfer coefficient for a wet surface \cite[Eq. 7.41]{incropera_fundamentals_2006}:

{eq_gbw}

$$g_{bw} = N_{Sh_L} D_{va}/L_l$$

where $N_{Sh_L}$ is the dimensionless Sherwood number and $D_{va}$ is the diffusivity of water vapour in air. If the convection coefficient for heat is known, the one for mass ($g_{bw}$) can readily be calculated from the relation \cite[Eq. 6.60]{incropera_fundamentals_2006}:

{eq_gbw_hc}

$$g_{bw} = \frac{a_s h_c}{\rho_a c_{pa} N_{Le}^{1-n}}$$

where $a_s$ is the fraction of one-sided transpiring surface area in relation to the surface area for sensible heat exchange, $c_{pa}$ is the constant-pressure heat capacity of air, $n$ is an empirical constant ($n=1/3$ for general purposes) and $N_{Le}$ is the dimensionless Lewis number, defined as \cite[Eq. 6.57]{incropera_fundamentals_2006}:

{eq_Le}

$$N_{Le} = \alpha_a/D_{va}$$

where $\alpha_a$ is the thermal diffusivity of air. The value of $a_s$ was set to 1 for leaves with stomata on one side only, and to 2 for stomata on both sides. Other values could be used for leaves only partly covered by stomata.

eq_gbw = g_bw == D_va*Sh/L_l
units_check(eq_gbw).simplify_full()
eq_gbw_hc = g_bw == a_s*h_c/(rho_a*c_pa*Le^(2/3))
print units_check(eq_gbw_hc).simplify_full()
eq_Le = Le == alpha_a/D_va
units_check(eq_Le).simplify_full()

meter/second == meter/second

1 == 1

Effect of leaf temperature on the leaf-air vapour concentration gradient

\label{sec_Tleaf} The concentration difference in Eq. eq_Elmol is a function of the temperature and the vapour pressure differences between the leaf and the free air. Assuming that water vapour behaves like an ideal gas, we can express its concentration as:

{eq_Cwl}

$$C_{wl} = \frac{P_{wl}}{R_{mol} T_l}$$

where $P_{wl}$ is the vapour pressure inside the leaf, $R_{mol}$ is the universal gas constant and $T_l$ is leaf temperature. A similar relation holds for the vapour concentration in free air, $C_{wa} = P_{wa}/(R_{mol} T_l)$. In this study the vapour pressure inside the leaf is assumed to be the saturation vapour pressure at leaf temperature, which is computed using the Clausius-Clapeyron relation \citep[Eq. B.3]{hartmann_global_1994}:

{eq_Pwl}

$$P_{wl} = 611 \exp \left( \frac{\lambda_E M_w}{R_{mol}} \left( \frac{1}{273} - \frac{1}{T_l} \right) \right)$$

where $\lambda_E$ is the latent heat of vaporisation and $M_w$ is the molar mass of water.

Note that the dependence of the leaf-air water concentration difference ($C_{wl} - C_{wa}$) in Eq. eq_Cwl is very sensitive to leaf temperature. For example if the leaf temperature increases by 5K relative to air temperature, $C_{wl} - C_{wa}$ would double, while if leaf temperature decreased by 6K, $C_{wl} - C_{wa}$ would go to 0 at 70% relative humidity.

 

 

eq_Cwl = C_wl == P_wl/(R_mol*T_l)
print units_check(eq_Cwl).simplify_full()
eq_Pwl = P_wl == 611*exp(lambda_E*M_w/R_mol*(1/273 - 1/T_l))
show(eq_Pwl)

mole/meter^3 == mole/meter^3

Concentration or vapour pressure gradient driving transpiration?

Note that $E_{l,mol}$ is commonly expressed as a function of the vapour pressure difference between the free air ($P_{wa}$) and the leaf ($P_{wl}$), in which the conductance ($g_{tw,mol}$) is expressed in molar units (mol~m$^{-2}$~s$^{-1}$):

{eq_Elmol_conv}

$$E_{l,mol} = g_{tw,mol} \frac{P_{wl} - P_{wa}}{P_a}$$

For $P_{wl} = P_{wa}$, Eq. eq_Elmol can still give a flux, whereas Eq. eq_Elmol_conv gives zero flux. This is because the concentrations of vapour in air (mol~m$^{-3}$) can differ due to differences in tempertaure, even if the partial vapour pressures are the same (see Eq. eq_Cwl). Therefore, the relation between $g_{tw}$ and $g_{v,mol}$ has an asymptote at the equivalent temperature. It can be obtained by combining Eqs. eq_Elmol and eq_Elmol_conv and solving for $g_{tw,mol}$:

{eq_gtwmol_gtw}

$$g_{tw,mol} = g_{tw}\frac{P_a(P_{wa} T_l - P_{wl} T_a)}{(P_{wa} - P_{wl}) R_{mol} T_a T_l)}$$

For $T_l = T_a$, the relation simplifies to:

{eq_gtwmol_gtw_iso}

$$g_{tw,mol} = g_{tw}\frac{P_a }{R_{mol} T_a}$$

which, for typical values of $P_a$ and $T_a$ amounts to $g_{tw,mol}\approx40$~mol~m$^{-3}g_{tw}$. For all practical purposes, we found that Eqs. eq_Elmol and eq_Elmol_conv with $g_{tw,mol} = g_{tw}\frac{P_a }{R_{mol} T_a}$ give similar results when plotted as functions of leaf temperature.

eq_Elmol_conv = E_lmol == g_twmol*(P_wl-P_wa)/P_a
print units_check(eq_Elmol_conv).simplify_full()
soln = solve([eq_Elmol.subs(eq_Cwl, eq_Cwl(C_wl = C_wa, P_wl = P_wa, T_l = T_a)), eq_Elmol_conv], g_twmol, E_lmol)
eq_gtwmol_gtw = soln[0][0]
print units_check(eq_gtwmol_gtw).simplify_full()
eq_gtwmol_gtw_iso = eq_gtwmol_gtw(T_l = T_a).simplify_full()
print units_check(eq_gtwmol_gtw_iso).simplify_full()

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

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

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

Model closure

Given climatic forcing as $P_a$, $T_a$, $R_s$, $P_{wa}$ and $v_w$, and leaf-specific parameters $a_s$, $a_{sH}$, $L_l$ and $g_{sw}$, we need to compute $C_{wa}$, $h_c$, $g_{bw}$ and a series of other derived variables, as described below.

The vapour concentration in the free air can be computed from vapour pressure analogously to Eq. eq_Cwl:

{eq_Cwa}

$$C_{wa} = \frac{P_{wa}}{R_{mol} T_a}$$

The heat transfer coefficient ($h_c$) for a flat plate can be determined using the non-dimensional Nusselt number ($N_{Nu_L}$):

{eq_hc}

$$h_{c} = k_a\frac{N_{Nu_L}}{L_l}$$

where $k_a$ is the thermal conductivity of the air in the boundary layer and $L_l$ is a characteristic length scale of the leaf.

For sufficiently high wind speeds, inertial forces drive the convective heat transport (forced convection) and the relevant dimensionless number is the Reynolds number ($N_{Re_L}$), which defines the balance between inertial and viscous forces \cite[Eq. 6.41]{incropera_fundamentals_2006}:

{eq_Re}

$$N_{Re_L} = \frac{v_w L_l}{\nu_a}$$

where $v_w$ is the wind velocity (m s$^{-1}$), $\nu_a$ is the kinematic viscosity of the air and $L_l$ is taken as the length of the leaf in wind direction.

In the absence of wind, buoyancy forces, driven by the density gradient between the air at the surface of the leaf and the free air dominate convective heat exchange (free'' or natural convection''). The relevant dimensionless number here is the Grashof number ($N_{Gr_L}$), which defines the balance between buoyancy and viscous forces \cite[Eqs. 9.3 and 9.65]{incropera_fundamentals_2006}:

{eq_Gr}

$$N_{Gr_L} =\frac{g(\frac{\rho_a - \rho_{al}}{\rho_{al}}) L_l^3}{\nu_a^2}$$

where $g$ is gravity, while $\rho_a$ and $\rho_{al}$ are the densities of the gas in the atmosphere and at the leaf surface respectively.

For $N_{Gr_L} \ll N_{Re_L}^2$, forced convection is dominant and free convection can be neglected, whereas for $N_{Gr_L} \gg N_{Re_L}^2$ free convection is dominant and forced convection can be neglected \cite[P. 565]{incropera_fundamentals_2006}. For simplicity, the analysis is limited to forced conditions, which is satisfied by considering wind speeds greater than 0.5~m~s$^{-1}$ for $5 \times 5$ cm leaves. %See leaf_chamber4.sws or leaf_capacitance_steady_state3.sws.

The average Nusselt number under forced convection was calculated as a function of the average Reynolds number and a critical Reynolds number ($N_{Re_c}$) that determines the onset of turbulence and depends on the level of turbulence in the free air stream or leaf surface properties \cite[P. 412]{incropera_fundamentals_2006}

{eq_Nu_forced_all}

$$N_{Nu_L} = (0.037 N_{Re_L}^{4/5} - C_1)N_{Pr}^{1/3}$$

with

{eq_C1}

$$C_1 = 0.037 C_2^{4/5} - 0.664 C_2^{1/2}$$

and

{eq_C2}

$$C_2 = \frac{N_{Re_L} + N_{Re_c} - |N_{Re_c} - N_{Re_L}|}{2}$$

Eq. eq_C2 was introduced to make Eq. eq_Nu_forced_all valid for all Reynolds numbers, and following considerations explained in our previous work \citep{schymanski_stomatal_2013}, we chose $N_{Re_c}=3000$ in the present simulations.

eq_Cwa = C_wa == P_wa/(R_mol*T_a)
print units_check(eq_Cwa)
eq_hc = h_c == k_a*Nu/L_l
print units_check(eq_hc)
eq_Re = Re == v_w*L_l/nu_a
print units_check(eq_Re)
eq_Gr = Gr == g*(rho_a - rho_al)/rho_al * L_l^3/nu_a^2
print units_check(eq_Gr)
C2 = (Re + Re_c - (abs(Re - Re_c))/2)
C1 = 0.037*C2^(4/5) - 0.664*C2^(1/2)
eq_Nu_forced_all = Nu == (0.037*Re^(4/5) - C1)*Pr^(1/3)
eq_Nu_forced_all.show()

mole/meter^3 == mole/meter^3

kilogram/(kelvin*second^3) == kilogram/(kelvin*second^3)

1 == 1

1 == 1

Thermodynamic variables

In order to simulate steady state leaf temperatures and the leaf energy balance terms using the above equations, it is necessary to calculate $\rho_a$, $D_{wa}$, $\alpha_a$, $k_a$, and $\nu_a$, while $L_l$, $Re_c$ and $g_{sv}$ are input parameters, and $P_{wa}$ and $v_w$ (vapour pressure and wind speed) are part of the environmental forcing. $D_{wa}$, $\alpha_a$, $k_a$ and $\nu_a$ were parameterised as functions of air temperature ($T_a$) only, by fitting linear curves to published data \cite[Table A.3]{Monteith_principles_2007}:

{eq_Dwa}

$$D_{wa} = (1.49\times 10^{-7})T_a - 1.96\times 10^{-5}$$

{eq_alphaa}

$$\alpha_a = (1.32\times 10^{-7})T_a - 1.73\times 10^{-5}$$

{eq_ka}

$$k_a = (6.84\times 10^{-5})T_a + 5.62\times 10^{-3}$$

{eq_nua}

$$\nu_a = (9\times 10^{-8})T_a - 1.13\times 10^{-5}$$

Assuming that air and water vapour behave like an ideal gas, and that dry air is composed of 79% N$_2$ and 21% O$_2$, we calculated air density as a function of temperature, vapour pressure and the partial pressures of the other two components using the ideal gas law:

{eq_rhoa_Pa_Ta}

$$\rho_a = \frac{n_a M_a}{V_a} = M_a \frac{P_a}{R_{mol} T_a}$$

where $n_a$ is the amount of matter (mol), $M_a$ is the molar mass (kg~mol$^{-1}$), $P_a$ the pressure, $T_a$ the temperature and $R_{mol}$ the molar universal gas constant. This equation was used for each component, i.e. water vapour, N$_2$ and O$_2$, where the partial pressures of N$_2$ and O$_2$ are calculated from atmospheric pressure minus vapour pressure, yielding:

{eq_rhoa_Pwa_Ta}

$$\rho_a = \frac{M_w P_{wa} + M_{N_2} P_{N_2} + M_{O_2} P_{O_2}}{R_{mol} T_a}$$

where $M_{N_2}$ and $M_{O_2}$ are the molar masses of nitrogen and oxygen respectively, while $P_{N_2}$ and $P_{O_2}$ are their partial pressures, calculated as:

{eq_PN2}

$$P_{N_2} = 0.79(P_a - P_{wa})$$

and

{eq_PO2}

$$P_{O_2} = 0.21(P_a - P_{wa})$$

eq_Dva = D_va == 1.49e-07*T_a - 1.96e-05
eq_alphaa = alpha_a == (1.32e-07)*T_a - 1.73e-05
eq_ka = k_a == 6.84e-05*T_a + 5.63e-03
eq_nua = nu_a == 9.e-08*T_a - 1.13e-05
eq_rhoa_Pwa_Ta = rho_a == (M_w*P_wa + M_N2*P_N2 + M_O2*P_O2)/(R_mol*T_a)
print units_check(eq_rhoa_Pwa_Ta)

kilogram/meter^3 == kilogram/meter^3

eq_Pa = P_a == P_N2 + P_O2 + P_wa
print units_check(eq_Pa)
eq_PN2_PO2 = P_N2 == 0.79/0.21*P_O2
soln = solve([eq_Pa, eq_PN2_PO2], P_O2, P_N2)
eq_PO2 = soln[0][0]
print units_check(eq_PO2)
eq_PN2 = soln[0][1]
units_check(eq_PN2)

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

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

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

eq_rhoa = eq_rhoa_Pwa_Ta.subs(eq_PN2, eq_PO2)
units_check(eq_rhoa).simplify_full()

kilogram/meter^3 == kilogram/meter^3

Example calculations

# Energy balance as a function of variables that are independent of leaf temperature
eqenbalTl = (eq_Rs_enbal - R_s).rhs().subs(eq_El, eq_Hl, eq_Rll).subs(eq_Elmol).subs(eq_Cwl).subs(eq_Pwl)
eqenbalTl.subs(cdict).args()

(C_wa, R_s, T_a, T_l, T_w, g_tw, h_c)

Below, we create a copy of cdict, which is a dictionary with general constants, generated during the definition of variables at the top of the worksheet, and some example values for external forcing and leaf properties needed to compute steady-state fluxes and leaf temperature:

# Input values
vdict = cdict.copy()
vdict[a_s] = 1.0    # one sided stomata
vdict[g_sw] = 0.01    
vdict[T_a] = 273 + 25.5
vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature
vdict[P_a] = 101325
rha = 1
vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)
vdict[L_l] = 0.03
vdict[Re_c] = 3000
vdict[R_s] = 600
vdict[v_w] = 1

Now, we compute all derived variables, using the equations described above.

# Nusselt number
vdict[nu_a] = eq_nua.rhs().subs(vdict)
vdict[Re] = eq_Re.rhs().subs(vdict)
vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)
vdict[Nu]

26.1863624980041

# h_c
vdict[k_a] = eq_ka.rhs().subs(vdict)
vdict[h_c] = eq_hc.rhs().subs(vdict)
vdict[h_c]

22.7362219510171

# gbw
vdict[D_va] = eq_Dva.rhs().subs(vdict)
vdict[alpha_a] = eq_alphaa.rhs().subs(vdict)
vdict[rho_a] =  eq_rhoa.rhs().subs(vdict)
vdict[Le] =  eq_Le.rhs().subs(vdict)
vdict[g_bw] = eq_gbw_hc.rhs().subs(vdict)
vdict[g_bw]

0.0209367439791525

# Hl, Rll
vdict[R_ll] = eq_Rll.rhs().subs(vdict)
print vdict[R_ll]
vdict[H_l] = eq_Hl.rhs().subs(vdict)
vdict[H_l]

(1.13400000000000e-7)*T_l^4 - 900.306522304087

45.4724439020342*T_l - 13573.5245047572

# El
vdict[g_tw] =  eq_gtw.rhs().subs(vdict)
vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)
vdict[P_wl] = eq_Pwl.rhs().subs(vdict)
vdict[C_wl] = eq_Cwl.rhs().subs(vdict)
vdict[E_lmol] = eq_Elmol.rhs().subs(vdict)
vdict[E_l] = eq_El.rhs().subs(eq_Elmol).subs(vdict)
vdict[E_l]

21932.0720390026*e^(-5304.00487246815/T_l + 19.4285892764401)/T_l - 386.319257853523

Now, we have a dictionary with $E_l$, $H_l$ and $R_{ll}$ as functions of leaf temperature ($T_l$) only. At steady state, eq_Rs_enbal must be satisfied, so we will substitute the above dictionary into eq_Rs_enbal, subtract $R_s$ and do a numerical search for the root of the equation to obtain steady-state leaf temperature:

# Numerical search for steady-state Tl
vdict[T_l] = find_root((eq_Rs_enbal - R_s).rhs().subs(vdict), 273, 373)
vdict[T_l]

305.6506484227355

Below, we define a function that performs all the above calculations given a dictionary with the instantaneous forcing:

def fun_SS(vdict1):
    '''
    Steady-state T_l, R_ll, H_l and E_l under forced conditions.
    Parameters are given in a dictionary (vdict) with the following entries:
    a_s, a_sh, L_l, P_a, P_wa, R_s, Re_c, T_a, g_sw, v_w
    ''' 
    vdict = vdict1.copy()

    # Nusselt number
    vdict[nu_a] = eq_nua.rhs().subs(vdict)
    vdict[Re] = eq_Re.rhs().subs(vdict)
    vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)
    
    # h_c
    vdict[k_a] = eq_ka.rhs().subs(vdict)
    vdict[h_c] = eq_hc.rhs().subs(vdict)
 
    # gbw
    vdict[D_va] = eq_Dva.rhs().subs(vdict)
    vdict[alpha_a] = eq_alphaa.rhs().subs(vdict)
    vdict[rho_a] =  eq_rhoa.rhs().subs(vdict)
    vdict[Le] =  eq_Le.rhs().subs(vdict)
    vdict[g_bw] = eq_gbw_hc.rhs().subs(vdict)   
    
    # Hl, Rll
    vdict[R_ll] = eq_Rll.rhs().subs(vdict)
    vdict[H_l] = eq_Hl.rhs().subs(vdict)   

    # El
    vdict[g_tw] =  eq_gtw.rhs().subs(vdict)
    vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)
    vdict[P_wl] = eq_Pwl.rhs().subs(vdict)
    vdict[C_wl] = eq_Cwl.rhs().subs(vdict)
    vdict[E_lmol] = eq_Elmol.rhs().subs(vdict)
    vdict[E_l] = eq_El.rhs().subs(eq_Elmol).subs(vdict)    

    # Tl
    try:
        vdict[T_l] = find_root((eq_Rs_enbal - R_s).rhs().subs(vdict), 273, 373)
    except: 
        print 'too many unknowns for finding T_l: ' + str((eq_Rs_enbal - R_s).rhs().subs(vdict).args())
    
    # Re-inserting T_l
    Tlss = vdict[T_l]
    for name1 in [C_wl, P_wl, R_ll, H_l, E_l, E_lmol]:
        vdict[name1] = vdict[name1].subs(T_l = Tlss)
    
    # Test for steady state
    if n((E_l + H_l + R_ll - R_s).subs(vdict))>1.:
        return 'error in energy balance: El + Hl + Rll - R_s = ' + str(n((E_l + H_l + R_ll - R_s).subs(vdict))) 
    return vdict

# Test
vdict = cdict.copy()
#vdict= {}
vdict[a_s] = 1.0    # one sided stomata
vdict[g_sw] = 0.01    
vdict[T_a] = 273 + 25.5
vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature
vdict[P_a] = 101325
rha = 1
vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)
vdict[L_l] = 0.03
vdict[Re_c] = 3000
vdict[R_s] = 600
vdict[v_w] = 1
resdict = fun_SS(vdict)

fun_dict_print(resdict)

C_wa            1.29441408346663
C_wl            1.91570361006325
D_va            0.0000248765000000000
E_l             185.424519010311
E_lmol          0.00420463761928142
H_l             325.157459266011
L_l             0.0300000000000000
Le              0.888469037042992
M_N2            0.0280000000000000
M_O2            0.0320000000000000
M_w             0.0180000000000000
Nu              26.1863624980041
P_a             101325         
P_wa            3212.56734153661
P_wl            4868.42309771766
Pr              0.710000000000000
R_ll            89.4180217236781
R_mol           8.31447200000000
R_s             600            
Re              1927.40122068744
Re_c            3000           
T_a             298.500000000000
T_l             305.650648423  
T_w             298.500000000000
a_s             1.00000000000000
a_sh            2              
alpha_a         0.0000221020000000000
c_pa            1010           
epsilon_l       1              
g               9.81000000000000
g_bw            0.0209367439791525
g_sw            0.0100000000000000
g_tw            0.00676759777734245
h_c             22.7362219510171
k_a             0.0260474000000000
lambda_E        2.45000000000000e6
nu_a            0.0000155650000000000
rho_a           1.16339248053449
sigm            5.67000000000000e-8
v_w             1              

Calculation using known $T_l$

If leaf temperature ($T_l$) is known (e.g. measured), but stomatal conductance $g_{sw}$ is not known, we can still calculate $H_l$ and $R_{ll}$, and then obtain $E_l$ from the energy balance equation:

eq_El_enbal = solve(eq_Rs_enbal, E_l)[0]
units_check(eq_El_enbal)

kilogram/second^3 == kilogram/second^3

# Test with known Tl but no g_sw
vdict = cdict.copy()
vdict[a_s] = 1.0    # one sided stomata
   
vdict[T_a] = 273 + 25.5
vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature
vdict[P_a] = 101325
rha = 1
vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)
vdict[L_l] = 0.03
vdict[Re_c] = 3000
vdict[R_s] = 600
vdict[v_w] = 1

vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)
vdict[nu_a] = eq_nua.rhs().subs(vdict)
vdict[Re] = eq_Re.rhs().subs(vdict)
vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)
vdict[k_a] = eq_ka.rhs().subs(vdict)
vdict[h_c] = eq_hc.rhs().subs(vdict)


vdict[T_l] = 305.65
vdict[H_l] = eq_Hl.rhs().subs(vdict)
vdict[R_ll] = eq_Rll.rhs().subs(vdict)

eq_El_enbal.subs(vdict)

E_l == 185.462402956757

Saving definitions to separate file

In the below, we save the definitions and variables to separate files in the /temp directory, one with the extension .sage, from which we can selectively load functions using %load fun_name filenam.sage and one with the extension .sobj, to be loaded elsewhere using load_session()

fun_export_ipynb('leaf_enbalance_eqs', 'temp/')
save_session('temp/leaf_enbalance_eqs')

jupyter nbconvert  --to=python 'leaf_enbalance_eqs.ipynb'
Exporting specified worksheet to .py file...
Checking if specified ipynb file was run with sage kernel...
Renaming .py file to .sage if notebook kernel was sage (to avoid exponent error)

nbconvert returned 0

Table of symbols

# 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)

Variable	Description (value)	Units
	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
	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
	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
	Gravitational acceleration (9.81)	m s
	Boundary layer conductance to water vapour	m s
	Boundary layer conductance to water vapour	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
	Characteristic length scale for convection (size of leaf)	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
	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
	Longwave radiation away from leaf	J m s
	Molar gas constant (8.314472)	J K mol
	Solar shortwave flux	J m s
	Stomatal resistance to water vapour, inverse of	s m
	Total leaf resistance to water vapour,	s m
	Density of dry air	kg m
	Density of air at the leaf surface	kg 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
	Wind velocity	m s