boundary_layer_conductance.cpp

Go to the documentation of this file.

#include <cmath>                          // for std::max, std::min, pow, log
#include "../framework/constants.h"       // for celsius_to_kelvin
#include "conductance_helpers.h"          // for g_to_mass
#include "../math/roots/onedim/dekker.h"  // for dekker
#include "water_and_air_properties.h"     // for saturation_vapor_pressure
#include "boundary_layer_conductance.h"
 
double leaf_boundary_layer_conductance_campbell(
    double air_temperature,  // degrees C
    double delta_t,          // degrees C
    double lw,               // m
    double windspeed,        // m / s
    double p                 // Pa
)
{
    // Set constants
    double constexpr coef_forced = 0.147;
    double constexpr coef_free = 0.055;
 
    // Calculate conductances
    double const gbv_forced = coef_forced * sqrt(windspeed / lw);           // mol / m^2 / s
    double const gbv_free = coef_free * pow(std::abs(delta_t) / lw, 0.25);  // mol / m^2 / s
 
    // The overall conductance is the larger one
    double const gbv_leaf = std::max(gbv_forced, gbv_free);  // mol / m^2 / s
 
    // Convert to a mass conductance and return
    return g_to_mass(p, gbv_leaf, air_temperature + delta_t);  // m / s
}
 
double leaf_boundary_layer_conductance_nikolov(
    double air_temperature,  // degrees C
    double delta_t,          // degrees C
    double ea,               // Pa
    double gsv,              // m / s
    double lw,               // m
    double windspeed,        // m / s
    double p                 // Pa
)
{
    double const leaftemp = air_temperature + delta_t;                             // degrees C
    double const Tak = air_temperature + conversion_constants::celsius_to_kelvin;  // K
    double const Tlk = leaftemp + conversion_constants::celsius_to_kelvin;         // K
 
    double const esTl = saturation_vapor_pressure(leaftemp);  // Pa
 
    // Set constants
    double constexpr cf = 4.322e-3;   // for broad leaves
    double constexpr ce = 1.6361e-3;  // for broad leaves
    double constexpr ct = 0.378;
 
    // Equation 29
    double const gbv_forced = cf * pow(Tak, 0.56) *
                              sqrt((Tak + 120) * windspeed / (lw * p));  // m / s
 
    // This lambda function equals zero only if gbv_free satisfies the Nikolov
    // model equations for free boundary layer conductance. Here, gbv_free
    // should be expressed in m / s.
    auto check_leaf_gbv_free = [=](double const gbv_free) {
        // Equation 35
        double const eb = (gsv * esTl + gbv_free * ea) / (gsv + gbv_free);  // Pa
 
        // Equation 34
        double const Tvdiff = Tlk / (1.0 - ct * eb / p) -
                              Tak / (1.0 - ct * ea / p);  // K
 
        // Equation 33
        double const new_gbv_free =
            ce * pow(Tlk, 0.56) * sqrt((Tlk + 120) / p) *
            pow(std::abs(Tvdiff) / lw, 0.25);  // m / s
 
        return gbv_free - new_gbv_free;  // m / s
    };
 
    // Run the Dekker method; check_leaf_gbv_free is always positive for
    // gbv_free = 0, but it is difficult to find a finite value where
    // check_leaf_gbv_free is guaranteed to be negative; here we just use a
    // very large value and hope for the best.
    root_finding::dekker solver(500, 1e-12, 1e-12);
 
    root_finding::result_t result = solver.solve(
        check_leaf_gbv_free,
        1e-4,  // first guess
        0,     // lower bound of initial bracket
        0.5    // upper bound of initial bracket
    );
 
    // Throw exception if not converged
    if (!root_finding::is_successful(result.flag)) {
        throw std::runtime_error(
            "gbv_free solver reports failed convergence with termination flag:\n    " +
            root_finding::flag_message(result.flag));
    }
 
    // Get final value
    double const gbv_free = result.root;  // m / s
 
    // The overall conductance is the larger one
    return std::max(gbv_forced, gbv_free);  // m / s
}
 
double canopy_boundary_layer_conductance_thornley(
    double CanopyHeight,    // m
    double WindSpeed,       // m / s
    double minimum_gbw,     // m / s
    double WindSpeedHeight  // m
)
{
    // Define constants used in the model. `kappa` is von Karmon's constant.
    // In the original text, the value of `dCoef` is reported as 0.64. In the
    // 2000 reprinting of the text, the authors state that this value should be
    // 0.77. See "Errata to the 2000 printing" on the page after the preface of
    // the 2000 reprinting of the 1990 text.
    constexpr double kappa = 0.41;      // dimensionless, Thornley and Johnson pgs 414 and 416.
    constexpr double ZetaCoef = 0.026;  // dimensionless, Thornley and Johnson 1990, Eq. 14.9o
    constexpr double ZetaMCoef = 0.13;  // dimensionless, Thornley and Johnson 1990, Eq. 14.9o
    constexpr double dCoef = 0.77;      // dimensionless, Thornley and Johnson 1990, Eq. 14.9o
 
    // Apply the height limit
    CanopyHeight = std::min(CanopyHeight, 0.98 * WindSpeedHeight);  // meters
 
    // Calculate terms that depend on the canopy height (Eq. 14.9o)
    const double Zeta = ZetaCoef * CanopyHeight;    // meters
    const double Zetam = ZetaMCoef * CanopyHeight;  // meters
    const double d = dCoef * CanopyHeight;          // meters
 
    // Calculate the boundary layer conductance `ga` according to Thornley and
    // Johnson Eq. 14.9n, pg. 416
    const double ga0 = pow(kappa, 2) * WindSpeed;                   // m / s
    const double ga1 = log((WindSpeedHeight + Zeta - d) / Zeta);    // dimensionless
    const double ga2 = log((WindSpeedHeight + Zetam - d) / Zetam);  // dimensionless
    const double gbv = ga0 / (ga1 * ga2);                           // m / s
 
    // Apply the minimum
    return std::max(gbv, minimum_gbw);  // m / s
}