2.1. Scintillometer Processing Chain
The functioning of a scintillometer consists of a transmitter that emits an electromagnetic beam through the atmosphere to a receiver located at some distance, which records the intensity of the signal. This electromagnetic wave has a certain wavelength that typically varies between the optical and microwave ranges. One of the scintillometers that work at optical wavelengths is the LAS, while scintillometers that work at microwave wavelengths are named Microwave Scintillometers (MWS). The combination of these two scintillometers is known as the Optical Microwave Scintillometer (OMS), which was the one used in this research. Note that acronyms, variables, and constants are defined in the Abbreviations section.
On their path length, the electromagnetic beam passes through turbulent eddies with differences in temperature, humidity, and, therefore, refractive index, which cause the signal to fluctuate. Thus, signal intensity fluctuations are proportional to the turbulence intensity. These fluctuations are measured on the receiver as signal intensity variances (
σ2ln(I)) [
13,
24]. When deriving these variances from raw data, some spectral cleaning is always performed to get rid of absorption and electronic noise with a band pass filter (BPF). Then, when combining theoretical principles of atmospheric turbulence with the physics of electromagnetic wave propagation (estimated indirectly from intensity fluctuations measured on the receiver as
σ2ln(I)), it is possible to obtain
H and
LvE [
25,
26]. The atmospheric turbulence is quantified by structure parameters, which are a measure of the turbulent energy present in the inertial range of the refractive index (
[m
−2/3]), temperature (
[K
2m
−2/3]), specific humidity (
[kg
2/kg
2 m
−2/3]), and temperature–humidity correlation (
[kg/kg Km
−2/3]). Structure parameters depend on the specific wavelength,
, of the light beam [
27]. Turbulence and wave propagation theory allow one to compute structure parameter
from the variances measured on the receiver, and once the refractive index structure parameter is obtained,
and
can be determined. Finally, using the Monin–Obukhov similarity theory, MOST, which relates structure parameters with heat fluxes,
H and
LvE can be obtained.
According to Wang et al. [
27], the analytical solutions that relate
σln2(I) and
for different wavelengths (
lopt and
lmw for the optical and microwave range, respectively) can be written as:
in which
,
, and
are constants [-], where the subindices
opt,
mw, and
opt,mw refer to optical, microwave, and optical-microwave;
L is the path length [m],
D is the optical scintillometer diameter [m], and
is the Fresnel length [m].
Then,
is estimated by the
σln2(I) measured in the receiver. However, heat fluxes are not directly related to
but to the temperature structure parameter
, for estimating
H, and to the humidity structure parameter
, for estimating
LvE. Hill [
16] derived empirical formulas that relate the
to the structure parameters
,
, and
through the following equation:
where
T is temperature [K],
q is specific humidity [kg/kg], and
Ai are dimensionless coefficients [-] dependent on atmospheric pressure, temperature, and specific humidity. Overbars denote mean values over an interval. In Equation (4),
,
, and
are independent of the electromagnetic beam wavelength, while the rest of the parameters depend on it.
The spectral region, where the LAS is sensitive to temperature (and, consequently, to
), is independent of humidity (
q), whereas for the MWS, there is a dependence on humidity (
) and the cross correlated function (
), and only a small correction is required using the Bowen ratio (e.g., Green et al. [
28,
29]). Since LAS is more sensitive to
and little information is usually available on
and
, additional information on the relationship between
T and
q fluctuations is used when using a standalone LAS [
30]. For this, Moene [
29] added a correction to the effect of
and
over
, known as the Bowen’s correction, which assumes that the correlation coefficient between the temperature and humidity fluctuations (
[-]) is ±1:
where
β is the Bowen ratio, and
cp is the specific heat at constant pressure.
When using an OMS, Equation (4) can be written for each operating wavelength. As a result, there are two equations (for and ) and three unknowns, which correspond to the meteorological structure parameters (, , and ). Thus, three approaches have been proposed to determine the structure parameters.
The first approach, also known as the two-wavelength method, was proposed by Hill [
16]. It consists of assuming that the three structure parameters are not independent and that the
, defined in Equation (6), is 0.8 for the unstable conditions that occur during the day and −0.6 for the stable conditions that occur during the night:
Many investigations have used the two-wavelength method [
11,
28,
31,
32]. The second approach, also known as the bichromatic method, was proposed by Lüdi et al. [
17]. This method is based on the covariance between the LAS and MWS signals to derive the
. The path-averaged
is found by cross-correlating the two electromagnetics signals at different wavelengths that pass through the same air volume [
17]. This approach has the advantage that no assumptions are made to find
. However, the bichromatic approach is less robust than the two-wavelength method, as sometimes it delivers chaotic behaviour or unrealistic
values [
17,
18]. In this method, the equations required to find the three structural meteorological parameters are the two versions of Equation (4), written for each operating wavelength, and Equation (6). Since both Hill [
16] and Lüdi et al. [
17] approaches are bichromatic methods, this research will refer to them as the Hill [
16] and the Lüdi et al. [
17] methods.
The third and most recent approach to determine the structure parameters, also known as the hybrid method, was developed by Stoffer [
18]. This approach detects and replaces unrealistic
values estimated from the Lüdi et al. [
17] method and applies them in the Hill [
16] method to compute heat fluxes. First,
values are estimated using the Lüdi et al. [
17] method, where unrealistic values of
> |±1| are detected and replaced by those proposed by Hill [
16] (0.8 for
> 1 and −0.6 for
< −1). Then,
erratic-sign fluctuations are detected and replaced by realistic
values. An erratic
sign change is defined as a data interval that has an
value of an opposite sign than those of its surrounding intervals. In this case, the unrealistic
value is replaced by that of the preceding data interval. The resulting
values are then used in the Hill [
16] method. Thus, this procedure reduces the chaotic behaviour of the
values obtained from the Lüdi et al. [
17] method, and it produces more robust results than those obtained in the Hill [
16] method [
18]. Note the importance of computing corrected
values, as they indicate the sign of the heat fluxes (
> 0 refers to unstable conditions and, thus, positive sensible heat fluxes).
Once the meteorological structure parameters are found,
and
can be directly related to
H and
LvE, respectively, by combining the MOST similarity functions with the relation between kinematic and dynamic heat fluxes [
13]:
where
ρ is the air density [kg/m
3],
u* is the friction velocity [m/s],
z is the beam effective height [m],
d is the zero-plane displacement [m], and
are the MOST universal similarity functions.
MOST functions depend on
z/, which represents the atmosphere stability.
LO is the Obukhov length [m], which represents the height above the surface where friction and convective turbulence production are equal.
LO is estimated as:
where
κ is the von Kármán constant (0.4),
g is gravity (9.8 m/s
2), and
is the scaling temperature [K].
The scintillometers’ type (LAS/MWS) does not provide information about friction velocity, so flux-profile-relationships must be used. Therefore,
u* can be estimated by the standard Businger–Dyer flux-profile relation [
33,
34]:
where
u [m/s] is the mean wind speed at height
z [m], and Ψ is a function also dependent on
z/. From Equation (10) it is clear that two wind speed measurements at elevations
z1 and
z2 are required to compute
u*. Nonetheless, when locating one elevation at the roughness length, in which the wind speed drops to zero, only one measurement of wind speed is needed to estimate
u*. In this study, this measurement was performed at the elevation of a meteorological station (details below). Therefore, given
,
u, and
z, Equations (7), (9), and (10) can be solved iteratively for
H. Once
and
have been solved,
LvE can be calculated directly from Equation (8).
Although
and
are assumed to be universal, it has been shown that different MOST functions can cause an up to 20% difference in
H estimations [
35]. Wyngaard et al. [
36] found general similarity functions for unstable (
z/LO < 0) and stable (
z/LO > 0) conditions:
where
refers to
for temperature and
for humidity, and
are empirical coefficients. As shown in
Table 1, there are many MOST functions in the scientific literature. Here, the Wyngaard et al. [
36] formulation was selected, as they were used by de Bruin et al. [
34] to assess scintillometry over a vineyard located in an arid region.
When using a standalone LAS, Bowen’s correction shown in Equation (5) allows calculation of
H [W/m
2] without making measurements or assumptions about humidity fluctuations [
10]. However,
LvE [W/m
2] cannot be estimated directly, since
is disregarded. Thus, the surface EBC is typically used to obtain
LvE after
H has been estimated [
39]:
where
Rn is net radiation [W/m
2], and
G is the ground heat flux [W/m
2]. As a convention, it is considered that
Rn is positive when directed into the surface, while the other fluxes are positive when directed away from it. This method will be referred to as the standalone LAS method.
2.2. Study Area, Instrumentation, and Collected Data
The study area is an irrigated vineyard of 12 ha located in Pirque [
40], an agriculture-based province located in the semi-arid Central valley of Chile surrounded by neighbouring land also used for agricultural activities or very dry wasteland (
Figure 1a). Vineyards are north–south oriented in a vertical trellis system with a vegetation height of ~1.6 m, with little variations throughout the study site. The space between rows is 2.45 m and between plants is 1.20 m. This region has a semi-arid climate with typical annual rainfall amounts of about 460 mm and temperatures in summer ranging between 12–30.4 °C versus 4.4–14.3 °C in winter. Our field campaign was performed in the summer period between 9 January 2019 and 17 January 2019, where no precipitation and no changes in the grapes’ growth stage were observed. An EC system (IRGASON, Campbell Scientific Inc., Logan, UT, USA) and a meteorological station were set up in the middle of the vineyard 4 m above ground. The average wind speed observed over the vineyard was 1.5 m/s, with a maximum of 4.5 m/s and a minimum of 0.1 m/s. The highest wind speeds came mainly from the west (
Figure 1a) over the afternoons. The 20-Hz EC data were processed to obtain 30-min fluxes using the EddyPro 6.2.2 software and activating the default recommended correction procedures [
41,
42]. The meteorological station included a temperature and relative humidity sensor (HMP45C, Campbell Scientific Inc., Logan, UT, USA), a net radiometer (NR Lite 2, Kipp & Zonen, Delft, The Netherlands), four soil heat flux plates (HFP01-L, Campbell Scientific Inc., Logan, UT, USA), two soil thermocouples (TCAV, Campbell Scientific Inc., Logan, UT, USA), and a water content reflectometer (CS616, Campbell Scientific Inc., Logan, UT, USA) All these sensors were connected to a datalogger (CR3000, Campbell Scientific Inc., Logan, UT, USA).
An OMS was installed at 3.02 m height in a NE–SW orientation over a 480-m path. The EC system and the meteorological station were located downwind near the center of the OMS path (
Figure 1b). Under these conditions, the footprint of both systems covered the same type of vegetation and did not extend beyond the vineyard (see
Appendix A). Therefore, the estimates of
H and
LvE from the EC and the OMS were comparable, and ET fluxes at the plot scale could be determined. The OMS was composed of a LAS (LAS MkII, Kipp & Zonen, Delft, The Netherlands), with an aperture of 149 mm and a wavelength of 850 nm, and an MWS (RPG-MWSC-160, RPG Radiometer Physics GmbH, Meckenheim, Germany), with an aperture of 300 mm and a wavelength of 1.86 mm. The zero-plane displacement height was 1.12 m, and the roughness caused by the crops was 0.2 m [
43]. The zero-plane displacement height was estimated from the crop averaged height (
as
[
44]. 1 kHz raw scintillometer data were saved to perform the spectral cleaning method.
As shown in
Figure 1, the vineyard is not completely homogeneous, as it has rows of plants and rows of bare soil. Nonetheless, considering the small horizontal scale of these heterogeneities and that the height of the plants is fairly constant, the experimental configuration of both the EC and OMS was believed to be at or above the blending height. The blending height is the elevation at which the turbulent signatures of each surface are mixed and above which MOST is generally accepted to be applicable [
11]. Moreover, our study site resembled that of Ezzahar et al. [
45], in which sensible and latent heat fluxes were successfully determined over olive trees. In their study, the ratio between the height of the trees and the instruments was 0.43, whereas in this study, this ratio was 0.53 for the OMS and 0.40 for the EC.
2.4. Spectral Cleaning of the Raw Signal Intensities
Unwanted contributions to the OMS intensity variances can result in unrealistic heat fluxes [
22,
46]. The proposed spectral data cleaning process is shown in
Figure 2. The first step in this process is to identify unwanted contributions to the 1-kHz scintillometer raw signal by performing a spectral analysis. The spectral analysis allows one to detect absorption at low frequencies, erratic spikes and tripod vibrations [
47,
48] within the scintillometer frequencies, and electronic noise in high frequencies. The second step is to apply the electronic noise and absorption filter presented in
Figure 2. To eliminate these unwanted contributions, Stoffer [
18] recommends applying BPFs to the original spectrum: a 0.1 Hz high-pass filter (HPF) and a 100 Hz low-pass filter (LPF). Thus, any frequency contribution below 0.1 Hz or above 100 Hz is cut off from both LAS and MWS signals. Here, we applied Stoffer’s [
18] BPFs and assessed the recovered signal variance to ensure that the signal contained only contributions from the scintillation properties related to the heat fluxes. Based on Clifford [
49] scintillation spectra, the application of the 0.1-Hz HPF with a minimum cross wind of ~0.1 m/s recovers more than 90% of the signal variance for the LAS and more than 90% for the MWS. Based on Clifford [
49] scintillation spectra, the application of the 0.1-Hz HPF with a minimum cross wind of ~0.1 m/s recovers more than 90% of the signal variance for the LAS and more than 90% for the MWS; whereas the 100-Hz LPF with a maximum cross wind of ~4 m/s recovers more than 99% of the LAS and MWS signal variances. Then, the third step is to apply a filtering process of erratic spikes and unwanted contributions from tripod vibrations. This step is not trivial, as typically, tripod vibration frequencies are within the scintillometer frequencies, which are relevant for computing heat fluxes. Therefore, it needs to be carried out carefully to avoid filtering the correct signal that should be used by the LAS and MWS. The method shown in
Figure 2 consists in dividing the spectrum in blocks of 5% of the total data and calculating the block’s median (MED) and mean absolute deviation (MAD).
Then, erratic peaks and tripod vibrations within the data block are defined when the following condition is met:
where
Sf is the value of the spectrum in the frequency range,
a and
b are empirical coefficients, and FT [Hz] is a frequency threshold. Hence, the values of a and b should be defined for each dataset, as they depend on the unwanted contributions to the scintillometer signal that appear in the spectra. The FT also needs to be selected for each set of data (e.g., 2 Hz is generally where tripod vibrations occur). When the erratic peaks of one data block are found, the next block to be analyzed is overlapped to ensure continuity, and the process is repeated until the unwanted contributions of all the blocks are found. To remove the unwanted contributions of erratic peaks and tripod vibrations, their Fourier coefficients are randomly varied by assuming that they distribute uniformly within each data block, so the relevant statistical properties are not modified. As this approach could produce new artificial erratic peaks, an iterative process is carried out until all the unwanted contributions are removed.
Figure 3 shows an example of the spectral data cleaning method. In this example, a random spectrum is used (
Figure 3a). A BPF removes absorption and electronic noise of the original spectrum. As a result, any frequency contributions below the HPF and above the LPF are removed (
Figure 3b). Then, the spectrum is evaluated in blocks to identify unwanted contributions from tripod vibrations and erratic spikes. This process is performed from low to high frequencies (
Figure 3c). As shown in Equations (15) and (16), the criteria to categorize an unwanted contribution will depend on the frequency threshold (FT). Thus, a corrected spectrum without electronic noise, absorption, erratic spikes, and unwanted contributions from tripod vibrations is obtained (
Figure 3d).
From the cleaned spectrum, the variance of the signal is calculated using the integral of the spectrum curve, as the spectrum is scaled such that the area under the spectrum is proportional to the signal variance. Furthermore, the corrected signal time series is estimated using the inverse fast Fourier transform (IFFT). The spectral cleaning method was performed using a time step of 30 min to obtain the clean variances and signals for computing the heat fluxes. The entire dataset (9 January 2019 and 17 January 2019) was considered with the aim of analyzing the proposed method effectiveness for removing unwanted contributions with statistical significance, for which the scintillometer heat fluxes are presented before and after the cleaning process using the Hill [
16] method. In the results presented below, we used different dates for the analysis, depending on the quality of the collected data.