Evaluating four gap-filling methods for eddy covariance measurements of evapotranspiration over hilly crop fields

Estimating evapotranspiration in hilly watersheds is paramount for managing water resources, especially in semiarid/subhumid regions. The eddy covariance (EC) technique allows continuous measurements of latent heat flux (LE). However, time series of EC measurements often experience large portions of missing data because of instrumental malfunctions or quality filtering. Existing gap-filling methods are questionable over hilly crop fields because of changes in airflow inclination and subsequent aerodynamic properties. We evaluated the performances of different gap-filling methods before and after tailoring to conditions of hilly crop fields. The tailoring consisted of splitting the LE time series beforehand on the basis of upslope and downslope winds. The experiment was setup within an agricultural hilly watershed in northeastern Tunisia. EC measurements were collected throughout the growth cycle of three wheat crops, two of them located in adjacent fields on opposite hillslopes, and the third one located in a flat field. We considered four gapfilling methods: the REddyProc method, the linear regression between LE and net radiation (Rn), the multi-linear regression of LE against the other energy fluxes, and the use of evaporative fraction (EF). Regardless of the method, the splitting of the LE time series did not impact the gap-filling rate, and it might improve the accuracies on LE retrievals in some cases. Regardless of the method, the obtained accuracies on LE estimates after gap filling were close to instrumental accuracies, and they were comparable to those reported in previous studies over flat and mountainous terrains. Overall, REddyProc was the most appropriate method, for both gap-filling rate and retrieval accuracy. Thus, it seems possible to conduct gap filling for LE time series collected over hilly crop fields, provided the LE time series are split beforehand on the basis of upslope–downslope winds. Future works should address consecutive vegetation growth cycles for a larger panel of conditions in terms of climate, vegetation, and water status.


Introduction
Hilly watersheds are widespread within coastal areas around the Mediterranean basin, as well as in eastern Africa, India, and China.They experience agricultural intensification since hilly topographies allow water-harvesting techniques that compensate for precipitation shortage (Mekki et al., 2006).Their fragility is likely to increase with climate change and human pressure, especially as water scarcity already limits crop production.In this context, understanding evapotranspiration processes within hilly watersheds is paramount for the design of decision support tools devoted to water resource management (McVicar et al., 2007).Indeed, evapotranspiration (or latent heat flux, LE) directly drives biomass production through intertwining with photosynthesis (Olioso et al., 2005), and it is a major term of water balance, up to twothirds of the annual water balance for semiarid/subhumid Mediterranean climates (Montes et al., 2014;Moussa et al., 2007;Yang et al., 2014).
Using evapotranspiration measurements for environmental and water sciences requires complete time series of LE at the hourly timescale, to be next converted into daily, monthly, or Published by Copernicus Publications on behalf of the European Geosciences Union.
annual values (Falge et al., 2001a, b).This is a prerequisite for long-term studies in relation to global change, but also for short-term studies in relation to agricultural issues and modeling challenges.However, common time series of eddy covariance (EC) measurements, which are today considered to be the reference method, include missing data because of experimental troubles such as power failures or instrumental malfunctions.Also, unfavorable micrometeorological conditions lead to the rejection of significant parts of data that do not fulfill theoretical requirements for EC measurements.Statistical studies based on long-term measurements suggest that missing data rates range from 25 to 35 % (Baldocchi et al., 2001;Falge et al., 2001a;Law et al., 2002), while data rejection rates through quality control range from 20 to 60 % (Papale et al., 2006).Therefore, gap-filling methods are necessary to obtain continuous time series of land surface energy fluxes.
Most existing gap-filling methods are devoted to carbon dioxide (CO 2 ) measurements (Aubinet et al., 1999;Falge et al., 2001a;Goulden et al., 1996;Greco and Baldocchi, 1996;Grünwald and Bernhofer, 1999;Moffat et al., 2007;Reichstein et al., 2005;Ruppert et al., 2006).Table 1 summarizes the few studies that addressed measurements of LE, along with underlying methodologies and resulting performances.Prior to gap filling, time series are often split in different ways according to the experimental conditions (e.g., nighttime or daytime, wind directions, vegetation phenology, weekly or monthly time windows), so that missing data are filled with observations collected in similar conditions for micrometeorology, vegetation phenology, and water status.Overall, gap-filling methods for LE time series have been evaluated over flat, hilly, and mountainous areas.However, the existing studies for hilly areas did not address their specific conditions (Hui et al., 2004), or were restricted to one gap-filling method (Zitouna-Chebbi, 2009).
Gap-filling methods for LE have to be designed in accordance with the terrain characteristics that impact evapotranspiration.Conversely to flat terrains that correspond to a slope lower than 2 % (Appels et al., 2016), solar and net radiations within sloping terrains change depending on slope orientation, with larger values for ecliptic-facing slopes (Holst et al., 2005).Over sloping terrains, the conditions of topography, atmospheric stability, and airflow within the atmospheric boundary layer differ between hilly areas and mountainous areas (Dupont et al., 2008;Hammerle et al., 2007;Hiller et al., 2008;Prima et al., 2006;Raupach and Finnigan, 1997).Most existing gap-filling methods rely on covariation of convective fluxes with meteorological variables, or temporal autocorrelation of the convective fluxes.Within hilly areas, these relationships are likely to change with wind direction and vegetation development because of changes in airflow inclination (Zitouna-Chebbi et al., 2012, 2015), and therefore changes in aerodynamic properties (Blyth, 1999;Rana et al., 2007).
In the context of obtaining continuous time series of evapotranspiration from EC measurements of LE, the current study aimed to examine and compare LE gap-filling methods over hilly crop fields.For this, we evaluated the performances of different methods before and after tailoring to the conditions of hilly crop fields.We used the following methodological framework.
The experiment was set within a Tunisian agricultural hilly watershed with rainfed crops.It included the data collection and preprocessing, the analysis of the experimental conditions, and the analysis of the dataset to be filled.
We considered several gap-filling methods that differ in the use of ancillary information, either micrometeorological data or energy flux data other than LE.Given the possible influence of airflow inclination, the gap-filling methods were tailored by splitting the dataset on the basis of airflow inclination as driven by wind direction.
We assessed the performances of the gap-filling methods by addressing (1) filling rate as compared to missing data after preprocessing, (2) retrieval accuracy on filled data, and (3) quality of gap-filled time series through energy balance closure.
2 The experiment: study site and materials

Experimental site
The Lebna watershed is located on the Cap Bon Peninsula, northeastern Tunisia.It extends from the Jebel Abderrahmane to the Korba Laguna, and it includes the Kamech watershed (outlet at 36 • 52 30 N, 10 • 52 30 E, 108 m above sea level -m a.s.l.) that has an area of 2.7 × 0.9 km 2 (Fig. 1).The El Gameh wadi crosses Kamech from the northeast to the southwest.A hilly dam (140 000 m 3 nominal capacity) is located at the watershed outlet.The Kamech watershed belongs to the environmental research observatory OMERE (French acronym for Mediterranean Observatory of Water and Rural Environment, http://www.umr-lisah.fr/omere,last access: 1 April 2018).
The climate of the Kamech watershed is subhumid Mediterranean.Over the 1995-2014 period, cumulated values at the yearly timescale for precipitation and Penman-Monteith reference crop evapotranspiration (Allen et al., 1998) are 624 and 1526 mm, respectively.Terrain elevation ranges from 94 to 194 m a.s.l., and terrain slopes range between 0 and 30 %, the quartiles being 6, 11, and 18 % (Zitouna-Chebbi et al., 2012).The soils have sandy-loam textures, and soil depth ranges from a few millimeters to 2 m according to both the location within the watershed and the local topography.These swelling soils exhibit shrinkage cracks under dry conditions during the summer (Raclot and Albergel, 2006).

Feed-forward (FF) artificial neural networks (ANN) with different inputs
Best performance with the following inputs: daily maximum and minimum temperature, daily solar radiation, day of the year.are winter cereals (barley, oat, triticale, wheat) and legumes (chickpeas, fava beans).Land use and parcels are strongly related to topography and soil quality.The watershed includes 273 plots, the sizes of which range from 0.08 to 13.65 ha (0.62 ha on average, with a standard deviation of 1.05 ha).

Measurement locations and experimental period
Three flux stations simultaneously collected measurements of energy fluxes and meteorological variables within three wheat crop fields (Fig. 1): two sloping fields (A, B) and a flat field (C).
Field A was located on the northern rim of the Kamech watershed.It had a fairly homogeneous terrain slope (6 • ) that faced south-southeast and a 1.2 ha area.Field B was adjacent to field A, on the opposite hillside.It also had a homogeneous slope (5.2 • ) that faced north, and it had a 1 ha area.Fields A and B were separated by the northwestern limit of the Kamech watershed.Field C was located in the southeastern part of the Kamech watershed.It had a flat terrain and a 5 ha area.A meteorological station (labeled M in Fig. 1) was located near the watershed outlet.
The experimental period started at the beginning of December 2012 and ended mid-June 2013.It thus covered the full growth cycle within the three wheat crop fields, from sowing (1 December) to harvest (15 May for field A, 19 June for fields B and C).

Instrumental equipment and data acquisition
On fields A, B, and C, each flux station collected measurements of the land surface energy fluxes: net radiation (Rn), soil heat flux (G), and sensible (H ) and latent heat (LE) fluxes.Table 2 displays the type of instruments used for each flux station along with acquisition and storage frequencies, and sensor accuracies according to manufacturers.The sonic anemometers, krypton hygrometers, and air temperature and humidity probes were installed at constant heights above ground level: 1.98 m for field A, 2.0 m for field B, and 2.2 m for field C. The verticality of the sonic anemometers was carefully checked during the experiment with a spirit level, as described by Zitouna-Chebbi et al. (2012).The latter reported a 1 • accuracy on sonic anemometer verticality, according to the experimental protocol and to the analysis of airflow inclination.To avoid water ponding on mirrors of the krypton hygrometers, we rotated them in their mounts so that the mirrors were vertical and the measurement path was horizontal.The net radiometers were installed at 1.7 m height above ground level and their horizontality was also checked regularly.For each flux station, three soil heat flux sensors were distributed a few meters around the station, and they were buried at 5 cm below the soil surface.
Measurements at the meteorological station included (1) solar irradiance with a SP1100 pyranometer (Skye, UK); (2) air temperature and humidity with a HMP45C probe (Vaisala, Finland); (3) wind speed with an A100R anemometer (Vector Instruments, UK); and (4) wind direction with a W200P wind vane (Vector Instruments, UK).The instruments were installed at 2 m above ground level (1 m for the pyranometer).All instruments were connected to a CR10X data logger (Campbell Scientific, USA) that calculated and stored the 30 min averaged values from the 1 Hz frequency measurements.
All instruments were manufacturer calibrated.Hereafter in the paper, we focus on daytime measurements since LE flux is insignificant (small) during nighttime.

Data processing: calculation of net radiation and soil heat flux
On fields A and B, the measurements of Rn were corrected for the effect of slope following the procedure proposed by Holst et al. (2005).Details are given in Zitouna-Chebbi et al. (2012, 2015).Only direct solar irradiance was corrected by accounting for the angle between solar direction and the normal to local topography.Direct solar irradiance was empirically derived from total solar irradiance measured at the flux station.We characterized local topography with slope (topographical zenith with nadir as origin) and aspect (topographical azimuth with north as origin), both derived from a 4 m spatial resolution digital elevation model (DEM) obtained with a stereo pair of IKONOS images (Raclot and Albergel, 2006).The correction for slope effect on Rn was about 50 W m −2 on average.For each flux station, G was estimated by averaging the measurements collected with the three soil heat flux sensors.We did not apply any correction for heat storage between the surface and the sensors for several reasons.First, the existing solutions are questionable when considering swelling soils that exhibit shrinkage cracks under dry conditions during the summer since they require detailed and stable experimental protocols (Leuning et al., 2012).Second, the experiment lasted throughout wheat growth cycles without any flood events that are critical for heat storage correction.Third, neglecting the heat storage in the soil above the heat flux plates induces errors on soil heat flux estimates that are not systematically large since they range between 20 and 50 W m −2 on average (20-50 % relative to measured value), as reported by Foken (2008).

Data processing: calculation of convective fluxes
H and LE were calculated from the 20 Hz data collected by the sonic anemometers and the krypton hygrometers, using the ECPACK library version 2.5.22 (Van Dijk et al., 2004).H and LE fluxes were calculated over 30 min intervals.

Flux calculation
Most of the instrumental corrections proposed in the aforementioned version of the ECPACK library were applied.These corrections addressed (1) the calibration drift of the krypton hygrometers using air humidity and temperature measured by the HMP45C probes, (2) the linear trends over the 30 min intervals, (3) the effect of humidity on sonic anemometer measurement of temperature, (4) the hygrometer response for oxygen sensitivity, (5) the mean vertical ve-  6) the corrections for path averaging and frequency response (spectral loss), and ( 7) the rotation correction for airflow inclination (see Sect. 2.5.2).

Coordinate rotations
When calculating energy fluxes with the EC method, it is conventional to rotate the coordinate system of the sonic anemometer (Kaimal and Finnigan, 1994).Coordinate rotations were originally designed to correct the vertical alignment of the sonic anemometer over flat terrains, and they are commonly used over non-flat terrains to virtually align the sonic anemometer perpendicularly to the mean airflow, in an idealized homogeneous flow.
The main potential drawback of the double rotation method is that a significant variability in rotation angles can be observed at low wind speeds (Turnipseed et al., 2003).Since our study area was typified by large wind speeds (Zitouna-Chebbi et al., 2012, 2015), we selected the double rotation method that is applied to each time interval over which the convective fluxes are calculated (30 min in our case).After a first rotation (yaw angle) that cancels the lateral component of the horizontal wind speed, a second rotation (pitch angle) is applied around a horizontal axis perpendicular to the main wind direction to cancel the mean vertical wind speed.Thus, it implicitly accounts for changes in wind direction and vegetation height that are likely to be constant over 30 min intervals.

Data quality assessment
Several quality criteria for flux measurements have been proposed in the literature.The most commonly used are the steady-state (ST) test and the integral turbulence characteristics (ITC) test (Foken and Wichura, 1996;Geissbühler et al., 2000;Hammerle et al., 2007;Rebmann et al., 2005).These tests verify that the theoretical requirements for the EC measurements are fulfilled.The ST test assesses the stationarity of the turbulence, while the ITC test assesses the good development of the turbulence, both tests being performed over each 30 min interval.Although established over flat terrains, they have been used for a long time over mountainous terrains (Hammerle et al., 2007;Hiller et al., 2008) and more recently over hilly terrains (Zitouna-Chebbi et al., 2012, 2015) because there is no specific test for relief conditions.
Quality classes were assigned to each half-hourly flux datum according to the results of the two tests.For this, we followed the classification proposed by Foken et al. (2005) and Rebmann et al. (2005).H and LE flux data belonging to the quality class I could be used for turbulence studies.H and LE flux data belonging to classes II to IV could be used for long-term flux measurements.Finally, we rejected H and LE flux data belonging to class V that correspond to both ST > 0.75 and ITC > 2.5.
Regarding footprint, the flux contributions were likely to originate from the target fields, regardless of wind direction and vegetation height.On the one hand, experimental conditions (measurement height, field size, vegetation height, and micrometeorology) were similar to those indicated in Zitouna-Chebbi et al. (2012, 2015).On the other hand, the latter reported that calculated flux contribution from the target fields was about 75-80 % throughout three 1-year duration experiments.In the next section, we address the vegetation and micrometeorological conditions, as well as the subsequent relevance of measurement height.

Climate forcing and wind regime
During the experiment that lasted from December 2012 to June 2013, the meteorological station (M) recorded a cumulative precipitation amount of 563 mm.Over the same period, the reference evapotranspiration ET 0 recorded by the meteorological station ranged between 1.1 and 5.8 mm day −1 , with a cumulated value of 510 mm.The wind speed value recorded during the experimental period by the meteorological station was 4 m s −1 on average.The averaged wind speed value recorded by the meteorological station was very close to those recorded by the sonic anemometers installed on the flux stations within fields A, B, and C, with differences lower than 0.4 m s −1 .The spatial homogeneity for wind speed was also observed in previous studies conducted at different locations within the same watershed (Zitouna-Chebbi, 2009;Zitouna-Chebbi et al., 2012, 2015).
The wind rose obtained from the data collected at the meteorological station depicted two prevailing directions (Fig. 2).The first direction corresponded to winds coming from the south (directions between 70 and 220 • , clockwise; north is 0 • ).The second direction corresponded to winds coming from the other directions, hereafter referred to as northwest winds.The topography induced downslope winds on field A and upslope winds on field B under northwest winds.The reverse was observed under south winds.
Micrometeorological conditions were analyzed using the atmospheric stability parameter ξ = (z − D)/L MO , where z is measurement height (in meters), D is displacement height (in meters), and L MO is Monin-Obukhov length (in meters).D was set as two-thirds of vegetation height, the latter being derived from in situ measurements (see Sect. 2.6.2).The atmospheric stability parameter ξ was negative most of the time, with few values larger than 0.1, mainly during sunrise or sunset.The ξ median values were −0.007, −0.011, and −0.010 for field A, B, and C, respectively.These values corresponded to conditions of forced convection (neutrality or low instability) induced by large wind speeds.We did not observe notable differences between northwest and south winds.Zitouna-Chebbi et al. (2012, 2015) obtained similar results with a dataset collected between 2003 and 2006 on different fields within the same study area.
Overall, the analysis of wind direction and micrometeorological conditions indicated that the wind regime did not stem from valley wind or land-sea breeze.Indeed, the wind direction did not depict any diurnal course in relation to anabatic or katabatic flows or to land-sea temperature difference, while the ξ parameter did not correspond to conditions of atmospheric stability with free convection.

Vegetation conditions
Throughout the experiment, the evolution of the wheat phenology was monitored using the so-called "BBCH scale improved" of Feekes and Large (Lancashire et al., 1991).Fields A, B, and C depicted similar phenological evolutions.The beginning of the tillering stage appeared on 15 January, and full tillering occurred on 19 February.The start of bolting was on 5 March, and full flowering was on 22 April.The seed maturity stage lasted from the beginning to the end of May, and the beginning of senescence was late May.
Vegetation height was measured on a weekly basis using a tape measure.For each date, 30 height measurements were performed within each field, and were then averaged at the field scale.Vegetation height reached its maximum on 22 April, and maximum averaged values were 1.00, 0.87, and 0.98 m, for fields A, B, and C, respectively.Vegetation height measurements were next interpolated on a daily basis by using a logistic function.
The vegetation height data indicated that the sonic anemometers and KH20 krypton hygrometers, set up around 2 m above soil surface, were located above the roughness sublayer.Indeed, the experimental conditions were typified by neutral or slightly unstable conditions that corresponded to a roughness sublayer extension from the ground up to 1.43 × vegetation height (Pattey et al., 2006).
Leaf area index of green leaves (GAI) was measured using a planimeter.Every 2 weeks, all leaves were collected within three 1 m long transects to derive a spatially averaged value.GAI reached its maximum on 11 April, and maximum values were 2.5, 2.3, and 2.3 m 2 m −2 for fields A, B, and C, respectively.

The dataset to be filled
Missing LE data stemmed from (1) total shutdowns of flux stations following battery discharges or vandalism acts, (2) malfunctions of KH20 krypton hygrometers after precipitation events when air humidity permeated the sensor because of seal degradation, and (3) rejection of LE data identified as class V data by ST and ITC tests (Sect. 2.5.3).
Table 3 displays the number of available data derived from EC measurements over the three fields, when considering Table 3. Summary of the available latent heat flux (LE) data derived from the eddy covariance measurements conducted on each of the three fields: dates of the beginning and ending of measurement periods, number of total daytime data over 30 min intervals for calculating the fluxes, numbers and percentages of data belonging to quality control classes (I-IV: good quality data; V: rejected data), numbers and percentages of the missing data due to dysfunctions of the KH20 sensors, and numbers and percentages of missing flux data due to total shutdowns of the flux stations.LE.It gives the beginning and end dates of the EC measurements, the number of daytime data over 30 min intervals, the numbers and proportions of data with good (classes I to IV) and bad quality (class V) according to ST and ITC tests, the number of missing data due to malfunctions of the Krypton hygrometer (KH20), and the number of missing data because of total shutdown of flux stations.

Number of
The ratio of filtered to original data ranged between 20 and 61 %.It was rather low compared to the ratios reported by former studies at the yearly timescale for worldwide flux networks such as FLUXNET (65 %), in which these ratios stemmed from system failures or data rejection (Baldocchi et al., 2000;Falge et al., 2001a, b).The low ratio we obtained in the current study was ascribed to KH20 malfunctions and total shutdown of flux stations.Furthermore, the KH20 sensor installed on field B was out of order from the end of March until the end of the experiment because of severe instrumental malfunctions.
The proportion of bad-quality data was low, with around 3 % of data belonging to class V.The results of the quality control tests did not exhibit any difference among the fields.For H , the percentages of data belonging to the high-quality classes (I to IV) were 85, 84, and 88 % for fields A, B, and C, respectively.
On the one hand, the rate of missing data for the current study, between 40 and 80 %, was much larger than that reported in former studies, i.e., between 25 and 35 % (Baldocchi et al., 2001;Falge et al., 2001a;Law et al., 2002).On the other hand, the rate of data rejected by quality control, between 2 and 4 %, was much lower for the current study compared to that reported in former studies, i.e., between 20 and 60 % (Papale et al., 2006).Therefore, the overall rate of data to be filled was comparable to that reported in former studies.

Rationale in choosing and implementing gap-filling methods
Amongst the existing LE gap-filling methods listed in the Introduction (Table 1), we selected some methods that differ in the use of ancillary information, either meteorological variables or energy fluxes.The meteorological data to be used were those provided by the meteorological station, while the flux data to be used were those collected at each of the three flux stations of interest (Sect.2.3).We did not select methods that involve measurements of soil water content or vegetation canopy since energy fluxes indirectly account for the latter at a spatial scale closer to that of the LE missing data (see results about footprint analysis in the last paragraph of Sect.2.5.3).Amongst the existing LE gap-filling methods listed in the Introduction (Table 1), we selected the commonly used REd-dyProc method that relies on lookup tables (LUTs) and mean diurnal variation (MDV) to fill missing flux data with those collected under similar meteorological conditions or with averaged values over adjacent days.We also selected methods that fill LE gaps by using multilinear regressions on other energy flux data (Rn, H , and G).We did not select methods based on artificial neural networks because ensuring the relevance of calibration, testing, and validation steps requires large datasets of at least 1 year (Abudu et al., 2010;Beringer et al., 2007;Eamus et al., 2013;Papale and Valentini, 2003).

REddyProc
For the REddyProc method, we selected the online tool available at https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWebRPackage (last access: 1 April 2018), and that is based on Reichstein et al. (2005).The REd-dyProc method combines the covariation of the convective fluxes with meteorological variables (Falge et al., 2001b) and the temporal autocorrelation of the convective fluxes (Reichstein et al., 2005).Gaps are filled in accordance with available information by considering three cases: (1) solar radi-ation (Rg), air temperature (T air ), and vapor pressure deficit (VPD) data are available; (2) Rg data only are available; and (3) none of the Rg, T air , and VPD data are available.
For Case (1), the missing LE value is replaced by the average value under similar meteorological conditions within a time window of ±7 days.Similar meteorological conditions correspond to Rg, T air , and VPD values that do not deviate by more than 50 W m −2 , 2.5 • C, and 5 hPa, respectively.If no similar meteorological conditions occur within the ±7day time window, the latter is extended to ±14 days.
For Case (2), a similar approach is taken.Similar meteorological conditions correspond to Rg deviation by less than 50 W m −2 , and the window size is not extended.
For Case (3), the missing value is replaced by an adjacent value within ±1 h, or by an averaged value at the same time of the day that is derived from the mean diurnal course over ±1 day.
In case the three steps do not permit filling the gaps, the whole procedure is repeated while increasing the window sizes until the value can be filled.Thus, the window size increases using 7-day steps until ±70 days for Cases 1 and 2, and until ±140 days for Case 3, which obviously results in a degradation of the quality indicator.

LE reconstructed from Rn
Initially proposed by Cleverly et al. (2002), this method was successfully tested on our study site by Zitouna-Chebbi (2009).It assumes the stability of the LE / Rn ratio over a given period that can be 1 day, 1 month, or 1 year (Table 1).We implemented the method by first calibrating the linear regression LE = a Rn + b on existing LE and Rn data, and next applying the regression to missing LE data for which Rn was actually measured.This method will be referred to as "LE-Rn method" hereafter.

LE reconstructed from multi-linear regression against other energy fluxes
This method is an extension of the LE-Rn method since LE is estimated as a linear combination of the other energy fluxes Rn, H , and G.As for the LE-Rn method, the multi-linear regression (MLR) method was implemented by first calibrating the multi-linear regression on existing LE, Rn, H , and G data (LE = a Rn + b G + c H + d ), and next applying the regression to missing LE data for which the three other fluxes were actually measured.This method will be referred to as the "MLR method" hereafter.Energy balance theoretically implies a = 1, b = −1, c = −1, and d = 0.However, this is not the case in practice because of the energy imbalance problem for EC measurements.This problem has been mentioned in the literature for vegetated canopies and bare soils, as well as for flat, mountainous, and hilly terrains (Foken, 2008;Hammerle et al., 2007;Leuning et al., 2012;Wilson et al., 2002;Zitouna-Chebbi et al., 2012, 2015).As reported by Leuning et al. (2012), the energy imbalance problem is that the sum of the convective flux (H + LE) underestimates available energy (Rn − G) because of theoretical assumptions (e.g., neglecting storage terms or lateral turbulent transfers) and because of experimental assumptions (e.g., neglecting measurement inaccuracies, neglecting differences in measurement spatial extensions).Thus, applying the energy balance equation LE = Rn− G − H would transfer energy imbalance onto LE estimates, which is not the case with the MLR method that involves a regression calibration (a = 1, b = −1, c = −1 and d = 0).

LE reconstructed from evaporative fraction
Evaporative fraction (EF) is defined as the ratio of LE over available energy (Rn − G) when assuming the latter equals the sum of convective fluxes (H + LE).Li et al. (2008) and Shuttleworth et al. (1989) showed that EF was almost constant during daytime hours.Although rebutted (Hoedjes et al., 2008;Van Niel et al., 2011), various studies stated that EF at midday (EF md ) is statistically representative of daily EF, and thus recommended using EF md for estimating LE (Crago and Brutsaert, 1996;Crago, 1996;Gentine et al., 2011;Li et al., 2008;Peng et al., 2013).
The estimation of missing LE data was twofold.In a first step, EF md was calculated on a daily basis by using the measured data over the 4 h centered on solar noon, provided that at least 75 % of the eight 30 min data were available between noon −2 h and noon +2 h for LE, Rn, and G.
In a second step, the missing LE data were estimated as LE = (Rn− G) EF md , when Rn and G were actually measured.This method will be referred to as the "EF method" hereafter.
Compared to the MLR method that implicitly accounts for the energy imbalance problem via the regression calibration, the EF method induced an overestimation of the convective fluxes by replacing H + LE with available energy Rn− G. Conversely, averaging EF around solar noon rather than over the diurnal course might induce an underestimation of EF at the daily timescale.Therefore, the EF method was likely to (1) induce some errors on LE estimates used for filling gaps and (2) increase energy imbalance for the reconstructed data because of the difference between H + LE and Rn − G.

Tailoring the gap-filling methods to the conditions of hilly crop fields
The gap-filling methods were tailored to the conditions of hilly crop fields by splitting the dataset on the basis of the airflow inclination that is driven by the combined effect of Table 4. Splitting of the dataset into three periods when implementing the LE-Rn and MLR gap-filling methods.The three periods are labeled green vegetation (GV), pre-senescent vegetation (PS), and senescent vegetation (SV).They are indicated along with the vegetation and climatic conditions.GAI stands for leaf area index of green leaves; ET 0 stands for the reference evapotranspiration.Minimum and maximum GAI values are averaged values over the three fields A, B, and C. Cumulative precipitation, mean ET 0 , and mean air temperature are derived from measurements at the meteorological station.wind direction, topography, and vegetation height.The analysis of the experimental conditions showed that the wind regimes were typified by two main wind directions, i.e., northwest and south, which induces upslope and downslope winds on fields A and B (Sect.2.6.1).Therefore, any of the three datasets for fields A, B, and C were split into two sub-datasets that correspond to northwest and south winds.We recall that (1) northwest winds correspond to downslope and upslope winds on fields A and B, respectively; (2) south winds correspond to upslope and downslope winds on fields A and B, respectively; and (3) field C was horizontal.

Main
Most existing gap-filling methods for LE measurements include a prior splitting of the time series to be filled (Table 1), so that missing data are filled with existing observations collected under similar conditions.REddyProc relies on time windows up to 280 days, and the EF method relies on a daily time window.For both the LE-Rn and MLR methods that assume stable regressions among energy fluxes, we split the dataset into three periods that differed in vegetation phenology.This accounted for changes in soil water content and vegetation height at monthly to seasonal timescales.The beginning and end of each period are given in Table 4, along with the vegetation and climatic conditions.The first period corresponded to active green vegetation, with moderate reference evapotranspiration, and with abundant and frequent precipitation events that supply plant transpiration and soil evaporation.It was typified by the absence of water stress, and therefore large values for both EF and the LE / Rn ratio.We labeled this first period "GV" for green vegetation.The second period preceded grain maturation and leaf senescence.It corresponded to the beginning of water stress that resulted from the combined effect of limited precipitation and large reference evapotranspiration.We labeled this second period "PS" for pre-senescent vegetation.The third pe-riod corresponded to leaf senescence and grain maturation.It corresponded to a pronounced water stress that resulted from the combined effect of no precipitation and large reference evapotranspiration.We labeled this third period "SV" for senescent vegetation.

Assessing the performances of the gap-filling methods
The performances of the three gap-filling methods were assessed on filling rate, retrieval accuracy, and quality of gapfilled time series through energy balance closure.In order to make the performances of the four methods comparable, we used the following procedure.Conversely to REddyProc, the LE-Rn, MLR, and EF methods were not able to fill gaps induced by total shutdowns of the flux stations.Therefore, we addressed the filling of the gaps that resulted from malfunctions of the KH20 sensors and quality filtering only.
For field B, the LE-Rn, MLR, and EF methods were not able to fill gaps induced by the shutdown of the KH20 sensor from the end of March (middle of the GV period) to the end of the experiment.Indeed, the EF method required Rn, G, and LE data on a daily basis, while the LE-Rn and MLR methods required data for each of the periods GV, PS, and SV, which excluded periods PS and SV.Therefore, we disregarded the time period in question (from the end of March to the end of the experiment) for field B.
The filling performances were given in accordance with the number of reconstructible data (LE missing data because of both KH20 malfunctions and quality filtering).They were expressed as the ratio of reconstructed to reconstructible data.
The prior splitting of the time series to be filled is a common procedure for most gap-filling methods (Table 1), but it is different from one method to another (Sect.3.2).Therefore, we did not assess the performances of the gap-filling methods on the basis of the time periods GV, PS, and SV.We discriminated the periods GV, PS, and SV for the regression calibrations only (LE-Rn and MLR methods).
To quantify retrieval accuracy, REddyProc provides estimates for each existing datum, where the estimate is derived independently of the corresponding data.Therefore, we implemented a leave-one-out cross-validation (LOOCV) procedure to evaluate the retrieval accuracy for the LE-Rn, MLR, and EF methods.For this, any estimate for retrieval accuracy was calculated by removing the corresponding reference value.
We evaluated the performances of the gap-filling methods before and after splitting the time series on the basis of wind direction (northwest or south).We separately considered the fields A, B, and C, where fields A and B are located on two opposite hillsides with upslope and downslope winds, and field C is located on a horizontal terrain.
The retrieval accuracy was quantified using absolute and relative root-mean-square error (RMSE and RRMSE) as well as mean absolute difference (MAE) and bias and coefficient of determination R 2 (Jacob et al., 2002;Moffat et al., 2007).
To evaluate the quality of the gap-filled time series, we compared the sum of the convective energy fluxes (H + LE) against available energy (Rn − G) before and after gap filling, where gap filling was conducted after splitting the time series on the basis of wind direction.Although energy balance closure analysis is questionable for assessing the consistency of flux measurements, it permits the comparison of independent measurements.

Filling performances of the gap-filling methods
For the three fields (A, B, C) and the two wind directions (northwest, south), Table 5 displays the number of reconstructible data (LE missing data because of KH20 malfunctions or LE data belonging to quality class V), as well as the number and percentage of reconstructed data by the four methods (REddyProc, LE-Rn, MLR, and EF).For each field, the total number of reconstructible data is also indicated, as well as the total number and corresponding percentage of reconstructed data.The total number of reconstructible data in Table 5 corresponds to that given in Table 3 (i.e., sum of LE missing data because of KH20 malfunctions and of LE data belonging to quality class V), apart from field B (2083 versus 3060) for which we restricted the time period to the GV period since no LE data were available on periods PS and SV because of the KH20 shutdown (second item in Sect.3.3).Columns 1 and 2 correspond to upslope and downslope winds, respectively.Lines 1, 2, and 3 correspond to the three periods (GV, PS, SV) that differed in vegetation phenology, soil water content, and climatic conditions.The dashed line is the 1 : 1 line, and the continuous line is the regression line.R 2 is the coefficient of determination.RMSE and RRMSE are absolute and relative root-meansquare errors, respectively.N is the number of flux data calculated over 30 min intervals.
With both the REddyProc and LE-Rn methods, all the missing LE data could be reconstructed.The MLR method permitted the reconstruction of 84, 86, and 90 % of the missing LE data, on fields A, B, and C respectively.The EF method permitted the reconstruction of 32, 19, and 70 % of the missing LE data, on fields A, B, and C, respectively.The reconstruction rates obtained with the MLR method were similar on fields A, B and C. Conversely, the reconstruction rate with the EF method was much larger on field C (flat terrain) than that on fields A and B (sloping terrains).Overall, the filling rate was the same for a given field, whether we split the time series on the basis of wind direction or not.

Accuracy of the gap-filling methods
The calibration of the LE-Rn method for the three periods (GV, PS, and SV) was similar for fields A (Fig. 3), B, and C Table 5. Filling rate performance for each of the four gap-filling methods (REddyProc, LE-Rn, MLR, EF), expressed as the number and percentage of reconstructible LE data that were actually reconstructed.The filling rates are given for each field (A, B, C) and each wind direction.S and NW stand for south and northwest winds, respectively.Up and down stand for upslope and downslope winds, respectively, for the sloping fields A and B (field C was flat).
(Fig. S1a and b in the Supplement).The LE / Rn ratio exhibited a notable temporal stability for each of the three periods, and we did not observe any distinct scatter for the period GV, even if the scattering was larger compared to the periods PS and SV.However, we observed significant differences in slope and offset from one period to another, with changes in slope between 90 and 170 % (relative to mean value), and changes in offset between 60 and 120 % (relative to mean value).We obtained similar LE-Rn regressions for fields A (Fig. 3) and B (Fig. S1a) when splitting the time series on the basis of south and northwest winds that correspond to upslope (respectively downslope) and downslope (respectively upslope) winds on field A (respectively B).Apart from the SV period with too few data on field A, we noted some differences in regressions between the two wind directions for any period, with changes in slope between 5 and 50 % (relative to mean value) and changes in offset between 40 and 80 % (relative to mean value).Conversely, the differences were lower on field C with a flat terrain (see Fig. S1b), with changes in slope between 0.5 and 10 % (relative to mean value) and changes in offset between 15 and 30 % (relative to mean value).A covariance analysis conducted on the regression coefficients showed that the changes in slope and offset were statistically significant in most cases (Table S1).
We quantified the retrieval accuracies of the four gapfilling methods by comparing reference data and gap-filling retrievals of LE over 30 min intervals for each field and each wind direction (Fig. 5 and Table S2).The retrieval accuracies were obtained using a LOOCV procedure (Sect.3.3).We observed the following trends.
The four methods provided similar retrieval accuracies, with differences among RMSE values lower than 20 W m −2 .Bias values were almost null, apart from the EF method.To a lesser extent, the RMSE values were lower with REd-dyProc, which also provided better R 2 values, and the EF method provided the larger RMSE and bias values, down to −20 W m −2 for bias.
Regardless of the gap-filling method, the retrieval accuracies were similar for fields A and C, whereas they were lower for field B.
The method performances could be either different or similar before and after the splitting of the time series on the basis of wind direction.For field A, the RMSE values were similar for upslope and downslope winds, and they were comparable to those obtained before the splitting.For field B, the RMSE values were much lower (respectively slightly larger) for downslope winds (respectively upslope winds) compared to those obtained before splitting the time series.For field C with a flat terrain, the statistical indicators were comparable before and after the splitting.

Energy balance closure analysis
We recall that the gap-filling retrievals we considered for energy balance closure analysis were those obtained by splitting the time series on the basis of wind direction (Sect.3.3).We obtained similar results for energy balance closure for field A (Fig. 4) and fields B and C (Fig. S2a and b).Before and after gap filling, the sum of the convective fluxes systematically underestimated available energy, apart from field B after gap filling with the EF method.On a field basis, change in energy balance closure from one gap-filling method to another was 15 % for field A, 65 % for field B, and 44 % for field C, according to changes in the H + LE versus Rn − G regression slope.On a method basis, energy balance closure varied from 5 % (MLR) to 32 % (EF) from one field to another, according to changes in the H + LE versus Rn − G regression slopes.Finally, energy balance closure could be better after gap filling, and energy balance closure on slop-  ing fields A and B was comparable to that on the flat field C.
When comparing energy balance closure after gap filling with the four methods, we could not identify any clear trend on the basis of the (H + LE) versus (Rn − G) linear regression.Gap filling with the LE-Rn method provided among the best energy balance closure, and gap filling with the REd-dyProc method provided among the worst energy balance closures.Energy balance closure was very similar for the LE-Rn and MLR methods, with changes in the regression .Accuracy of LE retrievals for the four gap-filling methods (REddyProc labeled as REP, LE-Rn, MLR, EF).Fluxes were calculated over 30 min intervals.Retrieval accuracy is given for each field (A, B, C) and each wind direction (All stands for all data.NW and S stand for northwest and south winds, respectively).Accuracy is quantified using statistical indicators (absolute RMSE, bias, coefficient of determination R 2 ).slope between 2.5 % (field C) and 4.5 % (field A).Further, the EF method could provide the worst (field A) or the best (field B) energy balance closure.The scattering around the (H + LE) versus (Rn − G) regressions was reduced after gap filling, either slightly with the REddyProc and EF methods, or a lot with the LE-Rn and MLR methods.

Filling performances of the gap-filling methods
The filling rate was maximal with the REddyProc and LE-Rn methods.Indeed, REddyProc relies on existing LE values within a given time window, either corresponding to similar meteorological variables or derived from averaged diurnal courses.Similarly, the LE-Rn method relied on contin-uous measurements of Rn.The MLR method was less efficient than the REddyProc and LE-Rn methods because of both missing H measurements and H data rejection by quality control.In this case, the filling rate was comparable to the percentage of available H data given in Sect.2.7 (84, 86 and 90 % versus 85, 84, and 88 % for fields A, B, and C, respectively).The worst efficiency of the EF-based gap-filling method was explained by the fact that Rn, G, and LE data around solar noon are required on a daily basis.
The filling rate was similar whether we split the time series on the basis of wind direction or not.For REddyProc, this was explained by the capability of the method to find LE data under similar meteorological conditions or to obtain averaged values from diurnal courses within a scalable time window.For the LE-Rn and the MLR methods, this was explained by existing data for regressions within the three periods GV, PS, and SV, when applicable.For the EF method, this was explained by the daily basis computation of EF and the subsequent filling at the daily timescale.Overall, the four methods were able to fill all gaps in the time series, in spite of larger gap occurrences induced by the splitting of the time series on the basis of wind direction.Also, it is important to note that conversely to the LE-Rn, MLR, and EF methods that relied on energy fluxes (Rn, G, and H ), REddyProc had the capability of filling gaps induced by total shutdowns of the flux stations, although we did not address these total shutdowns to make the performances of the four methods comparable.
We could not compare the filling rates we obtained in the current study against outcomes from the former studies listed in Table 1 for LE data owing to the absence of information on this issue.

Accuracy of the gap-filling methods
When calibrating the LE-Rn method, it was relevant to split the time series into the three periods GV, PS, and SV because of large changes in the LE-Rn regression from one period to another.The strong decrease in LE / Rn ratio from GV to SV was ascribed to the decrease in LE magnitude because of vegetation senescence that was combined with no precipitation and increasing reference evapotranspiration.This emphasizes the impact of changes in soil water content and vegetation canopy at monthly to seasonal timescales.When calibrating the LE-Rn method, it was also relevant to split the time series on the basis of northwest and south winds.Indeed, some differences were observed between the two wind directions for the periods GV and PS, and these differences were larger for sloping terrains (fields A and B) than for the flat terrain (field C).Compared to former studies listed in Table 1, these outcomes were consistent with those from Zitouna-Chebbi (2009).Indeed, the latter reported the need to split time series into distinct periods and wind directions, so that it was possible to take into account changes in aerodynamic conditions for measurements collected within the same study area, over other crop fields, and during other years.
The slightly better accuracies obtained with REddyProc indicates that this method was able to find appropriate LE values under similar meteorological conditions within a given time window, in spite of possible changes in soil water content.LE-Rn and MLR provided very similar accuracies.We expected that MLR would outperform LE-Rn because of the additional inclusion into the regression of G and H fluxes that are driven by vegetation canopy and soil water content.Then, the similar accuracies might result from too large time windows for periods GV, PS, and SV, and especially for period GV with large scattering around the regression line (see for instance Fig. 3 with the LE-Rn regression).The EF method provided lower accuracies.We expected better accuracies with the EF method that filled gaps on a daily basis, and the underperformance might result from the combination of (1) the EF underestimation at the daily timescale when computed between 10:00 and 14:00 solar time, and (2) the overestimation of H + LE by Rn − G as a result of energy imbalance.Overall, the method performances were driven by the temporal dynamics of the local conditions in terms of micrometeorology, vegetation canopy, and soil water content.For instance, large precipitations were likely to induce sharp changes in soil water content, thus giving the EF method an advantage that is based on a daily basis computation, and giving the REddyProc method a disadvantage that relies on similar meteorological conditions or average diurnal courses.
Overall, the RMSE values between reference data and gapfilling retrievals of LE ranged between 20 and 90 W m −2 , and almost two-thirds of these values were lower than 50 W m −2 .The retrieval accuracy was similar for the four gap-filling methods, and it was comparable to those reported by the previous studies listed in Table 1 (e.g., between 25 and 50 W m −2 for RMSE).
Finally, the performances could be better when splitting the time series on the basis of northwest and south winds, with much lower RMSE values for downslope winds.This was not systematic for the sloping fields, but was systematic for all methods when applicable, although these methods involved different information for the reconstruction of the missing data.Thus, our study confirmed that it may be necessary to discriminate upslope and downslope winds when implementing gap-filling methods.This is consistent with reports from Zitouna-Chebbi et al. (2012, 2015), who showed the need to discriminate upslope and downslope winds when correcting the influence of airflow inclination on EC measurements collected over hilly crop fields.

Energy balance closure analysis
For the LE-Rn and MLR methods, energy balance closures were similar, they varied little from one field to another, and they were better than those obtained with the REddyProc www.geosci-instrum-method-data-syst.net/7/151/2018/Geosci.Instrum.Method.Data Syst., 7, 151-167, 2018 and EF methods.This was ascribed to the constraint on energy balance closure when replacing gaps with LE estimates derived from regression between energy balance fluxes (LE versus Rn on the one hand, and LE versus Rn, G, and H on the other hand).Energy balance closure was worse with REddyProc, and it varied much from one field to another.This was ascribed to the lack of constraint on energy balance closure when replacing gaps with LE data collected at different times.For EF, energy balance closure varied much from one field to another, and especially on field B with (H + LE) overestimating (Rn − G).This might be explained by changes in compensation effects between (1) the EF underestimation at the daily timescale when computed between 10:00 and 14:00 solar time and (2) the overestimation of H + LE by Rn − G as a result of energy imbalance.
For the four gap-filling methods, energy balance closure after reconstruction of the LE data was comparable to that observed before gap filling, which showed the consistency of the gap-filled time series.Further, energy balance closure for the two sloping fields (A and B) was comparable to that obtained on the flat field (C), which showed the consistency of the reconstructed data after the splitting of the time series on the basis of upslope-downslope winds.We could not compare the energy balance closures we obtained in the current study against the outcomes from the former studies listed in Table 1 for LE data owing to the absence of information on this issue.Nevertheless, our values of energy balance disclosure ([15-35 %]) were comparable to those reported in the literature ([10-30 %]) for flat, hilly, and mountainous terrains (Foken, 2008;Hammerle et al., 2007;Li et al., 2008;Wilson et al., 2002;Zitouna-Chebbi et al., 2012, 2015).

Conclusions
For the four gap-filling methods we evaluated (REddyProc, LE-Rn, MLR, and EF), the retrieval accuracies were similar and comparable to instrumental accuracies.The filling rate was maximal for REddyProc and LE-Rn, whereas it was lower for MLR and EF.Therefore, the REddyProc and LE-Rn methods were the most appropriate for our study case, in terms of completing time series as much as possible while providing retrievals with good quality.This outcome applied even more for the REddyProc method that is able to fill gaps induced by total shutdowns, although a deeper analysis is needed beforehand to evaluate the retrieval accuracies in such situations.
Our results led us to recommend the splitting of LE time series on the basis of wind direction, prior to the implementation of the gap-filling methods.Indeed, the prior splitting of time series on the basis of wind direction might improve retrieval accuracies, although the benefit was not systematic.In addition, the obtained accuracies on LE estimates after gap filling were comparable to those reported in the literature for flat and mountainous areas, and the same applied for energy balance closure as a consistency indicator for the filled time series.Finally, the splitting of the time series did not impact the gap-filling rate, in spite of larger gap occurrences.Therefore, we conclude that it is possible to conduct gap filling for time series collected over hilly terrains, provided the prior splitting of the time series is applied in an appropriate manner by discriminating upslope and downslope winds.
Our study case is widespread within the Mediterranean basin because of orography and climate conditions within coastal areas across the Mediterranean shores.To a lesser extent, the outcomes of our studies are also of potential interest for hilly watersheds in eastern Africa, India, and China.Conversely, the experiment on which the current study relied lasted over one crop growth cycle only, and we offset this temporal restriction by simultaneously considering three locations that differed much in topographical conditions and resulting airflow inclination.Nevertheless, future works should strengthen the outcomes of the current study by addressing (1) a larger panel of environmental conditions in relation to climate, vegetation type, and water statuses and (2) consecutive vegetation growth cycles.
Data availability.The data used in this study are available according to the data policy of the Environmental Research Observatory OMERE (Molénat et al., 2018; http://www.umr-lisah.fr/omere,last access: 1 April 2018).
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Location of the Kamech watershed within the Cap Bon Peninsula, northeastern Tunisia (a).Kamech has a 0.9 km width and a 2.7 km length.Three-dimensional view of Kamech (b), including locations of the experimental fields (A, B, C) and of the standard meteorological station (M).

Figure 2 .
Figure 2. Distribution of the wind directions and wind speeds throughout the experimental period (December 2012-June 2013), as recorded by the meteorological station.

Figure 3 .
Figure3.Calibration of the LE-Rn gap-filling method on field A. Columns 1 and 2 correspond to upslope and downslope winds, respectively.Lines 1, 2, and 3 correspond to the three periods (GV, PS, SV) that differed in vegetation phenology, soil water content, and climatic conditions.The dashed line is the 1 : 1 line, and the continuous line is the regression line.R 2 is the coefficient of determination.RMSE and RRMSE are absolute and relative root-meansquare errors, respectively.N is the number of flux data calculated over 30 min intervals.

Figure 4 .
Figure 4. Energy balance (EB) closure for field A. Flux data are calculated over 30 min intervals.Statistical indicators correspond to the comparison of convective energy (H + LE) on the y axis against the available energy (Rn − G) on the x axis, before (top left subplot)and after (other subplots) reconstruction of LE data by the four gapfilling methods.The dashed line is the 1 : 1 line, and the continuous line is the regression line.Letters a and b are the slope and the intercept of the linear regression, respectively.R 2 is the coefficient of determination.MAE is the mean absolute error.RMSE and RRMSE are absolute and relative root-mean-square errors, respectively.N is the number of 30 min interval data.
Figure5.Accuracy of LE retrievals for the four gap-filling methods (REddyProc labeled as REP, LE-Rn, MLR, EF).Fluxes were calculated over 30 min intervals.Retrieval accuracy is given for each field (A, B, C) and each wind direction (All stands for all data.NW and S stand for northwest and south winds, respectively).Accuracy is quantified using statistical indicators (absolute RMSE, bias, coefficient of determination R 2 ).

Table 1 .
Summary of relevant studies that deal with the performances of different gap-filling methods for LE time series.Landform classification includes flat, mountainous, and hilly.Dataset splits are based on time window or on specific regimes.NR stands for "not reported".
-One site/NR/salt cedar trees -NR/NR -Random split for calibration/testing over X % of existing data with X =

Table 2 .
Details about experimental setup for each of the three flux stations: type of sensors used, acquisition and storage frequencies, and accuracies as given by manufacturers.HR and T stand for air relative humidity and temperature.