**Modelling Canopy Actual Transpiration in the Boreal Forest with Reduced Error Propagation**

#### **M. Rebeca Quiñonez-Piñón <sup>1</sup> and Caterina Valeo 2,\***


Received: 5 September 2020; Accepted: 22 October 2020; Published: 27 October 2020

**Abstract:** The authors have developed a scaling approach to aggregate tree sap flux with reduced error propagation in modeled estimates of actual transpiration (*Tplot*) of three boreal species. The approach covers three scales: tree point, single tree trunk, and plot scale. Throughout the development of this approach the error propagated from one scale to the next was reduced by analyzing the main sources of error and exploring how some field and lab techniques, and mathematical modeling can potentially reduce the error on measured or estimated parameters. Field measurements of tree sap flux at the tree point scale are used to obtain canopy transpiration estimates at the plot scale in combination with allometric correlations of sapwood depth (measured microscopically and scaled to plots), sapwood area, and leaf area index. We compared the final estimates to actual evapotranspiration and actual transpiration calculated with the Penman–Monteith equation, and the modified Penman–Monteith equation, respectively, at the plot scale. The scaled canopy transpiration represented a significant fraction of the forest evapotranspiration, which was always greater than 70%. To understand climate change impacts in forested areas, more accurate actual transpiration estimates are necessary. We suggest our model as a suitable approach to obtain reliable *Tplot* estimates in forested areas with low tree diversity.

**Keywords:** climate change; actual evapotranspiration; modified Penman–Monteith; sap flow; scaling methods; allometric correlations; sapwood depth; sapwood area; leaf area index

#### **1. Introduction**

Climate change is reflected in almost every ecosystem by rising temperatures and changing precipitation patterns [1]. Higher seasonal temperaturesin the boreal forest willincrease evapotranspiration rates, decrease groundwater recharge, and affect water runoff levels [2–4]. Recent findings project a decrease in boreal forest biomass due to climate change, especially in the dominant conifer species [5], and these factors will likely change the boreal forest composition [6] and therefore change overall forest evapotranspiration. In vegetated areas, reliable water balance component estimates are of great importance for water management, sustainability, wildlife conservation, and nowadays, to understand and define plant species' challenges to climate adaptation [7]. Evaporation and transpiration are two water balance components whose estimates are well known to carry uncertainty due to the use of equations and models that (1) lump the components into a single estimate of evapotranspiration (ET), and (2) focus on estimating ET under ideal atmospheric conditions and with an unlimited source of water (i.e., potential ET) [8,9]. This paper focuses on estimating actual transpiration (*Tplot*) and actual ET (*Ea*) as a first step to reduce the uncertainty introduced by assumptions of ideal conditions. Recent studies have proven that modeling actual evapotranspiration and using site-specific calibrated variables greatly reduces the error and improves estimates of actual evapotranspiration [10,11].

Scaling a single tree's point sap flow measurements to the watershed level is a complex and challenging task, not to mention the error propagated during the scaling process. Thus, most studies tend to focus on refining single components of the scaling process. For instance, direct ways to estimate actual transpiration of a single tree include the heat thermal dissipation method developed by Granier [12]. Still, sap flow radial variability should be accounted for and several researchers focus their effort on modelling and reducing the error associated with tree sap flow fluctuations [13–16].

Initial efforts to account for radial variation was to measure sap flow at different depths [17–20] of the tree trunk. The problem with measuring sap flow at different depths is that none of the thermal techniques are sensitive enough to determine the boundary between the sapwood and heartwood (i.e., sapwood depth) of a tree, and some radial flow and moisture transfer between sapwood and heartwood can be confounded with transversal sap flow. These radial variations in sap flow are also related to the species vascular structure, specifically the radial distribution and length of sapwood depth around the tree circumference [21–24]. Thus, the influence of a single tree's sapwood depth variability on estimating sap flow velocity should be acknowledged. The key point for most of the thermal techniques is that accurate measures of sapwood depth or sapwood area are required to take into account the sap flow radial pattern. Most studies measure sapwood depth using visual methods, which carry uncertainty [25], and studies suggest that more accurate methods to measure sapwood depth and sapwood area should be included in the scaling process in order to reduce error propagation at the tree scale [25–27]. The use of more accurate sapwood depth values to estimate sapwood area and tree sap mass flux should be a first step towards reduction of error propagation to larger scales (i.e., canopy).

Based on our previous work [25,26], the importance of accurate estimates of sapwood area and sapwood depth to model robust allometric correlations [28] cannot be underestimated. Allometric correlations to scale up sap flow to sap flux and to canopy transpiration normally require estimates or measurements of the structural characteristics of trees, such as diameter at the breast height, leaf area, and leaf area index. These findings indicate that there is not always a positive, direct correlation between the structural characteristics [29–38], and it is necessary to create species specific allometric correlations in the area of study.

Canopy heterogeneity is another source of uncertainty when estimating canopy and basin transpiration estimates [20,37–44]. Due to the differences in vegetation structure, each type of plant differs in its physiological process and therefore, in the amount of water required for transpiration [45]. For instance, coniferous trees require less water than deciduous trees because of their more conservative vascular structure and their tolerance to growing in xeric–mesic environments [46]. At the same time, trees' transpiration rates change according to micrometeorological conditions, solar energy, and water availability [47]. This combination of physiological, micrometeorological and energy factors generate spatial and temporal heterogeneity. Indeed, it is expected that the total transpiration of a forested area with mixed vegetation will be the aggregation of each tree's transpiration at a given time and under specific conditions.

In this work, we hypothesize that an accurate scaling approach can be reached by using the current, and most robust techniques to measure the necessary scaling parameters and in situ sap flow. We also consider that the main sources of error while scaling transpiration are (1) leaf area and sapwood area estimates, (2) sap flow radial variability, (3) canopy heterogeneity and density (interspecific and intraspecific variability), and (4) forest fragmentation. In addition to these, it is necessary to validate the final estimates, and also to estimate the propagated error.

Our previous work demonstrates that accurate estimates of canopy sapwood area [25,26,28] are predicated on obtaining sapwood depth values with a negligible error. Hence, sapwood areas of single trees can produce accurate regression models for estimating canopy sapwood area while considering the canopy's heterogeneity. This paper presents the computation of sap flow radial variations for each tree species included into the scaling approach. We use the allometric models—sapwood area and leaf area—reported in [28] to scale tree sap flow to canopy actual transpiration in a five-day ()

period. These canopy actual transpiration estimates are compared to actual evapotranspiration (*Ea*) and actual transpiration values calculated with the Penman–Monteith equation, and the modified Penman–Monteith equation, respectively.

#### **2. Methodology**

#### *2.1. Scaling Approach Concept*

Our scaling approach is based on the concept that calibrated values of tree sap flow that are aggregated to the tree scale by using allometric correlations with low uncertainty, will provide estimates of tree sap flux that can be used to compute low error estimates of canopy transpiration when plot's vegetation heterogeneity is integrated into the mathematical model. For a more accurate estimate of canopy transpiration, we included canopy heterogeneity into our model by counting the total number of individuals of each species within our 60 <sup>×</sup> 60 m<sup>2</sup> delimited plots, and we measured each tree's circumference and estimated outer bark diameter, *DOB*. Based on allometric correlations between *DOB* and sapwood depth, we computed each plot's total sapwood area. These allometric correlations are reported in [26,28]. Figures 1 and 2 explain graphically the scaling approach to estimate canopy actual transpiration.

**Figure 1.** Sap flow measured at a single tree point, can provide actual tree sap flux estimates when **Figure 1.** Sap flow measured at a single tree point, can provide actual tree sap flux estimates when the tree's sapwood area is known or accurately estimated. Each tree species sap flux can provide accurate estimates of a plot's actual transpiration by using reduced error allometric correlations of canopy structure parameters such as Leaf Area or Sapwood Area. The plot scale can provide estimates of ET and T comparable to estimates from Penman–Monteith (P-M) and Modified P-M method, respectively. Appendix A.4 details the calibration of tree sap flow measurements.

**Figure 2.** Thermal Dissipation Probes (TDPs) installed in a boreal, coniferous tree with the isolation material (upper part of the picture) ready to cover the sensors and protect them from direct exposure to solar radiation that could affect the sensors' temperature.

In summary, the expected outcomes are:


#### *2.2. Study Site*

௧ A forest region located in the Sibbald Areas of Kananaskis Valley in Alberta, Canada, was the study site for field measurements to support this work. The Kananaskis Valley is a Montane closed forest formation [48] within the Rocky Mountains [49,50]. This type of forest has ridged foothills and a marked rolling topography. The Montane forest is classified as an ecoregion within the Cordilleran eco-province and experiences unique climatic conditions arising from the combination of physiography and air masses [48]. Within the province of Alberta, the Montane forest maintains the warmest temperatures during the winter than any other forested ecosystem.

௧ ௧ Two plots of 60 × 60 m were delimited and used to scale up sap mass flow and to calculate the total rate of transpiration per plot. One plot was a pure coniferous site of *Pinus contorta* (Lodgpole pine) mixed with *Picea glauca* (White spruce), while the other was a pure deciduous site composed of *Populus tremuloides* (Trembling aspen). They are henceforth to be referred to as Coniferous site and Deciduous site, respectively. These two plots were part of the plot samples used to create allometric regression models to scale up five boreal species sapwood area to the plot level (*SAplot*) [28].

Field campaigns conducted in 2004 and 2005 were used to collect all the material and biometrics required for the entire scaling process. Field data collected at each plot included: each plot's number of trees per species, each tree's outside bark diameter at breast height (*DOB*), leaf area index (*LAIplot*), soil moisture, and sap flow velocity (*J<sup>i</sup>* ). We measured *LAIplot* within the 60 × 60 m plots using the Tracing Radiation and Architecture Canopies device (TRAC, 3rd Wave Engineering Co., Nepean, ON, Canada).

#### *2.3. Measuring Single Trees' Sap Flow*

 ை ை Sap flow in a single tree was measured using the heat dissipation technique [51]. At each plot, a group of four trees (for each species inside the plot) were set up with thermal dissipation probes (TDP-30, Dinamax, Inc., Houston, TX, USA) for periods of 48 h. The thermal dissipation probes (TDP's) were installed in the North side of the trees and covered with a special insulating material (Figure 2) to avoid direct solar incidence and overheating of the sensors that might alter the logger readings. At the same time, a set of soil moisture sensors (six sensors) was placed in the soil (below the litter) to observe the changes in soil moisture content and to later compare with the trees' water uptake. The soil moisture values are also used in the empirical calculation of *Ea*. After 48 h, another group of four trees was set up with the TDP's and soil moisture sensors. This task was performed at least four

times within each plot. We were interested in capturing the sap flow patterns of trees with different diameters (i.e., interspecific heterogeneity); therefore, the trees selected for measuring daily sap mass flow were chosen in order to cover the range of trees' outside bark diameter at breast height (*DOB*) in each plot (i.e., the largest, the smallest, the mean, and other intermediate *DOB* values of each species found inside the plot).

Sap mass flow measurements were corrected by applying the original calibration presented by Granier [52], and presented here in Appendix A. The canopy transpiration estimates were computed after sap mass flow data were corrected for radial patterns of sap flow. Trembling aspen individuals were excluded from the radial correction, since it had been proven that diffuse-porous tree radial sap flow does not vary significantly [18,53,54]. The method to compute canopy transpiration from sap mass flow data is detailed in Section 2.2.

In order to validate the scaled transpiration values, the actual forest evapotranspiration (*Ea*) and plant actual transpiration (*Tplant*) were estimated for both sites. The former estimate was computed using the Penman–Monteith equation, while the latter was computed using a modified version of the Penman–Monteith equation [47]. The mathematical theory behind the three models' computations is detailed in Sections 2.3 and 2.4.

The meteorological data required to compute *E<sup>a</sup>* and *Tplant* were collected using a HOBO meteorological station (Onset Computer Corporation, Bourne, MA, USA), which was set up in a 25 m radius clearing located inside the Barrier Lake forestry trails nearby. The installed sensors measured temperature, relative humidity, dew point, rainfall, atmospheric pressure, wind speed, gust speed, wind direction, solar radiation, and photosynthetically active radiation. The sensors were placed at height of about 3 metres above the ground level. All the variable data were collected every minute.

The sap flow values were assessed by observing the order of magnitude and their agreement between some meteorological variables and the sap flow trends. It was expected that sap flow rates would be greater in sunny, calm days than in rainy, cold, cloudy days, with a plateau at night. There are periods of the day when sap flow decreases to avoid desiccation, and some other periods in which it is known that all trees reach their maximum sap flow rates.

#### *2.4. Scaling Actual Canopy Transpiration*

#### 2.4.1. Radial Patterns of Sap Flow

The acropetal sap transport rate has a radial gradient that decreases from the outermost part of the sapwood towards the pith. Since there is enough evidence of the significance of the sap flow radial gradient while scaling up sap flux density from a single point to the entire tree [17,18,24,55,56] a sap flow radial profile function developed by [57] was used to calculate the sap flow velocity along the entire sapwood depth of each tree. The radial profile function accounts for the fractional changes in sap flow as a function of the maximum sap flow rate, the sapwood depth at which this rate occurs, the total sapwood depth and the rate at which the sap flow velocity decreases from the outer to the inner sapwood:

$$f(\mathbf{x}) = a \exp(-0.5 \left[ \frac{\mathbf{x} - \mathbf{x\_0}}{\mathcal{B}} \right]^2) \tag{1}$$

where *f*(*x*) is the sap flow rate index (expressed as a fraction), ̟ is the maximum sap flow rate (equal to 1.0) occurring at the *x<sup>o</sup>* sapwood depth, and 1/β is the rate at which the sap flow radially decreases towards the pith's trunk. In order to calculate sap flow velocity changes instead of fractional changes, Equation (1) was modified to the following form:

$$V\_{0-3}/V\_{\text{max}} = \frac{1}{\beta \sqrt{2\pi}} \int\_0^3 \exp(-0.5\left[\frac{\mathbf{x} - \mathbf{x}\_0}{\beta}\right]^2) \,d\mathbf{x} \tag{2}$$

where *V*0−<sup>3</sup> is the sap flow velocity in the first three centimeters of sapwood, and the maximum sap flow velocity is *Vmax*. Studies in variations of radial sap flow have found that in conifers, the maximum velocity or the largest portion of sap flow occurs in the first centimetre [17], the first 2 cm and 3 cm [58] of sapwood depth (from cambium to pith). Research outcomes [56] reported graphs showing that maximum sap flow occurs at 20% of the depth (from cambium to pith as well). It seems that the depth at which the maximum sap flow occurs is a standard pattern independent of the tree size. Based on these previous results, here it is assumed that *Vmax* occurs somewhere between the first two centimeters; thus, *x*0= 2 cm. Other studies have reported that the rate of decrement in radial sap flow is about 20–24% in conifers; thus, β has been assumed to equal 4 (i.e., a 25% of decrement). As *V*0−<sup>3</sup> is known; that is, it is calculated from the field measurements, *Vmax* can be estimated from Equation (2). Finally, *Vmax* is used to estimate the sap flow velocity along the entire sapwood depth (*sd*) at a specific time:

$$V\_{0-\overline{sd}} = V\_{\text{max}} \frac{1}{4\sqrt{2\pi}} \int\_0^{\overline{sd}} \exp(-0.5\left[\frac{\chi - 2}{4}\right]^2) \,d\chi \tag{3}$$

Note that *<sup>V</sup>*0−*sd* is *<sup>J</sup><sup>i</sup>* , the original symbol used by [53] to define sap flow velocity; thus,

$$F\_s = SAI\_l = SA\_{tree}V\_{0-\overline{sd\_l}} \tag{4}$$

where the sap flow velocity, *<sup>V</sup>*0−*sd<sup>i</sup>* (cm s−<sup>1</sup> ), is converted into the total trunk's mass flow, *F<sup>s</sup>* , (cm<sup>3</sup> s −1 ), by scaling from the point of measurement, to the total sapwood cross-sectional area, *SAtree* (cm<sup>2</sup> ). Values of *<sup>V</sup>*0−*sd<sup>i</sup>* were computed with Equation (3) at each time step (5 min) and then used to estimate each tree's daily *F<sup>s</sup>* .

To estimate an average canopy transpiration rate, *Tplot*, single tree transpiration values where scaled up to the whole plot. First, a diurnal average sap flow per species was estimated, and multiplied by the total sapwood area of that species inside the plot, obtaining the total average canopy water mass flow:

$$\overline{J}\_{sp} = \frac{1}{m} \sum\_{i=1}^{m} \overline{V}\_{0-\overline{sd}} \tag{5}$$

and

$$
\overline{F}\_{sp} = \overline{f}\_{sp} \, \mathcal{S} A\_{sp} \tag{6}
$$

where *Jsp* is the average sap flow velocity of the species *sp* obtained by the summation of the diurnal average sap flow velocity of each *i*th individual and divided by the total *m* individuals of the same species whose sap flow was measured. *Fsp* is the average of the species *sp* total mass flow, and *SAsp* is the total sapwood area of the species *sp*. present in a specific plot The calculation of the average canopy water mass flow (*Fplot*) is through the summation of each plot's species total water mass flow:

$$\overline{F}\_{plot} = \sum\_{i=1}^{n} \overline{F}\_{sp\_i} \tag{7}$$

The average sap flow within the plot is:

$$
\overline{J}\_{\text{plot}} = \sum\_{i=1}^{n} \overline{J}\_{sp\_i} \tag{8}
$$

To estimate *Tplot*, the *Fplot* is normally divided by a unit area of ground (i.e., 1 ha). This division allows one to observe the agreement between canopy transpiration and actual forest evapotranspiration (*Ea*), or actual forest transpiration (*Tplant*). Here, we assessed three different units of surface ground area: Plot sapwood area (*SAplot*), the actual leaf area (*LA*), and the effective leaf area (*LAe f f*). Finally, *Tplot* values were compared with an average of *E<sup>a</sup>* and *Tplant*.

#### 2.4.2. Azimuthal Sap Flow Variation

Previous work [27] where we reported that we measured sapwood depth in four different sides of the tree to account for the sapwood depth variation around the tree trunk, and therefore the azimuthal variation in sap flow. The sapwood depth was measured under the microscope, and for each tree, an average sapwood depth was computed. We consider that this is sufficient to account for azimuthal sap flow variation.

#### 2.4.3. Water Storage Capacity Estimates

It is assumed that the water stored in the tree trunk equals the amount of water replenished at night. The assumption is based on previous research focused on the contribution of a tree trunk's stored water to transpiration, under dry and wet conditions. Reported results show that on average, of the daily amount of water transpired by a tree, 14.8–20.0% corresponds to the trunk's stored water [59–62]. In Trembling aspen, the water trunk provided 11.6% of the mean daily transpiration [61] and most of the time, full replenishment for the tree trunk occurs at night [60], which creates a water balance between the tree water lost during the day and the water recharged at night. In addition, it has been determined that for scaling purposes, the error associated with water storage capacity is practically null if between individuals, the sap flux variability is low [60].

#### *2.5. Estimating Actual Forest Evapotranspiration*

The Penman–Monteith equation estimates the actual evapotranspiration of vegetated surfaces by accounting for all the micrometeorological factors that influence evapotranspiration as well as the influence of the canopy conductance and aerodynamic resistance in the rates of vegetation transpiration:

$$
\lambda E\_d = \frac{\Delta (R\_{\rm ll} - G) + \rho\_d c\_p (e^{\rho} - e\_d) / r\_d}{\Delta + \gamma \left[1 + \frac{r\_d}{r\_d}\right]} \tag{9}
$$

where λ*E<sup>a</sup>* is the latent heat of actual evapotranspiration, ∆ is the slope of the saturation vapor pressure curve (kPa◦C −1 ), *R<sup>n</sup>* is the net solar radiation, and *G* is the soil heat flux (all these terms in units of (J m−<sup>2</sup> s −1 )). The air density is denoted by ρ*<sup>a</sup>* (kg m−<sup>3</sup> ) and *c<sup>p</sup>* is the specific heat of air at constant pressure (i.e., 1010 Jkg−1◦C −1 ). The term (*e <sup>o</sup>* <sup>−</sup> *<sup>e</sup>a*) is the vapor pressure deficit (*VPD*) calculated by the difference between the saturation vapor pressure (*e o* , (kPa)) and the actual vapor pressure (*ea*, (kPa)). The psychrometric constant γ is in units of (kPa◦C −1 ). The aerodynamic terms, *r<sup>a</sup>* and *r<sup>c</sup>* are the aerodynamic resistance to vapor and heat transfer, and the bulk canopy resistance (both expressed in sm−<sup>1</sup> ). To convert the latent heat of evapotranspiration to actual evapotranspiration (*Ea*), *E<sup>a</sup>* = λ*Ea* λ in units of mms−<sup>1</sup> was employed. The equations used to solve the aerodynamic and energy parameters of the Penman–Monteith equation are detailed in Appendix A.

#### *2.6. Estimating Actual Canopy Transpiration*

Liu et al. [47] presented a modified version of the Penman–Monteith equation in order to estimate actual canopy transpiration at large scales. According to Liu et al. [47], a model such as Penman–Monteith should be adjusted by separately estimating the transpiration of shaded and sunlit leaves as follows (stratified model):

$$T\_{\text{plant}} = T\_{\text{sun}}LAI\_{\text{sun}} + T\_{\text{shade}}LAI\_{\text{shade}} \tag{10}$$

where *Tsun* and *Tshade* are the actual transpiration of sunlit and shaded leaves, respectively; *LAIsun* and *LAIshade* are the leaf area indices for sunlit and shaded leaves as well. The Penman–Monteith equation is then used by Liu et al. [47] to estimate *Tsun* and *Tshade*:

$$
\lambda T\_{\text{sun}} = \frac{\Delta(\mathcal{R}\_{\text{n,sum}}) + \rho\_a c\_p (e^{\rho} - e\_a) / r\_a}{\Delta + \gamma [1 + r\_s / r\_a]} \tag{11}
$$

and

$$
\Delta T\_{\text{slude}} = \frac{\Delta \left( R\_{\text{n,shude}} \right) + \rho\_a c\_p \left( e^o - c\_a \right) / r\_a}{\Delta + \gamma \left[ 1 + r\_s / r\_a \right]} \tag{12}
$$

where *Rn*,*sun* and *Rn*,*shade* are the net solar radiation available for sunlit and shaded leaves (Jm−<sup>2</sup> s −1 ), respectively, and *r<sup>s</sup>* is the stomatal resistance (sm−<sup>1</sup> ). The rest of the parameters and units remain the same as in Equation (9).

The boreal ecosystem productivity simulator (BEPS) provides a set of equations to calculate *Rn*,*sun* and *Rn*,*shade* (Liu et al., [47,63]). The equations compute the shortwave solar radiation for sunlit and shaded leaves as well. The net longwave solar radiation is assumed to behave equally for sunlit and shaded leaves; therefore, a single equation is used to calculate net longwave solar radiation. Thus, *Rn*,*sun* and *Rn*,*shade* are respectively given by:

$$R\_{\rm n,sum} = R\_{\rm s,sum} + R\_{\rm nl,sum} \tag{13}$$

and

$$R\_{n, \text{shape}} = R\_{\text{s,shape}} + R\_{n \text{l,s} \text{dade}} \tag{14}$$

where *Rs*,*sun* and *Rs*,*shade* are the shortwave solar radiation for sunlit and shaded leaves, respectively, and *Rnl*,*sun* and *Rnl*,*shade* are the net longwave solar radiation for sunlit and shaded leaves, respectively. The solution to these equations is formulated in Appendix A, and also provided by Liu et al. [47].

#### **3. Results**

Not all of the instrumented trees provided credible data and after all the sap flow data collection were checked for quality, only four Trembling aspen, five Lodgepole pine, and four White spruce trees provided credible sap flow measurements adequate for scaling up to the plot scale. In the case of the Coniferous site, eight days of sap flow measurements were used to calculate *Fplot* and *Tplot*. The Deciduous site provided four days of sap flow measurements and meteorological data. For the same dates, *E<sup>a</sup>* and *Tplot* were computed. Each site's daily values were averaged and compared with the average *Tplot* obtained for their respective time periods.

#### *3.1. Scaling Canopy Transpiration*

The Deciduous plot's ratio of *SAplot* to the plot's basal area was 0.57, while in the Conifer site, the ratio was 0.54 for the Lodgepole pine trees and 0.38 for the White spruce trees. Thus, the Trembling aspen showed a larger sapwood area per unit basal area at the plot scale than the conifer species. That was expected since diffuse-porous trees have larger sapwood areas in order to meet their water demand (i.e., they are less efficient at transporting water). As it is shown in the following sections, the Deciduous site drew larger mass flow per plot than the Coniferous site.

Figures 3–5 exemplify the diurnal sap flow patterns in Lodgepole pine, White spruce, and Trembling aspen. In each plot, the dashed line is *R<sup>s</sup>* and the solid line is *J<sup>i</sup>* . Two individuals of different *DOB* are presented in order to exemplify the differences in *J<sup>i</sup>* due to the tree size. Notice that the Lodgepole pine *Ji* is somewhat tempered in comparison to *R<sup>s</sup>* .

,௦௨ ,௦ௗ

,௦௨ ,௦ௗ

௧

௦,௦௨ ௦,௦ௗ

 ௧ ሜ ௧

௦

௦

,௦௨ = ௦,௦௨ + ,௦௨

,௦ௗ = ௦,௦ௗ + ,௦ௗ

ሜ

௧ ሜ

௧

௦ **Figure 3.** An example of a Lodgepole pine tree diurnal sap flow. Tree's diameter at breast height was 17 cm. Day of the year: 216, in 2004. Dashed line is *R<sup>s</sup>* and the solid line is *J<sup>i</sup>* .

௦ **Figure 4.** An example of a White spruce tree diurnal sap flow. Tree's diameter at breast height was 18 cm. Day of the year: 232, in 2004. Dashed line is *R<sup>s</sup>* and the solid line is *J<sup>i</sup>* . ௦ 

௦ **Figure 5.** An example of a Trembling aspen tree diurnal sap flow. Tree's diameter at breast height was 31 cm. Day of the year: 228, in 2004. Dashed line is *R<sup>s</sup>* and the solid line is *J<sup>i</sup>* .

௦

ሜ

ሜ

௦ ሜ

௦ ሜ

௧

௧

Each tree's sap flow pattern was analyzed in order to determine the times of initial and final daily transpiration activity. The transpiration patterns of the sampled trees showed activity starting early in the morning (around 500 and 545 h) and finishing between 1700 and 1900 h. Variations in the time at which the tree stopped transpiring and starting transpiring again were related to the meteorological changes.

The radial profile function to correct the sap flow velocity showed that the sap flow velocity values could have an underestimation of 12.5% in trees with a relatively small sapwood depth (3.5 ± 1.5 cm). The average *sd* in conifers ranged between 3.10 cm and 3.50 cm. In this particular case, if the radial profile correction could not be applied, the sap velocity will be underestimated when scaled to the entire tree. Each species, *Fsp* and *Fplot* are reported in Table 1. The Coniferous site total mass flow is the summation of the two species populating the site.


**Table 1.** *Fsp* and *Fplot* (m3day−<sup>1</sup> ) at each site (*n* is the number of individuals used per plot).

<sup>1</sup> Conifer-4 *Fplot* is the sum of *Fsp* of both species within the plot.

#### *3.2. Actual Evapotranspiration Results*

All the mathematical models to calculate *E<sup>a</sup>* are detailed in Appendix A. The most complex parameter to obtain is *rc*. A series of reduction functions were used, and the assumptions made provided half-hourly *r<sup>c</sup>* values that are in reasonable agreement with the values listed by [64–66]. The other parameter that was estimated in an uncommon way was *Rn*. This was done by integrating parameters that take into account the influence of *LAI*, gap fraction and emmisivity of understory and overstory. Since the determination of *LAI<sup>u</sup>* and Ω*<sup>u</sup>* was essentially based on previous reports, which at the same time are based on a few assumptions, it was necessary to observe the influence of *LAI<sup>u</sup>* and Ω*<sup>u</sup>* values on the calculation of *Ea*. Thus, a sensitivity analysis of *E<sup>a</sup>* was performed by using different *LAI<sup>u</sup>* and Ω*<sup>u</sup>* values. The range of values to test *LAI<sup>u</sup>* and Ω*<sup>u</sup>* were 0.6–1.5 and 0.5–0.9, respectively. The obtained estimates of *<sup>E</sup><sup>a</sup>* with respect to the initial *<sup>E</sup><sup>a</sup>* differ in the range of <sup>−</sup>2.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> to 9.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> mmd−<sup>1</sup> . When *LAI<sup>u</sup>* and Ω*<sup>u</sup>* are set up as 0.6 and 0.9, respectively, *E<sup>a</sup>* estimates are practically the same as when *LAI<sup>u</sup>* and Ω*<sup>u</sup>* are set up as 1.0 and 0.5 (the values used here), respectively. The sensitivity analysis was also performed to see the impact on the average of *E<sup>a</sup>* (i.e., *Ea*) per day. The analysis showed differences among values in the range of 2.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> to 6.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> . In conclusion, the assumed *LAI<sup>u</sup>* and Ω*<sup>u</sup>* values were considered adequate. Final estimates of *E<sup>a</sup>* are listed in Table 2. The *E<sup>a</sup>* values are shown per date and sorted by the type of site that was set up for sap flow measurements in the same dates.


**Table 2.** Penman–Monteith *E<sup>a</sup>* and *E<sup>a</sup>* estimates during the same days that sap flow was measured at each site (in 2004). *E<sup>a</sup>* is the daily *E<sup>a</sup>* average.

Liu et al. [47] reported that Canadian boreal forest evapotranspiration values range between 100 and 300 mm year−<sup>1</sup> . Additionally, Liu et al. [47] estimated that a coniferous land cover could have a yearly transpiration of 123 mm with an *s* = 55 mm; and deciduous and mixed forests land covers were reported with yearly transpiration values of 327 mm and 244 mm, respectively. On examination of the previous results, it would seem that there is an overestimation of *Ea*; however, 2004 had a particularly wet and hot summer, that exceeded reported rainfall normals [67] by a factor of 0.75 in July and 2.27 in August. Moreover, daily maximum temperatures during the months of July and August were greater than the daily maximum values reported in Environment Canada's Climatic Normals [68]. That is, July and August maximum temperatures varied between 24 and 29 ◦C, respectively, while the Climatic Normals reported maximum temperatures of 21.5 and 21.1 ◦C, respectively. Thus, the conditions for large evapotranspiration amounts that are greater than normal maximums could be considered reasonable for this wet and hot summer.

#### Variation in the Soil Field Capacity

The soils in the field plots were comprised of a sandy loam. The field capacity of a sandy loam soil varies between 0.16 and 0.22, and its wilting point is 0.073 [68]. The reported *E<sup>a</sup>* was calculated using an average value of the soil field capacity. Calculations of θ*<sup>e</sup>* , *g<sup>s</sup>* , and *E<sup>a</sup>* were made using the lower and upper bounds of the soil's θ*f c*.

Results showed that in days when θ*<sup>e</sup>* ≤ 0.00, the function limiting *E<sup>a</sup>* was *g*(θ*sm*), causing *g<sup>s</sup>* to become practically null, and making *r<sup>c</sup>* reach its maximum value. In these days, there was no difference in the final *E<sup>a</sup>* since the computation of θ*<sup>e</sup>* will always be zero or negative, no matter the θ*f c* value. Of course, in those days the factor limiting *E<sup>a</sup>* was soil moisture to the point that observed *E<sup>a</sup>* values were lower than 1 mmd−<sup>1</sup> (e.g., days 213 and 235, Coniferous site).

When 0.16 ≥ θ*sm* ≥ 0.22, soil moisture is not limiting at all, and other environmental factors drive *Ea*. In these cases, there was no variation in the final *E<sup>a</sup>* estimate. It was noticed as well that the immediate limiting factor was *VPD*, and then *R<sup>s</sup>* (e.g., days 231, Coniferous site).

Finally, if θ*sm* ≈ θ*wp*, there is variation in the estimates of *Ea*. This was noticeable for just two days in the whole data set used here (days 215 and 216, set up in the Conifer site). When θ*sm* varied from 0.0750 to 0.0795, the changes in θ*f c* generated *E<sup>a</sup>* to vary between 2.54 mmd−<sup>1</sup> and 3.73 mmd−<sup>1</sup> , when θ*f c* was set up as 0.22 and 0.16 respectively (day 215). When θ*sm* varied from 0.0735 to 0.0743, the changes in θ*f c* caused *E<sup>a</sup>* of 0.90 mmd−<sup>1</sup> , either θ*f c* was 0.22 or 0.16, respectively (day 216). The reported *E<sup>a</sup>* values for these two days are 3.01 mmd−<sup>1</sup> and 1.68 mmd−<sup>1</sup> . In those two days, it could be said that there is a variation in the *E<sup>a</sup>* estimates between 0.47 mmd−<sup>1</sup> and 0.78 mmd−<sup>1</sup> .

#### *3.3. Actual Canopy Transpiration Results*

Appendix A details the mathematical models, and the computation of *Tsun* and *Tshade* is very similar to the one applied for computing *Ea*. The main changes rely on substituting *R<sup>n</sup>* by either *Rsun* or *Rshade* and the use of *r<sup>s</sup>* instead of *rc*. Tables 3 and 4 show the obtained transpiration estimates for shaded, sunlit leaves, and the total canopy transpiration in the Conifer and Deciduous sites, respectively. It is worth mentioning that for the Deciduous site, the *Tplant* estimates were based on the estimation of *g*(VPD) computed with *KVPD* = 0.84 kPa.

**Table 3.** Modified Penman–Monteith *Tplant* estimates during the same days that sap flow was measured at the Coniferous site. *Tplant* is the summation of *Tshade* and *Tsun*. *Tplant* is the average of the daily *Tplant*. All estimates are in mm/d.


**Table 4.** Modified Penman–Monteith *Tplant* estimates during the same days that sap flow was measured at the Deciduous site (measured in 2004). *Tplant* is the summation of *Tshade* and *Tsun*. *Tplant* is the average of the daily *Tplant*. All estimates are in mm/d.


#### Variation of the Soil Field Capacity

As in the *E<sup>a</sup>* estimates, there is variation in the estimates of *Tplant* if θ*sm* ≈ θ*wp*. At the Coniferous site, days 215 and 216 showed the variations at the Coniferous site. When θ*sm* varied from 0.0750 to 0.0795, the changes in θ*f c* generated *Tplant* to vary between 2.10 mmd−<sup>1</sup> and 3.22 mmd−<sup>1</sup> , (keeping θ*f c* equal to 0.22 and 0.16 respectively; day 215). When θ*sm* varied from 0.0735 to 0.0743, the changes in θ*f c* caused *Tplant* of 0.87 mmd−<sup>1</sup> , either θ*f c* was 0.22 or 0.16, respectively (day 216). The reported *Tplant* values for these two days are 2.53 mmd−<sup>1</sup> and 1.67 mmd−<sup>1</sup> . In those two days, it could be said that there is a variation in the *Tplant* estimates between 0.43 mmd−<sup>1</sup> and 0.80 mmd−<sup>1</sup> .

#### *3.4. Agreement between Methods*

The Coniferous site's daily average estimates of *E<sup>a</sup>* and *Tplant* (Tables 2 and 3) are practically the same (1.56 mm/d). For the Deciduous site *Tplant* > *E<sup>a</sup>* by 0.13 mm/d. Hence, both Equations (9) and (10) draw very similar estimates. Such close similarity could mean that the wet, hot summer conditions of the studied area made the evaporation component negligible. Nevertheless, this should be part of future studies that could observe the agreement between the original Penman–Monteith equation and the stratified model developed by Liu et al. [47].

The comparison between *E<sup>a</sup>* and *Tplot* is shown in Table 5, while Table 6 shows the comparison between *Tplant* and *Tplot*. For these comparisons, the transpiration values are expressed as the average of the *Jsp* [mm<sup>3</sup> sapmm−<sup>2</sup> SAd −1 ] measured in the trees inside of each plot per unit ground area. This unit ground area was estimated as the ratio of *SAplot* per 1 ha (from now on referred to as "*SAplot* as unit ground area"). Additionally, *E<sup>a</sup>* and *Tplant* were averaged (i.e., *E<sup>a</sup>* and *Tplant*) on the same days for which *Fplot* was computed.

**Table 5.** *E<sup>a</sup>* and *Tplot* at the Coniferous (8 days average) and Deciduous (4 days average) sites. *SAplot* was used as the unit ground area to estimate *Tplot*. All estimates in mm/d.


**Table 6.** *Tplant* and *Tplot* at the Coniferous (8 days average) and Deciduous (4 days average) sites. *SAplot* was used as the unit ground area to estimate *Tplot*. All estimates in mm/d.


The agreement between the Coniferous *E<sup>a</sup>* and *Tplot* is acceptable and showed that *Tplot* is about 97% of the total forest evapotranspiration. The remaining 3% of *E<sup>a</sup>* may be attributed to the other sources of forest evapotranspiration such as surface evaporation and understory transpiration. The contribution of understory evapotranspiration varies and it could be fairly large during the growing season; however, [69] listed different sources that measured understory *ET* in stands of different pinaceas, and percentages range from 6% to 60% as understory contribution to forest *ET*. Thus, it is reasonable to attribute the difference between both methods to understory *ET*. Very similar results were seen in the comparison between *Tplant* and *Tplot* where *Tplot* is 97% of the *Tplant* estimates. Although both values are quite similar, the *Tplant* is greater than *Tplot* by 0.04 mm/d. The agreement is acceptable as well; however, it was expected that both values would be equal (i.e., *Tplant* = *Tplot*).

The Deciduous plot results showed a better agreement with *E<sup>a</sup>* when *KVPD* was set as 0.84 kPa and *VPD<sup>c</sup>* = 1.0 kPa. In this case, *Tplot* is about 73% of *Ea*, and about 71% of *Tplant*. This value can also be considered acceptable as well, since the days when the *Jsp* was measured, the soil moisture was not limiting, and *VPD* was the driving factor. As it has been shown in other works [52], when this situation happens, the sap flow reaches a plateau and becomes quasi constant along the day. Just as when water is limiting, *Jsp* decreases. Thus, the remnant 28% of the *E<sup>a</sup>* may be attributed to the understory transpiration and other surfaces of evaporating water.

#### *3.5. Leaf Area Indices as Scaling Factors*

Assessing other unit areas that could be helpful in transforming sap flux density values into a canopy transpiration rate, effective leaf area (*LAe f f*) and actual leaf area (*LA*) were used as unit areas as well:

$$
\overline{T}\_{\rm plot} = \overline{f}\_{\rm plot} \times \text{SAI}\_{eff} \tag{15}
$$

where

$$SAI\_{eff} = \frac{SA\_{sp}}{LAI\_{eff} \times A\_{plot}} = \frac{SA\_{sp}}{LA\_{eff}}\tag{16}$$

or

$$
\overline{T}\_{\text{plot}} = \overline{f}\_{\text{plot}} \times SAI\_{\text{actual}} \tag{17}
$$

where

$$SAI\_{actual} = \frac{SA\_{sp}}{LAI \times A\_{plot}} = \frac{SA\_{sp}}{LA} \tag{18}$$

Leaf area values are the same values as those used to create the regression model with *SAplot* in [25]. Results are shown in Tables 7 and 8. As it is appreciated, *LAplot* and *LAe f f* as unit ground areas describe the canopy transpiration of the Coniferous site as 48% and 67% of *Ea*, respectively; and the same agreements are shown with *Tplant*. In the case of the Deciduous site, *Tplot* is described as 64% and 83% of *Ea*. On the other hand, *Tplot* is 62% and 80% of *Tplant*. The *LAe f f* as a unit area describes the Deciduous *Tplot* as a larger proportion of *E<sup>a</sup>* than the unit ground area when these values are based on *KVPD* = 0.84 kPa).

**Table 7.** *LAplot*, *LAe f f* , and site average canopy transpiration over eight days at the Coniferous site.




#### **4. Discussion**

The main objectives of this study involve scaling issues in transpiration: firstly, to identify those parameters influencing transpiration at different scales in order to use them as scaling parameters if adequate models can be developed [25,26,28]; and secondly (but no less important), the improvement of the final transpiration estimates at larger scales. This is a complex task since there often exists large intra- and interspecific variability that, at the same time, is controlled by biophysical characteristics. In this study, these problems were faced and addressed by using more accurate methods to estimate the scaling factors in order to avoid large uncertainty in the final estimates. The effectiveness of using more accurate methods is proven through the validation of *Tplot* estimates. That is, *Tplot* will be reasonable result if first, *Tplot* = *Tplant*, or at least *Tplot* ≈ *Tplant*; and second, *Tplot* will be a significant proportion of *Ea*.

Each site's *Tplot* shows an acceptable agreement with the computed actual forest evapotranspiration and the actual canopy transpiration—Equations (9) and (10). In the Deciduous site case, the obtained *Tplot* motivates one to speculate if the agreement is good enough. In this particular case there is an issue worth mentioning here (in case the reader considers the *Tplot* fraction to be small). The days in which the *Tplot* was calculated, showed large *E<sup>a</sup>* and *Tplant* values because θ*sm* was not limiting, and *VPD* was driving *E<sup>a</sup>* and *Tplant* transpiration as well. In this case, the empirical factor *KVPD* was adjusted as much as possible by respecting previous reports on the influence of *VPD* on *g<sup>s</sup>* . The actual *r<sup>c</sup>* of the Trembling aspen individuals could go beyond the empirical estimates, but there is no field data that could evince this and allow modification to *KVPD*. Moreover, BEPS results suggest that at larger scales, a deciduous forest's transpiration is about 67% of the annual actual forest evapotranspiration [47]. Therefore, the Deciduous *Tplot* are considered reasonable estimates.

The three different area factors that were used to calculate actual canopy transpiration drew dissimilar results. Still, the three estimated *Tplot* values always met the expected agreements with *E<sup>a</sup>* and *Tplant*. The Coniferous site *Tplot* estimated by means of *SAplot* as a unit ground area (i.e., using *SAplot*<sup>10</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> <sup>m</sup><sup>2</sup> ) implies that there is a significant contribution of canopy transpiration to the total *ET* of the studied sites. Moreover, the estimated *Tplot* ≈ *Tplant*. Thus, the canopy transpiration rates are in good agreement with previous works when using *SAplot* as the unit ground area.

However, using any leaf area as a unit area factor, it seems that canopy transpiration is underestimated (in the Coniferous site), and overestimated (in the Deciduous site) in comparison with the obtained *Tplot* using *SAplot* as unit area. The *LAe f f* as a unit ground area that defines the Deciduous site's *Tplot* as a larger fraction of *E<sup>a</sup>* (than that estimated with *SAplot* as unit ground area). Thus, *Tplot* becomes a larger proportion of *Tplant* than using *SAplot* as a unit ground area. Conversely, *LAe f f* as a unit ground area defines the Coniferous site's *Tplot* as a smaller fraction. Results suggest that *SAplot* as unit ground area gives adequate *Tplot* estimates for the Coniferous site and the *SAIe f f* gives adequate *Tplot* estimates for the Deciduous site. Hence, the chosen unit ground area considerably influences the *Tplot* estimates.

With respect to the results from Equations (9) and (10), it was expected that *Tplant* will be a significant proportion of *E<sup>a</sup>* (i.e., between 70% and 90%). However, the daily values of *E<sup>a</sup>* and *Tplant* slightly differ; and most days showed that indeed *E<sup>a</sup>* > *Tplant* (i.e., Days 212, 215, 216, 234, 225, 226, and 228). As it was expected, *Tplant* is always a significant proportion of *E<sup>a</sup>* for most days (i.e., *Tplant* > 90% of the *Ea*. There were also days when *Tplant* > *E<sup>a</sup>* (Days 231, 232, 235, and 227). These results suggest that there were humid days causing some water condensation. Indeed, in those days, the *E<sup>a</sup>* morning estimates (i.e., *E<sup>a</sup>* hourly values) are negative. Moreover, approximately half of each day drew *VPD* ≤ 0.5 kPa.

If the volumetric soil moisture approximates its wilting point, there could be significant variations in *E<sup>a</sup>* (as well as *Tplant*). Therefore, the influence of θ*<sup>e</sup>* variations on *E<sup>a</sup>* was studied. During the days that this research was conducted, the soil was either extremely dry (below its θ*wp*) or very wet (0.16 ≥ θ*sm* ≥ 0.22), causing just two days of transition between dryness and wetness to affect *E<sup>a</sup>* values.

Even though the authors' methods help reduce the error carried out by scaling parameters, there is certainly opportunity to improve sap flow measurements. The choice to capture inter-specific variability, that is, changing the thermal probes after 48 h to individuals of different diameters, may have introduced uncertainty into the final data due to the short sampling time. The two summers in which sap flow data were collected were not the most favorable in terms of capturing sap flow patterns (such as those shown in Figures 4 and 5). Overcast days were constantly present, and we were not able to observe expected diurnal sap flow patterns (e.g., Figure 3). Another challenge was that on several occasions, the probes had to be moved to different trees because the probes were not capturing sap flow activity, which was attributed to tree infestation that was not initially visible or obvious. Thus, the sample size and temporal replication of sap flow were reduced to only a few trees and only up to four and eight days of sap flow data for the Coniferous and Deciduous site, respectively. In addition, the calibration method used to account for radial variation was estimated, and many other authors could argue that in situ calibration methods are the most effective methods, but they also indicate that this is a rare practice [15]. The available calibration methods—using potometers—are invasive, and there are concerning factors with these calibration techniques. Among those concerns is the fact

that the base of the tree is damaged by the removal of its outermost bark, which means that most likely the tree will suffer from embolism and most likely will die. Injecting pressurized water into a tree could also cause severe damage to the secondary xylem tissue, which is known to be extremely delicate and easily broken by slight changes in pressure. Amid climate change and the global biodiversity crisis, the intention with this research was to propose a scaling approach that will conserve and protect the forest trees and their environment; thus, the authors used those techniques that were the least invasive and did not require sacrifice or damage to the integrity of the trees.

#### **5. Conclusions**

The scaling approach proposed here was shown to be an appropriate way to quantify the variation of scaling factors [25,26,28] and to prove their correlation at large scales. The use of these scaling factors and the careful formulation of the scaling approach were fundamental in obtaining canopy transpiration estimates that closely agreed with the estimated actual evapotranspiration using the Penman–Monteith and the modified Penman–Monteith equations.

The canopy transpiration values calculated using the *LAe f f* as a unit ground area factor are meaningful due to the close relationship between the total amount of leaves that fully operate during transpiration. Thus, we suggest that a deeper understanding and testing of this canopy transpiration number would be a significant contribution to the study of the efficiency of trees in water use. It is also recommended that prior characterization of the intraspecific biometrics' variations be made in order to further develop the scaling approach. In addition, future work should focus on observing the behavior of this scaling approach at larger scales.

Many canopy transpiration studies disregard the impact that regression models have in the final estimates of transpiration. Previous research has demonstrated the constant over and underestimations of sapwood depth, sapwood area, leaf area, and leaf area index by using general assumptions, and the concerning level of error that these values carry to larger scales. However, the authors acknowledge that the sap flow data collected for this work would be further improved by increasing the temporal scale and sample size per species. Our suggested approach for capturing sap flow interspecific variability requires further study and thus, for future work, it is advised to measure sap flow simultaneously on a series of trees with different diameters.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4433/11/11/1158/s1. Table S1: E\_p estimates during the same days that sap flow was measured at each site. Field campaign 2004.

**Author Contributions:** Conceptualization, M.R.Q.-P. and C.V.; methodology, M.R.Q.-P.; software, M.R.Q.-P.; validation, M.R.Q.-P.; formal analysis, M.R.Q.-P.; investigation, M.R.Q.-P.; data curation, M.R.Q.-P.; writing—original draft preparation, M.R.Q.-P.; writing—review and editing, M.R.Q.-P. and C.V.; supervision, C.V.; funding acquisition, M.R.Q.-P. and C.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Alberta Ingenuity Fund, NSERC, Consejo Nacional de Ciencia y Tecnología, México, and Universidad Autónoma Metropolitana-Iztapalapa, México.

**Acknowledgments:** Thanks are due to the Biogeoscience Institute of the University of Calgary, Alberta. We would also like to thank the anonymous reviewers for their insightful comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A Computing Actual Evapotranspiration and Canopy Transpiration**

#### *Appendix A.1 Actual Evapotranspiration*

Since the direct estimation of transpiration is complex, it is more common to estimate evapotranspiration (*ET*) of forested areas as a close estimate of transpiration. For dense, homogenous vegetated areas, transpiration is usually considered the largest portion of total evapotranspiration in forested areas. In Canada, it is estimated that forest transpiration has a large proportion of the total *ET* varying between 45% and 67% of total *ET*), while the rest of the water lost is through soil evaporation or evaporation of water on surfaces (e.g., leaves and trunks) and sublimation [47]. These statements are reinforced with detailed studies of *ET* in the boreal forest that demonstrate the large activity and amounts of energy and mass fluxes [70].

In this study, the Penman–Monteith equation [71] is used to estimate the actual evapotranspiration of the vegetated areas under study. These evapotranspiration estimates will be used to validate the daily transpiration rate estimates at the plot scale. The Penman–Monteith equation estimates the actual evapotranspiration of vegetated surfaces by accounting for all the micrometeorological factors that influence evapotranspiration as well as the influence of the canopy conductance and aerodynamic resistance in the rates of vegetation transpiration:

$$
\Delta E\_{\rm a} = \frac{\Delta (R\_{\rm nl} - G) + \rho\_d c\_p (e^{\rho} - e\_a) / r\_a}{\Delta + \gamma \left[1 + \frac{r\_c}{r\_a}\right]} \tag{A1}
$$

where λ *E*<sup>a</sup> is the latent heat of evapotranspiration ∆ is the slope of the saturation vapor pressure curve (kPa ◦ C −1 ), λ *E*<sup>a</sup> is the latent heat of actual evapotranspiration ( Jkg−<sup>1</sup> ), *R<sup>n</sup>* is the net solar radiation, and *G* is the soil heat flux (all these terms in units of Jm−<sup>2</sup> s −1 ). The air density, ρ*<sup>a</sup>* is in (kgm−<sup>3</sup> ); *c<sup>p</sup>* is the specific heat of air at constant pressure (i.e., 1010 Jkg−<sup>1</sup> ◦ C −1 ). The term (*e <sup>o</sup>* <sup>−</sup> *<sup>e</sup>a*) is the vapor pressure deficit (*VPD*) calculated by the difference between the saturation vapor pressure (*e* o , (kPa)) and the actual vapor pressure (*ea*, (kPa)). The psychrometric constant, γ, is in units of (kPa ◦ C −1 ). The aerodynamic terms, *r<sup>a</sup>* and *r<sup>c</sup>* are the aerodynamic resistance to vapor and heat transfer, and the bulk canopy resistance (both expressed in sm−<sup>1</sup> ). The following paragraphs explain in detail the calculation of each Penman–Monteith equation's parameter. To convert the latent heat of evapotranspiration to actual evapotranspiration (*E*<sup>a</sup> ), use *E*<sup>a</sup> = λ *E*a/ λ in units of mms−<sup>1</sup> .

#### Appendix A.1.1 Aerodynamic Parameters

To calculate the *VPD* term in the Penman–Monteith equation, the saturation vapor pressure was initially calculated using two different equations:

$$e^o = a \circ + a\_1 T\_a + a\_2 T\_a^2 + a\_3 T\_a^3 + a\_4 T\_a^4 + a\_5 T\_a^5 + a\_6 T\_a^6 \tag{A2}$$

and

$$e^{\rho} = \exp\left(\frac{16.78 \, T\_d - 116.9}{T\_d + 237.3}\right) \tag{A3}$$

In both equations, *T<sup>a</sup>* is the air temperature (◦C, field weather station measurements). The first equation is the resultant of a Chebyshev fitting procedure used by [72]. The polynomial coefficients (i.e., *a*◦ to *a*6) are reported in Lowe's paper [72] and *e* o is calculated in mbar units. The latter equation calculates *e o* in kPa was derived by [73] and its estimates are considered of high reliability [64]. The average difference between *e <sup>o</sup>* values calculated with both equations was of 0.00017 kPa. Thus, for further estimations, Equation (A3) is applied. The actual vapor pressure is calculated using the estimated *e <sup>o</sup>* and the relative humidity (*RH*, (%)) that was measured in the field [74]:

$$e\_d = \frac{RHe^o}{100} \tag{A4}$$

The air density, ρ*a*, can be derived from [64]:

$$
\rho\_d = \frac{1000 \text{ } P}{T\_v \text{ } R} \tag{A5}
$$

where *P* is the daily mean atmospheric pressure calculated with the field measurements (barometer, units of kPa), *R* is the specific gas constant (287 Jkg−1K −1 ). *T<sup>v</sup>* is the virtual temperature in degrees Kelvin, calculated as [64]:

$$T\_v = \frac{\overline{T}\_a}{1 - (0.378 \ \overline{e}\_a P^{-1})} \tag{A6}$$

where *e<sup>a</sup>* and *T<sup>a</sup>* are taken as the daily average of *e<sup>a</sup>* and *T<sup>a</sup>* respectively. A sensitivity analysis was performed to observe how *T<sup>a</sup>* values affect ρ*<sup>a</sup>* or the evapotranspiration estimates. There were no significant changes in the values. Thus, *T<sup>a</sup>* was used in the equation. This analysis was performed since [64] did not specify if an average temperature or temperature at each hourly time-step value should be used.

The psychrometric constant can be expressed as [75]:

$$\gamma = \frac{c\_p}{\varepsilon \lambda}^{P} \tag{A7}$$

where γ is given in units of kPa◦ C −1 , *c<sup>p</sup>* is entered as 1.010 kJ kg−<sup>1</sup> ◦ C −1 , P is in kPa. The water vapor ratio molecular weight (ε) is a constant value equal to 0.622, and λ is calculated using the following equation [64]:

$$
\lambda = 2.501 - 2.361 \times 10^{-3} \, T\_a \tag{A8}
$$

where λ is given in units of MJ kg−<sup>1</sup> (i.e., multiply by 1000 to match units of *cp*).

The slope of the saturation vapor pressure curve (∆) is derived from the following equation:

$$
\Delta = \frac{4098 \, e^{\rho}}{\left(T\_a + 237.3\right)^2} \tag{A9}
$$

The aerodynamic resistance to vapor and heat flux, *ra*, is estimated with the following equation [64,76]:

$$r\_a = \left( \left[ \ln \frac{z\_{\iota} - d}{z\_{om}} \right] \left[ \ln \frac{z\_{\iota} - d}{z\_{oh}} \right] \right) \div k^2 u\_z \tag{A10}$$

where *k* is von Karman's constant (0.40), *z<sup>u</sup>* is the height (m) at which the wind speed *u<sup>z</sup>* (ms−<sup>1</sup> ) has been recorded (12.19 m in this particular case), *d* is the zero-plane displacement (m) that is assumed as 67% of the canopy height (i.e., *d* = 0.67 *hc*) for vegetation with *LAI* > 2.0. Here, the average canopy height is 15 m, which is the same height used in previous estimations. The parameters *zom* and *zoh* are the roughness lengths for the momentum and heat transfer, respectively. Allen et al. [64] suggested applying *zoh* = 0.1 *zom*. In this study, the fact that *zom* varies with cover has been taken into account; thus, *zom* is calculated differently for the Deciduous and the Coniferous sites. For the Deciduous sites, whose vegetation is considered dense and homogeneous, the equation suggested by [76] is applied:

$$z\_{om} = \frac{1}{e}(h\_c - d) = 0.37(h\_c - d)\tag{A11}$$

For the Coniferous sites, the equation suggested by [64] is applied:

$$z\_{om} = \varsigma(h\_{\mathfrak{c}} - d) \tag{A12}$$

where ς is an empirical factor that is independent of vegetation height [77]. Based on their calculated values of *zom* and *d* for conifers, [77] determined ς = 0.22. Table 1 lists the constant terms of the aerodynamic resistance equation. The ratio *zom*/*h<sup>c</sup>* = 0.7 calculated for Coniferous sites concurs with the mean value reported by [64] for this ratio. The Deciduous' sites *zom* value is between the range of values listed for deciduous trees by [64].


**Table A1.** Steady parameters in the calculation of the aerodynamic resistance to heat and vapor transfer, *ra*. All parameters are reported in meters, with exception of ς, which is unitless.

The canopy resistance is more complicated to estimate since it varies along the day and it is a function of several atmospheric parameters [78]:

$$\mathbf{g\_c} = \mathbf{g\_{cmax}}[minimum(\mathbf{g(LAI)}, \mathbf{g(R\_s)}, \mathbf{g(VPD)}, \mathbf{g(T\_4)}, \mathbf{g(\theta\_{sm})})] \tag{A13}$$

This equation implies that the canopy conductance (g<sup>c</sup> ) is a function of the environmental parameters: *LAI*, *R<sup>s</sup>* (Wm−<sup>2</sup> ), *VPD* (kPa), *T<sup>a</sup>* ( ◦C), and volumetric soil moisture (θ*sm*, in m3m−<sup>3</sup> ). The parameter that reaches its minimum at a specific time (genv), drives the canopy conductance. The lower the value of the environmental parameter reduction function, the lower the value of g<sup>c</sup> , therefore the higher the *rc*. Each parameter is represented by a reduction function that computes the value of the function between zero and one (i.e., 0 ≤ g<sup>c</sup> ≤ 1). Different authors have developed and calibrated reduction functions for calculating each one of the parameters in Equation (A13). Allen et al. [64] suggested that these equations can be replaced in the function above. Here, a set of equations was chosen and presented below. Most of the equations and empirical factors are taken from [79]; otherwise, the author is cited. Stewart [79] developed and calibrated these functions for Scots pine. This is the closest species to the species studied in this work with reported functions. In the case of the Deciduous site, the empirical factors were adjusted according to the response of *r<sup>c</sup>* or g<sup>c</sup> to the environmental parameters. This task was performed based on previous results and results that were obtained in this study.

The gcmax is the reciprocal of the minimum canopy or surface resistance (*rcmin* ). Typical values reported for coniferous forests *rcmin* range from 30 sm−<sup>1</sup> to 60 sm−<sup>1</sup> [64]. Here, an average value of the reported ranges was taken for the Coniferous site (i.e., 51 sm−<sup>1</sup> ). Ref. [80] reported maximum values of canopy conductance for Trembling aspen (31 ms−<sup>1</sup> ) and it is the one applied here for the Deciduous site. To compute g(*LAI*):

$$\log(LAI) = \frac{LAI}{LAI\_{\text{max}}} \tag{A14}$$

where *LAImax* is the maximum LAI Along the year. Since data collection occurred during the peak of the summer (July and August), it is assumed that g(*LAI*) ≈ 1.0 for both the Coniferous and the Deciduous sites. The g(*Rs*) is calculated with

$$\log(R\_s) = \frac{R\_s(1000 + K\_R)}{1000(R\_s + K\_R)}\tag{A15}$$

where *R<sup>s</sup>* is in Wm−<sup>2</sup> and *K<sup>R</sup>* is an empirical factor that was set up as 104.4 Wm−<sup>2</sup> .

The VPD function is established based on the two following equations:

$$\log(VPD) = 1 - K\_{VPD} VPD \text{ for } 0 < \text{VPD} < \text{VPD}\_{\mathbb{C}} \tag{A16}$$

and

$$\log(VPD) = 1 - K\_{VPD\_\mathbb{C}} VPD \text{ for } VPD \ge VPD\_\mathbb{C} \tag{A17}$$

with *KVPD* = 0.5 kPa. The *VPD<sup>c</sup>* is called the "threshold vapor pressure deficit" and is set up as 1.5kPa for the Coniferous site. For the Deciduous site, [52] reported the sap flow trend of four hardwood species in relation to *VPD*. One of the species studied is from the genus *Populus*. For that result, it was reported that the *Populus* sap flow did not significantly vary when *VPD* was greater than 1 kPa, unless the soil moisture content was limiting. The results presented by [52] perfectly concur with our study results. Thus, the threshold for the Deciduous site was assumed as 1 kPa. Since a *KVPD* factor was not found in the literature, its value was determined by using previously reported trends of g<sup>c</sup> versus *VPD*. Thus, the value was assumed as *KVPD* = 0.79 kPa initially. This decision was somehow conservative and based on the fact that deciduous *r<sup>c</sup>* reported values have reached 160 sm−<sup>1</sup> [64]. Therefore, *KVPD* was set up to make the reciprocal of <sup>g</sup>maxgenv to quasi match *<sup>r</sup><sup>c</sup>* to 160 sm−<sup>1</sup> when *VPD<sup>c</sup>* is greater than 1 kPa and becomes the driving environmental parameter of *rc*. Using graphs by [80] of half-hourly changes in g<sup>c</sup> and *VPD*, it was observed that *r<sup>c</sup>* can change from 81 sm−<sup>1</sup> to 200 sm−<sup>1</sup> as *VPD* reaches values greater than 1 kPa. In this case, a second run for *E*<sup>a</sup> was performed assuming *KVPD* = 0.84 kPa, to make *<sup>r</sup><sup>c</sup>* <sup>≈</sup> <sup>200</sup> sm−<sup>1</sup> when *VPD* <sup>&</sup>gt; <sup>1</sup> kPa. Values of *<sup>E</sup>*<sup>a</sup> obtained with both parameters are presented here.

For calculating g(*Ta*), a maximum and a minimum temperature (*T<sup>M</sup>* and *TN*, in ◦C) is required that constrain the stomas process, plus another empirical factor, *K<sup>T</sup>* (called the "optimum conductance temperature"):

$$\mathbf{g}(T\_d) = \frac{(T - T\_N)(T\_M - T)^\rho}{(K\_T - T\_N)(T\_M - K\_T)^\rho} \tag{A18}$$

where

$$
\rho = \frac{T\_M - K\_T}{K\_T - T\_N} \tag{A19}
$$

and *K<sup>T</sup>* is 18.35 ◦C for the Coniferous site. In the case of the Deciduous site, reported half-hour Trembling aspen g<sup>c</sup> and temperature values [80] were used to estimate the optimum conductance temperature for Trembling aspen g<sup>c</sup> . An average optimum temperature of 18.29 ◦C was obtained.

Finally, to estimate the g(θ*sm*), a function reported by [64], which is a slightly modified version of the one suggested by [79], was used:

$$\log(\theta\_{\rm sm}) = 1 - e^{-K\_{\theta}\theta\_{\varepsilon}} \tag{A20}$$

where *K*<sup>θ</sup> (*K*<sup>θ</sup> = 6.7) is the empirical factor used to calculate g(θ*sm*); and θ*<sup>e</sup>* is the fraction available for transpiration, also called the "effective fraction of available soil moisture" [64]:

$$
\Theta\_{\varepsilon} = \frac{\Theta\_{sm} - \Theta\_{wp}}{\Theta\_{fc} - \Theta\_{wp}} \tag{A21}
$$

where θ*sm* is the volumetric soil moisture (field measurements, m3m−<sup>3</sup> ), θ*wp* is the soil wilting point and θ*f c* is the soil field capacity. The values of θ*f c* and θ*wp* are obtained based on the soil texture. Direct studies of the soil type and texture in the area of Kananaskis [81–84] were used to define the soil texture in the Coniferous and Deciduous sites. The soil texture, generally defined as fine sandy loam (for both areas), drew a soil field capacity ranging between 0.16 and 0.22, while the soil wilting point was estimated as 0.07 (all values in volumetric fraction).

#### Appendix A.1.2 Energy Parameters

The soil heat flux is calculated using a "universal relationship" developed by [85]:

$$G = 0.4 \text{(}e^{-0.5LAI} \text{)} \text{R}\_{\text{ll}} \tag{A22}$$

*G* has the units of *Rn*. The net solar radiation is derived from the following equation [64]:

$$R\_{\rm ll} = (1 - \alpha)R\_{\rm s} + R\_{\rm nl} \tag{A23}$$

where *R<sup>s</sup>* is the shortwave solar radiation (measured in the field with a pyranometer), *Rnl* is the net outgoing longwave solar radiation, and α is the surface's albedo value. The term (1 − α) helps to calculate the fraction of incident net shortwave solar radiation that is absorbed by a specific surface. For coniferous forests, mean α values are in the range of 0.09–0.15 [66,76], and deciduous forests are in the range of 0.15–0.25 [76]. Monthly albedo values for mid-latitude forests are of 0.14 during the months of July and August [86–88]. The net longwave solar radiation is calculated based on the emissivity values of four different surfaces and the air temperature, *T<sup>a</sup>* [47]:

$$\begin{array}{rcl} R\_{\text{ul}} & = \left\{ \varepsilon\_{0} \left[ \varepsilon\_{a} \, \sigma\_{sb} \, T\_{a}^{4} + \varepsilon\_{u} \, \sigma\_{sb} \, T\_{a}^{4} \left( 1 - e^{-0.5 \, \text{AL}\_{\text{ul}} \, \Omega\_{\text{ul}} / \cos \, \overline{\theta}\_{\text{u}}} \right) + \varepsilon\_{\text{g}} \, \sigma\_{sb} \, T\_{a}^{4} \left( e^{-0.5 \, \text{AL}\_{\text{u}} \, \Omega\_{\text{ul}} / \cos \, \overline{\theta}\_{\text{u}}} \right) \right] \\ & & \quad - 2 \varepsilon\_{o} \, \sigma\_{sb} \, T\_{a}^{4} \left[ \left( 1 - e^{-0.5 \, \text{AL}\_{\text{b}} \, \Omega\_{\text{L}} / \cos \, \overline{\theta}\_{\text{o}}} \right) \end{array} \right\} \tag{A24}$$

where <sup>σ</sup>*sb* is the Stefan–Boltzmann constant (5.675 <sup>×</sup> <sup>10</sup>−<sup>8</sup> Jm−2K −4 s −1 ), *T<sup>a</sup>* is the air temperature (units of K). *LAI<sup>o</sup>* and *LAI<sup>u</sup>* are the Leaf Area indices of the overstory and understory respectively; Ω*<sup>E</sup>* and Ω*<sup>u</sup>* are the clumping indices of the overstory and understory; *cos* θ*<sup>u</sup>* and *cos* θ*<sup>o</sup>* are estimations of the transmission of diffuse radiant energy through the understory and overstory. The emissivity of the overstory, the ground, the understory, and the atmosphere are respectively represented by ε*o*, ε*g*, ε*u*, and ε*a*. Emissivity values for the first three surfaces are assigned from [47,89] as 0.98, 0.95 and 0.98, respectively. These emissivity values concur with values reported by [64]. Emissivity from the atmosphere is calculated with the following equation [76]:

$$
\varepsilon\_a = 1.24 \left( \frac{e\_a}{T\_a} \right)^{1/7} \tag{A25}
$$

where *e<sup>a</sup>* is in [mba] and *T<sup>a</sup>* is in degrees Kelvin. The transmission of diffuse radiant energy through the understory and overstory is given by the following two equations that were derived by [47]:

$$
\cos\bar{\theta}\_{\text{\textquotedblleft}} = 0.537 + 0.025LAI\_{\text{\textquotedblleft}}\tag{A26}
$$

$$
\cos\bar{\theta}\_0 = 0.537 + 0.025LAI\_0 \tag{A27}
$$

*LAI<sup>o</sup>* was measured for every coniferous and deciduous site (i.e., *LAIplot* = *LAIo*); *LAI<sup>u</sup>* is more complex to measure directly and it was derived from previous reports of understory NDVI and *LAI* values. Buerman et al. [90] used the reflectance values to estimate the understory NDVI and calculate *LAI* indices based on understory NDVI-*LAI* scatterplots developed by [91]. The *LAI<sup>u</sup>* values reported by [90] range between 0.6 and 1.0 (being the largest values for Black spruce and the smallest for Jack Pine). Conifers understory NDVI (NDVIu) values reported by [90] were compared with the studied Coniferous sites NDVI<sup>u</sup> calculated from the understory spectral reflectance that was recorded in the 2003 field campaign at two Coniferous and two Deciduous sites [92]. For both Coniferous and Deciduous sites, the average NDVI<sup>u</sup> is 0.8, which is 0.3 larger than the values reported by [90] in 2002 (their NDVI<sup>u</sup> range is 0.35–0.50). Using information reported by [91], ref. [90] established that an NDVI<sup>u</sup> of 0.5 corresponded to an *LAI<sup>u</sup>* of 1.0. On the other hand, ref. [92] established a standard *LAI<sup>u</sup>* value of 0.5 for broadleaf and needle-leaf forests.

Therefore, based on these previous results, *LAI<sup>u</sup>* for the Coniferous sites in Kananaskis is assumed 1.0, and for Deciduous sites, 0.6. The latter value is also in the *LAI<sup>u</sup>* range reported by [69] for deciduous stands in a boreal forest. Figure A1 is the typical understory spectral response at a Coniferous and a Deciduous site in Kananaskis Field Station. It is convenient to stress the fact that these *LAI<sup>u</sup>* values are approximate; however, the main objective is to acknowledge the importance of understory in the overall evapotranspiration estimates. Thus, as [47] thought, it is convenient to somehow include the understory evapotranspiration based on assumptions about its *LAIu*.

**Figure A1.** Typical understory spectral reflectance in Kananaskis Field Station study sites during the summer of 2003.

Ω<sup>௨</sup> The understory clumping index Ω*u*, was derived by modifying the former Chen's equation:

$$LAI = (1 - \alpha\_l)LAI\_{eff} \text{ } \mathcal{V}\_\mathcal{E} \Omega\_\mathcal{E} \tag{A28}$$

ாΩா = 1Ωா Ωா where γ*E*Ω*<sup>E</sup>* = 1Ω*<sup>E</sup>* in vascular vegetation [93]. Thus, for understory vegetation Ω*<sup>E</sup>* does not have to be partitioned into fractions that account for the shoot effect. At the same time, the α*<sup>l</sup>* value is zero since there is no fraction of wood to account for in the understory vegetation present at the study sites. Thus,

$$\text{LAI}\_{\mathfrak{u}} = \text{LAI}\_{eff} \,\Omega\_{\mathfrak{u}} \tag{A29}$$

<sup>௨</sup> <sup>௨</sup> As *LAI<sup>u</sup>* is known, *LAIe f f* can be approximated as 50% of *LAI<sup>u</sup>* as suggested by [94] for grasses (the closest that can be found to a forest understory). Hence,

$$
\Omega\_{\rm u} = \text{LAI}\_{eff} / \text{LAI}\_{\rm u} = 0.5 \text{LAI}\_{\rm u} / \text{LAI}\_{\rm u} = 0.5 \tag{A30}
$$

#### *Appendix A.2 Potential Evapotranspiration*

୮ = 0 The potential evapotranspiration results are provided in the Supplementary Material section for the reader to compare the great disparity between potential evapotranspiration and actual evapotranspiration estimated with the Penman–Monteith equation. The Penman combination equation estimates the potential evapotranspiration, or also, the free water evaporation. Potential rates of evapotranspiration assume that the water is never a limiting factor, the plant completely shades the ground (thus, there is no soil evaporation) and it has the optimal environmental conditions to transpire at its maximum rate (there is no canopy resistance). Two versions of the Penman–Monteith equation are used here to estimate the Potential Evapotranspiration (*E*p), the combination equation for free water evaporation [95,96], and the Penman–Monteith equation that includes the aerodynamic parameter but sets *r<sup>c</sup>* = 0 [97]. The former equation is computed in the following form:

$$E\_{\rm p} = \frac{\Delta (R\_n - G) + \rho\_a c\_p (e^\circ - e\_a) \mu\_2}{\lambda (\Delta + \gamma) \rho\_w} \tag{A31}$$

ଶ ∘ = ଶ ∘

3 m 2 m

<sup>ଶ</sup> <sup>ଶ</sup> = 2 m <sup>∘</sup>

where ρ*<sup>w</sup>* is the water density in units of kg m−<sup>3</sup> , and *u*<sup>2</sup> is the wind speed at 2 m height. Wind speed measured at 3 m height was scaled down to 2 m using the aerodynamic function [98]:

$$\frac{\mu\_2}{\mu\_0} = \frac{z\_2}{z\_0} \tag{A32}$$

where *u*<sup>2</sup> is the wind speed to be estimated at height *z*<sup>2</sup> = 2 m and *u*◦ is the wind speed at the reference height *z*◦ (in this case, 3 m). Wind differences of ±6 cm were registered between the two heights. The rest of the parameters were already defined. The *G* parameter is not included in the original equation; however, it was decided to slightly modify the method and include *G*. Equation (A31) gives *E*<sup>p</sup> in units of ms−<sup>1</sup> . The second one is Equation (A1), making *r<sup>c</sup>* = 0, and *E*<sup>p</sup> is given in mms−<sup>1</sup> :

$$E\_{\rm p} = \frac{\Delta (R\_{\rm n} - G) + \rho\_a c\_p (e^{\rho} - e\_a) / r\_a}{\lambda [\Delta + \gamma]} \tag{A33}$$

The obtained *E*<sup>a</sup> and *E*<sup>p</sup> daily values were averaged along the eight days (for the Coniferous site) and the four days (for the Deciduous site) and compared with the average *Tplot* value obtained for their respective period of time.

#### *Appendix A.3 Canopy Transpiration Using Modified Penman–Monteith Equation*

Liu et al. [47] used a slightly modified version of the Penman–Monteith equation in order to estimate actual canopy transpiration at large scales. According to [47], a model such as Penman–Monteith should be adjusted by separately estimating the transpiration of shaded and sunlit leaves as follows (stratified model):

$$T\_{\text{plant}} = T\_{\text{sun}}LAI\_{\text{sun}} + T\_{\text{shade}}LAI\_{\text{shade}} \tag{A34}$$

where *Tsun* and *Tshade* are the actual transpiration of sunlit and shaded leaves respectively; *LAIsun* and *LAIshade* are the leaf area indexes for sunlit and shaded leaves as well. The Penman–Monteith equation is then used by [47] to estimate *Tsun* and *Tshade*:

$$
\Delta T\_{\text{sun}} = \frac{\Delta (R\_{\text{n,sum}}) + \rho\_d c\_p (e^{\rho} - e\_a) / r\_a}{\Delta + \gamma [1 + r\_s / r\_a]} \tag{A35}
$$

and

$$
\Delta T\_{\text{slude}} = \frac{\Delta \left( R\_{\text{n, slude}} \right) + \rho\_a c\_p \left( e^o - e\_a \right) / r\_a}{\Delta + \gamma \left[ 1 + r\_s / r\_a \right]} \tag{A36}
$$

where *Rn*, *sun* and *Rn*, *shade* are the net solar radiation available for sunlit and shaded leaves in Jm−<sup>2</sup> s −1 , and *r<sup>s</sup>* is the stomatal resistance in units of sm−<sup>1</sup> . The rest of the parameters and units remain the same as in Equation (A1). The boreal ecosystem productivity simulator (BEPS) sets up a set of equations to calculate *Rn*, *sun* and *Rn*, *shade* [47,63]. The equations compute the shortwave solar radiation for sunlit and shaded leaves as well. The net longwave solar radiation is assumed to behave equally for sunlit and shaded leaves; therefore, a single equation is used to calculate net longwave solar radiation. Thus, *Rn*, *sun* and *Rn*, *shade* are respectively given by

$$R\_{\text{n, sum}} = R\_{\text{s, sum}} + R\_{\text{nl, sum}} \tag{A37}$$

and

$$R\_{\rm n, shade} = R\_{\rm s, shade} + R\_{\rm nl, shade} \tag{A38}$$

where *Rs*, *sun* and *Rs*, *shade* are the shortwave solar radiation for sunlit and shaded leaves, respectively, and *Rnl*, *sun* and *Rnl*, *shade* are the net longwave solar radiation for sunlit and shaded leaves, respectively. The shortwave solar radiation terms are calculated by the following equations:

$$R\_{\rm s,sun} = (1 - a\_{\rm L}) \left( R\_{\rm s,dir} \cos \alpha\_{\rm sa} / \cos \theta \right) + R\_{\rm s,slade} \tag{A39}$$

where α*<sup>L</sup>* is the leaf scattering coefficient (constant that equals 0.25); α*sa* is the mean leaf–sun angle, which is taken as 60◦ [47]; θ is the solar zenith angle; and *Rs*, *dir* is the direct shortwave solar radiation. *Rs*, *shade* is calculated with

$$R\_{\rm s, shade} = \left(\mathbb{R}\_{\rm s, diff} - \mathbb{R}\_{\rm s, diff-under}\right) / \text{LAI}\_o + \mathbb{C} \tag{A40}$$

where *<sup>R</sup>s*, *di f* is the diffuse shortwave solar radiation; *<sup>R</sup>s*, *di f*−*under* is the diffuse shortwave solar radiation under the overstory; and *C* accounts for the multiple scattering of direct radiation, which is calculated by

$$\mathcal{C} = a\_{\rm L} \,\Omega\_{\rm E} \,\mathcal{R}\_{\rm s, dir} (1.1 - 0.1 LAI\_0) \, e^{-\cos\theta} \tag{A41}$$

*Rs*, *dir* is a function of *R<sup>s</sup>* and *Rs*, *di f* :

$$R\_{\rm s, dir} = R\_{\rm s} - R\_{\rm s, div} \tag{A42}$$

and *Rs*, *di f* can be estimated using the following cases:

$$\frac{R\_{\text{s, diff}}}{R\_{\text{s}}} = \begin{cases} 0.13 & \text{if } \text{r} \sim 2.0\\ 0.943 + 0.734\tilde{r} - 4.9\tilde{r}^2 + 1.796\tilde{r}^3 + 2.058\tilde{r}^4 & \text{if } \text{r} \sim 0.8 \end{cases} \tag{A43}$$

where <sup>e</sup>*<sup>r</sup>* is calculated as a function of the solar constant (*SC* <sup>=</sup> 1367 Wm−<sup>2</sup> ), *R<sup>s</sup>* and θ:

$$
\widetilde{r} = \frac{R\_{\rm s}}{\rm SC \, \cos \theta} \tag{A44}
$$

and finally, *<sup>R</sup>s*, *di f*−*under* can be calculated as a function of *<sup>R</sup>s*, *di f* , <sup>Ω</sup>*E*, *LAIo*, and the angle for diffuse radiation (θ*o*):

$$R\_{\rm s, \, diff-under} = R\_{\rm s, \, diff} \left( e^{-0.5 \Omega\_{\rm E} L A l\_0 / \cos \overline{\theta}\_0} \right) \tag{A45}$$

where *cos*θ*<sup>o</sup>* is calculated using Equation (A27). The Ω*<sup>E</sup>* is of course the clumping index of the overstory, which is taken as 0.83 and 0.64 for the Coniferous and Deciduous site respectively (values obtained in situ with the TRAC optical device). As mentioned, the net longwave radiation terms are considered to behave the same for sunlit and shaded leaves. Thus, *Rnl*, *sun* = *Rnl*, *shade*, and their value is calculated by

$$R\_{\rm nl,sum} = R\_{\rm nl, shade} = \frac{R\_{\rm nl}}{LAI\_o} \tag{A46}$$

and Equation (A24) calculates *Rnl*. As it is noticed, Equations (A35) and (A36) include the term *r<sup>s</sup>* instead of *rc*. The stomatal resistance is calculated based on the *r<sup>c</sup>* values obtained with the set of reduction functions that resolve *g<sup>c</sup>* (Equations (13)–(21)) and with the *LAIo*:

$$r\_{\rm s} = LAI\_{\rm o} r\_{\rm c} \tag{A47}$$

Allen et al. [94] reported the previous equation using a *LAI* value which is standardized for crops and relatively tall grasses (i.e., 0.5 *LAI*). Here, the equation is modified to make it applicable to overstory. In addition, it is considered that shaded and sunlit leaves have similar stomatal resistances responses.

#### *Appendix A.4 Calibration of Tree Sap Flow Measurements with TDPs*

In theory, when the sap flow is constant, [52] assumed that the sap's velocity is

$$J\_i = \frac{1}{\alpha} \left[ \frac{\Delta T\_m - \Delta T}{\Delta T} \right] \tag{A48}$$

where ∆*T<sup>m</sup>* is the maximum temperature difference given when the sap flow is null (i.e., *J<sup>i</sup>* = 0), ∆*T* is the difference in temperature between the two probes at a specific time. The ratio between the temperature differences becomes the calibrated constant *K* (flux index) in Granier's 1985 technique [52]. With a sample size of 53 trees of three different species and diameters, [52] determined that the flux index has an exponential relationship with the velocity of sap flow:

$$K = 0.0206 f\_l^{-0.8124} \tag{A49}$$

with a *R* <sup>2</sup> = 0.96, and units of *J<sup>i</sup>* are 10−<sup>6</sup> ms−<sup>1</sup> . *J<sup>i</sup>* is expressed in the same way as sap flux density; that is, flow rate of sap volume per unit of sapwood area (i.e., 10−6*m*<sup>3</sup> *sap m*<sup>2</sup> *SA s* −1 ). Substituting *K* into Equation (A48) by Equation (A49), and considering the α term independent of the experimentation the sap flow velocity is estimated by

$$J\_l = 0.0119 \left[ \frac{\Delta T\_m - \Delta T}{\Delta T} \right]^{1.231} \tag{A50}$$

where *J<sup>i</sup>* is in units of cms−<sup>1</sup> . Granier [52] validated his results with the Penman equation outcomes, finding a good agreement between both set of results (of course, Penman potential evapotranspiration estimates were greater than the ones obtained with the Granier method). The studied species were Douglas fir (*Pseudotsuga menziesii*), European black pine (*Pinus nigra*), and Oak tree (*Quercus pedunculata*). In this study, *J<sup>i</sup>* radial flow is correct and azimuthal sap flow variation was corrected by measuring sapwood depth in four sides of the tree to obtain an average sapwood depth and capture sapwood depth variation around the tree trunk [25,26,28].

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
