Estimation of the Latent Heat Flux over Irrigated Short Fescue Grass for Different Fetches

: Often in agrometeorology the instrumentation required to estimate turbulent surface ﬂuxes must be installed at sites where fetch is not sufﬁcient for a sector of wind directions. For different integrated ﬂux-footprints (IFFP) thresholds and taking as a reference the half-hourly latent heat ﬂuxes (LE) measured with a large weighing lysimeter (LE Lys ), the eddy covariance (EC) method and two methods based on surface renewal (SR) analysis to estimate LE were tested over short fescue grass. One method combined SR with the ﬂux-gradient (proﬁle) relationship, SR-P method, and the other with the dissipation method, SR-D method. When LE was estimated using traces of air moisture, good performances were obtained using the EC and the SR-P methods for samples with IFFP higher than 85%. However, the closest LE estimates were obtained using the residual method. For IFFP higher than 50%, the residual method combined with the sensible heat ﬂux estimates determined using the SR-P method performed close to LE Lys and using the SR-D method good estimates were obtained for accumulated LE Lys . To estimate the sensible heat ﬂux, the SR-D method can be recommended for day-to-day use by farmers because it is friendly and affordable.


Introduction
Evapotranspiration, or latent heat flux (LE), is required as input to improve weather forecasting, fire index, irrigation, dynamics of ecosystems and many other applications. In agriculture, most of the complexity to estimate LE is greatly reduced in lands where homogeneous crops grow in an open and flat terrain scenario [1][2][3]. Two crucial simplifications to estimate LE for such scenario are the following. First, the simplified surface energy balance equation holds. That is, the net radiation, R n , minus ground heat flux, G, often is close to the sum of sensible heat, H, and LE; (R n − G) = (H + LE). This equation is involved in different methods to estimate LE, such as in the Penman-Monteith equation, the Bowen ratio-surface energy balance method and the residual method [1,[4][5][6][7][8][9][10]. Second, measurements can be taken in the inertial sublayer (i.e., it is not required to install tall masts) where the Monin-Obukhov Similarity Theory (MOST) relationships hold. As a rule of thumb, the basis of the inertial sublayer is estimated as at least two times the canopy height and crucial canopy parameters involved in MOST relationships, such as the zero-plane displacement and the aerodynamic roughness lengths for momentum and heat, may be estimated as a proportion of the canopy height [1,2,11,12].
Unfortunately, the latter scenario is rarely met. Often, crops grow in small plots which may be surrounded by different crops, drier fields or bare soils. From a micrometeorological point of view, small plot refers to a surface which dimensions may compromise the method used to estimate LE. Measurements must be taken within the internal surface boundary layer, and the gold rule used as a guide in the field is to deploy the instrumentation at where k is the von Kármán constant, z is the measurement height (Z) above the zero-plane displacement (d), ζ is the stability parameter ζ = z/L 0 where L 0 is the Obukhov length defined as L 0 = − u 2 * kgT * T (T and g are the temperature of the air (virtual) expressed in Kelvin and acceleration due to gravity, respectively), ∅ h (ε) is the flux-gradient stability function for heat transfer [33] , and A s and τ s are the amplitude and period, respectively, of the mean ramp dimensions observed in the halfhourly trace of the scalar (typically, measured at 10-20 Hz). For traces of the temperature of the air, the ramp amplitude and ε has opposite sign (see Appendix A). The friction velocity may be estimated through the wind log-law which expresses the friction velocity (u * WP where the index WP denotes wind profile) as [1] u * WP = k u where u is the mean wind speed, z 0m is the roughness length for momentum and Ψ is the integrated Businger-Dyer relationship for momentum to correct for buoyancy, Ψ = ln (0.5(1 + x)) 2 + ln 0.5 1 + x 2 − 2 tan −1 x + 0.5π ξ ≤ 0 −5ξ ξ > 0 being x = (1 − 16ξ) 1/4

The SR-D Method
For steady and horizontally homogeneous flow, the budget equation for 0.5s 2 with negligible local source (sink) is [11] − w s ∂s ∂z The first term on the left-hand side of Equation (3) represents the production rate of 0.5s 2 associated with turbulent motions occurring within a mean scalar gradient. The second term represents the mean turbulent transport and s , namely mean dissipation rate, represents molecular dissipation of the scalar variance. For simplicity, the transport was neglected [1,28,34]; however, it can not necessarily be small either for stable [35,36] and unstable cases [37][38][39][40][41][42][43]. Therefore, given that it is difficult to relate w s and s through Equation (3), an original dimensional analysis involving z, the turbulent standard deviation of s, σ s , w s and s was proposed to relate w s and (z s /σ s ). The following dissipation, D, method was proposed [42], is the normalized variance similarity relationship [11]) that performed fairly constant regardless of the stability conditions, k D~1 .66 [43], and δ (ε) is the normalized production: dissipation rates ratio, δ (ε) = where ϕ (ε) is the MOST-based relationship that normalizes s . For temperature of the air, the production and dissipation terms in Equation (3) were compared, and it was shown that δ (ε) can be expressed as [22,23]. Assuming surface-layer similarity for the parameter k D and δ (ε) , the latter relationship allows rewriting the D method proposed in [42,43] in a generalized form In Equation (4) the sign of the kinetic surface flux is given by the sign of the ramp amplitude.

The Site and Climate
The experiment was carried out from 15 July to 8 September 2009 in the Alameda del Obispo experimental station (37 • 51 N, 4 • 51 W, 110 m a.s.l) which is located at the IFAPA Agricultural Training and Research Center in Córdoba (Spain). The landscape is flat and agriculture is the main activity. The plot was surrounded by short irrigated crops and bare soils corresponding to other experiments ( Figure 1). resents molecular dissipation of the scalar variance. For simplicity, the transport was neglected [1,28,34]; however, it can not necessarily be small either for stable [35,36] and unstable cases [37][38][39][40][41][42][43]. Therefore, given that it is difficult to relate and through Equation (3), an original dimensional analysis involving z, the turbulent standard deviation of s, σs, and was proposed to relate and (z /σs). The following dissipation, D, method was proposed [42], * * = k ( ) , where kD is a similarity based parameter, is the normalized variance similarity relationship [11]) that performed fairly constant regardless of the stability conditions, kD ~1.66 [43], and ( ) is the normalized production: dissipation rates ratio, ( ) where ( ) is the MOSTbased relationship that normalizes . For temperature of the air, the production and dissipation terms in Equation (3) were compared, and it was shown that ( ) can be expressed as ( ) = ( ) [22,23]. Assuming surface-layer similarity for the parameter kD and ( ) , the latter relationship allows rewriting the D method proposed in [42,43] in a generalized form * * = 1.66 (4) In Equation (4) the sign of the kinetic surface flux is given by the sign of the ramp amplitude.

The Site and Climate
The experiment was carried out from 15 July to 8 September 2009 in the Alameda del Obispo experimental station (37°51′ N, 4°51′ W, 110 m a.s.l) which is located at the IFAPA Agricultural Training and Research Center in Córdoba (Spain). The landscape is flat and agriculture is the main activity. The plot was surrounded by short irrigated crops and bare soils corresponding to other experiments ( Figure 1). In summer, the climate is characterized by the stagnation of the Azores high pressure. Regional advection of sensible heat flux is often severe along the Guadalquivir Valley [44]. The latter understood as that unstable cases are typically observed from dawn to about 14.00 h (local time). Afterwards (i.e., around 2-3 h after noon), positive LE exceeding the available net surface energy is observed until dusk [45,46]. For this case, during the experiment the stability parameter ranged between 0 and 0.22. The prevailing wind direction is In summer, the climate is characterized by the stagnation of the Azores high pressure. Regional advection of sensible heat flux is often severe along the Guadalquivir Valley [44]. The latter understood as that unstable cases are typically observed from dawn to about 14.00 h (local time). Afterwards (i.e., around 2-3 h after noon), positive LE exceeding the available net surface energy is observed until dusk [45,46]. For this case, during the experiment the stability parameter ranged between 0 and 0.22. The prevailing wind direction is west. The weather was typical of the climate which is characterized by clear skies, high temperatures, high evapotranspiration rates and a few short convective rainy events. One short rainy event (about 10 min) was observed during the experiment, and the Table 1 shows the mean maximum and minimum daily values of the maximum and minimum temperatures and the relative humidity of the air, the wind speed, the solar radiation and the evapotranspiration measured with a weighing lysimeter (described next) observed during the experiment. Tx and Tn, maximum and minimum temperature of the air, respectively; HRx and HRn, maximum and minimum relative humidity of the air, respectively; U, wind speed; Rs, solar radiation; ET, evapotranspiration.

Experimental Set Up
A weighing lysimeter was located at the center of a 1.3 ha (110 × 120 m 2 ) rectangular plot of well-watered (sprinkler irrigated) tall-fescue grass (Festuca arundinacea) clipped at about 0.10-0.15 m [45]. The grass was maintained to preserve homogeneity and a fairly constant canopy height. The plot was oriented in the prevailing wind direction and the edge-lysimeter distance was 81 m. The weighing lysimeter consisted of a steel box (3 m × 2 m × 1.5 m deep) measuring half-hourly weight changes with an accuracy of 0.2 kg. Thus, the measurement error is slightly smaller than 50 Wm −2 . Soil moisture in and out of the lysimeter were measured at different points using a FDR Diviner probe (Sentek). No trends or spikes were observed along the campaign.
The net radiation was measured half-hourly using a four components net radiometer, CNR1 (Kipp & Zonen). It was deployed adjacent to the lysimeter at 1.5 m height. Half-hourly G was determined using two HFP01 soil heat flux plates (Hukseflux Thermal Sensors, Delft, The Netherlands) and four TCAV soil temperature averaging sensor (Campbell Scientific, Inc., Logan, UT, USA). Below the radiometer, the heat flux plates (placed at 0.08 m depth and separated 1.5 m) were sampled every 10 s, and the mean reading was recorded half hourly. Thermocouples were placed at 0.06 and 0.02 m depth to obtain the half-hourly mean soil temperature above the plates. The soil bulk density and the specific heat capacity was 1.3 g m −3 and 900 J kg −1 K −1 , respectively, and the soil heat flux was estimated summing the measured heat flux at 0.08 m and the ground heat storage above the soil heat flux plate [47,48].
A 3-D sonic anemometer and a krypton hygrometer (CSAT3 and KH20, respectively, Campbell Scientific, Inc., Logan, UT, USA) and a fine wire thermocouple (FW05 Campbell Scientific, Inc., Logan, UT, USA) were deployed at 0.9 m and 145 m in the prevailing wind direction. The minimum fetch (North) was 17 m. The distance between the sonic anemometer and the hygrometer measuring paths was 0.11 m. The three wind speed components, sonic and air temperature and humidity fluctuations were measured at 20 Hz and recorded in a memory card using a CR5000 (Campbell Scientific, Inc., Logan, UT, USA) datalogger. The analytical flux-foot print, FFP, model [13] was taken as a guide for instrumentation positioning. The FFP for unstable cases peaked at distances (mean streamwise) between the EC system and the lysimeter. For stable cases, it peaked at distances closer to the lysimeter. Regardless of the stability case, the FFP peaked within the field of view of the outgoing longwave radiation emitted by the surface. Other standard measurements were the half-hourly mean wind speed and direction at 2.0 m (wind monitor RM Young 05103), humidity and temperature of the air at 1.5 m (HMP45C probe, Vaisala, Vantaa, Finland) and precipitation (tipping bucket rain gauge ARG100 deployed at 0.8 m).

Flux Estimates
The measured sonic temperature, which is close to the virtual temperature [5,15], the temperature of the air and the relative humidity were used to estimate the water vapor pressure, the atmospheric pressure, the density of the air and the specific humidity (q). The latter was used to determine the isobaric specific heat capacity of the air, C p = (1 + 0.84q) C pd where C pd is the isobaric specific heat capacity of dry air (C pd = 1005 J/(kg K)). The heat of vaporization was calculated as L= −0.0000614342 T 3 + 0.00158927 T 2 − 2.36418 T + 2500.79 kJ/(kg) where T is the temperature of the air expressed in degrees Celsius [1,5,15]. The von Kármán constant was set to 0.4 [11,28]. The zero-plane displacement and the roughness length for momentum were estimated as a portion of the canopy height, h, d = 0.7h and z o = 0.12h [1]. The ramp dimensions were determined according to [56] as described in [31] (see Appendix A). The sensible heat flux using the SR-P method was determined iterating Equations (1) and (2) because the stability parameter is required as input. Iterations started setting neutral conditions and the criterion for convergence was that the change in friction velocity (i.e., between two consecutive iterations) was smaller than 0.01 m/s [1,22,57]. Once the sensible heat flux was determined, the stability parameter obtained was used as input to determine the first proxy for the latent heat flux, LE 0 , which was corrected for density fluctuations [50,55] For the SR-D method, once the sensible heat flux was determined, the latent heat flux obtained using Equation (4) was corrected for density fluctuations. Hence, while the mean wind speed and the trace of the temperature of the air are required as input to estimate the sensible heat flux, in addition, traces of air moisture are required to estimate the latent heat flux.

Procedure for Comparison
The LE determined using the EC method, LE EC , the SR-P method, LE SR-P , the SR-D method, LE SR-D , and the residual method using different sensible heat flux estimates, (R n − G − H EC ), (R n − G − H SR-P ) and (R n − G − H SR-D ), were compared against the measured using the lysimeter, LE Lys (reference). The performance was evaluated using linear regression analysis (slope, intercept and coefficient of determination, R 2 ), the root mean square error, RMSE, and the ratio, RD, defined as, RD = ∑ y ∑ x where variables x and y corresponds to the reference and estimated values, respectively. The coefficient RD is an accumulated ratio used to determine the percentage (p) being over-or underestimated (p = 100(1 − RD)) which also gives an integrated evaluation of the bias by averaging out random errors in the half-hourly estimates (i.e., the bias is (RD − 1) times the mean reference value determined from the observations). Given that for well-watered crops the net radiation and the latent heat flux dominate the surface energy balance, the Figure 2 suggests that LE Lys was more reliable than LE EC . Figure 2 shows the closures (Rn − G − HEC − LELys) and (Rn − G − HEC − LEEC) obtained for all the samples with IFTT ≥ 85% along the experiment. In general, (Rn − G − HEC − LELys) closed the surface energy balance better than (Rn − G − HEC − LEEC). Given that for wellwatered crops the net radiation and the latent heat flux dominate the surface energy balance, the Figure 2 suggests that LELys was more reliable than LEEC. The Table 2 shows the number of samples, N, the slope, a, intercept, b, and the coefficient R 2 of the linear regression analysis, the RMSE and the coefficient RD to compare the latent heat flux estimates versus LELys for each dataset and for all the data. Regardless of the IFFP threshold, Table 2 shows that the RD values obtained using LEEC were the worst. In principle, the high coefficients of determination obtained for LEEC are not surprising because the EC method directly measures the turbulence (i.e., EC fluxes are expected highly correlated with the actual eddy fluxes). Therefore, regardless of the underestimation (RD values were smaller than one), LEEC and LELys are expected highly correlated. The Table 2 shows that LESR-P was highly correlated with LELys; in fact, LESR-P did the highest coefficient of determination for the 50% ≤ IFFP ≤ 75% dataset. The slopes of the linear regression analysis obtained using LESR-P were close to the coefficient RD because the intercepts were small and, though LESR-P underestimated the reference, it consistently gave the lowest RMSE. Therefore, LESR-P did the best performance in estimating the latent heat flux in half-hourly basis.  The Table 2 shows the number of samples, N, the slope, a, intercept, b, and the coefficient R 2 of the linear regression analysis, the RMSE and the coefficient RD to compare the latent heat flux estimates versus LE Lys for each dataset and for all the data. Regardless of the IFFP threshold, Table 2 shows that the RD values obtained using LE EC were the worst. In principle, the high coefficients of determination obtained for LE EC are not surprising because the EC method directly measures the turbulence (i.e., EC fluxes are expected highly correlated with the actual eddy fluxes). Therefore, regardless of the underestimation (RD values were smaller than one), LE EC and LE Lys are expected highly correlated. The Table 2 shows that LE SR-P was highly correlated with LE Lys ; in fact, LE SR-P did the highest coefficient of determination for the 50% ≤ IFFP ≤ 75% dataset. The slopes of the linear regression analysis obtained using LE SR-P were close to the coefficient RD because the intercepts were small and, though LE SR-P underestimated the reference, it consistently gave the lowest RMSE. Therefore, LE SR-P did the best performance in estimating the latent heat flux in half-hourly basis.

Results
Regardless of the IFFP threshold, Table 2 shows that LE SR-D had the lowest correlations, the highest RMSE and the closest RD values to one. This suggested that LE SR-D is useful to estimate half-hourly accumulations of latent heat flux.
The disagreement between the LE estimates (LE EC , LE SR-P and LE SR-D ) with respect to LE Lys increased with the IFFP threshold (Table 2). Comparing the 85% ≤ IFFP and the 50% ≤ IFFP ≤ 75% datasets, the impact in the RD values was of about 30%. It highlights the importance to take measurements within the internal boundary layer. In our case, it reflected the contribution of the drier surroundings to the traces of air moisture.
The impact of insufficient fetch was minimized when the latent heat flux was estimated using the residual method. Though some improvement with respect to LE EC , LE SR-P and LE SR-D could be expected because the latter require traces of air moisture as input and because (R n − G) mostly weights to LE for watered surfaces, the close agreement between the estimates and the reference ( Table 2) was surprising. The performance for the IFFP < 85% datasets was similar to the IFFP ≥ 85% dataset. In particular, the performance shown in Table 2 for (R n − G − H EC ) and (R n − G − H SR-P ) is similar to the obtained in other experiments over grass with sufficient fetch [58]. Figure 3 compares (R n − G − H SR-P ) against LE Lys for all the data. It shows that (R n − G − H SR-P ) scattered the full range of LE Lys and Table 2 shows that the RMSE was similar to the measurement error of LE Lys (about 50 Wm −2 ), regardless of the IFFP threshold.   Table 2 shows that (Rn − G − HEC) performed slightly closer to LELys than (Rn − G − HSR-P), but in practice such improvement may not justify the use of the EC method. That is, in terms of cost, the difference between the EC and the SR-P method is significant. In relation to (Rn − G − HSR-D), it consistently gave the closest RD values to one, but the worst R 2 and RMSE values. Table 2 shows that (Rn − G − HSR-D) performed similarly to (Rn − G − HEC) and (Rn − G − HSR-P) for the 50% ≤ IFFP ≤ 75% dataset but not as close for the other datasets. The latter was as a consequence that most samples in the 50% ≤ IFFP ≤ 75% dataset were collected near neutral conditions. In any case, (Rn − G − HSR-D) consistently gave RD values close to one; and therefore, it can be recommended to estimate accumulations of latent  Table 2 shows that (R n − G − H EC ) performed slightly closer to LE Lys than (R n − G − H SR-P ), but in practice such improvement may not justify the use of the EC method. That is, in terms of cost, the difference between the EC and the SR-P method is significant. In Atmosphere 2021, 12, 322 9 of 13 relation to (R n − G − H SR-D ), it consistently gave the closest RD values to one, but the worst R 2 and RMSE values. Table 2 shows that (R n − G − H SR-D ) performed similarly to (R n − G − H EC ) and (R n − G − H SR-P ) for the 50% ≤ IFFP ≤ 75% dataset but not as close for the other datasets. The latter was as a consequence that most samples in the 50% ≤ IFFP ≤ 75% dataset were collected near neutral conditions. In any case, (R n − G − H SR-D ) consistently gave RD values close to one; and therefore, it can be recommended to estimate accumulations of latent heat flux.

Discussion
Both the flux-gradient relationship and the budget for the turbulent variance of a scalar are principles of physics. Therefore, under the same simplifications (i.e., steady and horizontally homogeneous flow), it must be explained why the SR-P method performed closer to LE Lys than using the SR-D method. Partly, it can be explained as a consequence that MOST relationships are semi-empirical, being ∅ (ε) more robust to dissimilarity than ∅ σ(ε) [59]. The latter implying that the parameter k s may not be expected fairly constant.
It was surprising that LE SR-P performed closer to LE Lys than LE EC and that the best closures of the surface energy balance were obtained using the sensible heat flux determined with the SR-P method. Though it is difficult to explain, for rangeland grass and cotton fields with sufficient fetch, the available net surface energy was slightly better explained using the sum of the sensible heat and the latent heat fluxes determined with the SR-P method than the EC method [57,60]. This explanation is a matter of pending research. Both methods are based on different grounds; the EC method is a simplification of the Navier-Stokes equations, and the SR method is a solution of the mass conservation. Hence, inputs and calculations required are different for each method. In any case, it was shown that the SR-P method has potential for short irrigated crops growing in small plots because it can operate close to the ground, which helps to better match the flux-footprints [14]. Based on MOST, given that similarity relationships tend to become independent of the stability parameter when measurements are taken at low heights, perhaps LE SR-D could perform closer to LE Lys if the instrumentation was deployed at lower heights. Likely, it deserves further research because (according to author's knowledge) this is the first application of LE SR-D .
The combination of the residual method with SR-P sensible heat flux estimates is convenient. The latter not only because the instrumentation required is affordable. Here, the results obtained allow stating that this method can be recommended for irrigated small plots and prior studies have shown excellent performance (i.e., (R n − G − H SR-P ) was close to the LE determined using the EC method and weighing lysimeters) over crops with sufficient fetch [22,23,27,58,60]. It follows that (R n − G − H SR-P ) appears useful for precision agriculture and may allow validating methods or strategies for irrigation. In addition, larger fetches are required to estimate LE when the surface energy balance is forced using the Bowen ratio determined from gradients of temperature and moisture of the air [14,61]. Likely, a two-dimensional sonic anemometer, provided that it allows recording the sonic temperature at high frequency, is the most convenient instrument to estimate H SR-P because it is robust and measures the wind speed with high precision. It allows estimating the ramp dimensions and the friction velocity using different methods [62,63]. When half-hourly latent heat flux accumulations are of interest, (R n − G − H SR-D ) appeared friendly because H SR-D performed as a method exempt of calibration (the errors nearly balanced the accumulated estimates of latent heat flux) and, currently, free-software is available to estimate the ramp dimensions [64]. According to [27,65], likely (R n − G − H SR-D ) is the most affordable and friendly method to estimate latent heat fluxes on a time basis such as daily and weekly. It is worth mentioning that research on the SR-D method to estimate the sensible heat flux operating close to the ground over short canopies is pending.

Data Availability Statement:
The datasets analyzed in this study will be available upon request to Pedro Gavilán (pedrod.gavilan@juntadeandalucia.es).

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

Appendix A. The Surface Renewal Method to Estimate Surface Fluxes
Visual examination of the trace of a scalar measured over a plant canopy or bare soil reveals a ramp-like pattern which has been associated to an organized low-frequency motion (coherent structure). For the trace of the temperature of the air measured under unstable conditions, a ramp-like pattern is shown in Figure A1 that can be idealized as follows [30]. Consider a macroparcel of air (i.e., a parcel which volume per unit area is large enough to contain all the sources) uniformly heated traveling at a given height above the surface. At some instant, it suddenly moves down to the surface where it remains in contact with the sources for a period of time until, by continuity, it is renewed by another parcel sweeping in from above. While the parcel of air remains in contact with the surface heat transfers from the sources to the parcel. Thus, after a nearly steady period (namely quiescent period), the temperature of the air parcel shows a gradual increase (namely warming period) followed by a sudden drop (namely microfront period) to a temperature baseline representing the temperature of a fresh (descending) parcel of air [30,31,56]. The ramp-like pattern is associated to the role of a coherent structure (namely the signature of a coherent structure) where the ejection phase represents an injection of heat into the surface layer. Given that the air parcel is presumed uniformly heated, the ramp amplitude ( T A ) represents its net increase of temperature. Therefore, given that the volume (V) of the air parcel per unit surface (Su) covers all the sources, the amount of heat (QSR) injected to the surface layer above the canopy can be expressed as QSR = (ρCp αZ T A ) where ρ and Cp are the density and the isobaric specific heat capacity of the air, respectively, and (αZ) is the volume of the air parcel per unit area (i.e., V = Su αZ) where Z is the measurement height above the ground and α a correction parameter. A steady coherent structure will inject a regular amount of heat QSR with a frequency , where k is the von Figure A1. Time course of the temperature of the air measured at high frequency (dot) for an unstable case. T b is a baseline temperature representing the mean temperature of a parcel of air that descended to the surface. Quiescent (q), gradual warming (l) and ejection (f) phases are distinguished while the air parcel remains in contact with the surface for a period τ T . Ideally, the mean temperature time course of the air parcel during each phase is assumed linear showing a ramp-like shape in the trace (signature of a coherent structure characterized by a ramp amplitude, A T , and period τ T ).(A T + T b ) represents the mean temperature of the air parcel when it was renewed.
Given that the air parcel is presumed uniformly heated, the ramp amplitude (A T ) represents its net increase of temperature. Therefore, given that the volume (V) of the air parcel per unit surface (S u ) covers all the sources, the amount of heat (Q SR ) injected to the surface layer above the canopy can be expressed as Q SR = (ρC p αZA T ) where ρ and C p are the density and the isobaric specific heat capacity of the air, respectively, and (αZ) is the volume of the air parcel per unit area (i.e., V = S u αZ) where Z is the measurement height above the ground and α a correction parameter. A steady coherent structure will inject a regular amount of heat Q SR with a frequency 1/τ T ; and therefore, the sensible heat flux (H) can be estimated as H SR = C p ( α Z A T /τ T [30]. When the measurements are taken in, the inertial sublayer (αZ) can be estimated as, αZ = k π z u * τ T ∅ h −1 (ε) 1/2 , where k is the von Kármán constant, z is the measurement heat above the zero-plane displacement, u * is the friction velocity and ∅ h(ε) is the flux-gradient stability function for heat transfer [22]. Invoking similarity, the latter expression to estimate (αZ) is valid for any scalar. Thus, for latent heat flux, LE SR = L (αZ) A w /τ w where L is heat of vaporization and A w and τ w are the amplitude and period in the trace of water vapor density. Clearly, the ramp amplitude in the temperature trace is positive, zero and negative when it is measured under unstable, neutral and stable conditions, respectively. The SR method to estimate surface fluxes requires the analysis of the scalar trace to determine the ramp dimensions; however, a description of these methods is beyond the scope of this paper. A method widely used is based on a ramp model that assumes an instantaneous ejection-sweep phase, that is, a ramp model that neglects the microfront period in Figure A1 [31,56,64]. Accordingly, the mean ramp dimensions can be determined using structure functions of multiple orders (S n ). For the temperature trace where S denotes structure function, n is the order of the structure function (2nd, 3rd, and 5th), m is the number of data points and j is the sampling lag corresponding to the time lag (r), r = j/ f were f is the sampling frequency. The ramp amplitude is obtained by solving for the real roots of the following cubic equation where p = 10S 2 (r) − S 5 (r) S 3 (r) and q = 10S 3 (r). Once the ramp amplitude is solved from Equation (A2), the ramp period is determined where in Equations (A2) and (A3) r is the time lag that maximizes S 3 (r)/r [30], r x . Given that Equations (A2) and (A3) are only valid for time lags r ≥ r x and that the determination of r x depends on the sampling frequency (the slope of S 3 (r)/r is steep until the global maximum is reached), it is recommended to solve Equations (A2) and (A3) for a time lag slightly higher than r x (i.e., such as the next time lag).