Evaporation and CO2 fluxes in a coastal reef: an eddy covariance approach

ABSTRACT Introduction: We conducted season-long observations of evaporation and carbon flux at the Gulf of Aqaba coast, northern Red Sea. We used the eddy-covariance method with a two-tower setup to measure evaporation rates over land and sea and the advection between them. Using a three-dimensional mass balance approach, we calculated total evaporation as the sum of two main components in our site: horizontal advection and turbulent vertical flux, with half-hourly change of water vapor storage and horizontal flux divergence found to be negligible. Outcomes: Average evaporation rates were 11.4 [mm/day] from April through May (early summer) and 10.5 [mm/day] from June through August (summer). The coastal reef was a CO2 sink over the period of measurements, significantly higher in June through August than in April through May. The main environmental drivers of CO2 flux were humidity, water temperature, sensible heat flux, and wind speed. Discussion: The rates of evaporation near the shore were considerably higher than values reported in other studies typically used to represent the mean for the whole Gulf area. We found that evaporation rates computed by common bulk models approximate the mean values of evaporation but have poor representativeness of the intra-daily temporal variation of evaporation. There was a significant correlation between CO2 flux and evaporation attributed to common environmental drivers of gas diffusion, turbulent fluxes, and horizontal transport. Conclusion: We conclude that observations of fluxes in coastal waters need to use at least a two-tower system to account for the effect of horizontal advection on the total flux.


Introduction
Accurate estimates of evaporation and CO 2 flux over coastal waters are important for understanding biophysical processes at the focus of disciplines such as oceanography, lake hydrology, and marine ecology, among others. Coastal systems are among the most productive ecosystems, fixing considerable amounts of carbon, typically named "blue carbon" (Siikamaki et al. 2013). In coastal marine ecosystems, CO 2 flux can be strongly linked to evaporation because of their common environmental drivers of air-sea exchange such as temperature and salinity (Woolf et al. 2016). In particular, the colder temperature in the top millimeters of the ocean, the "thermal skin," which is a result of evaporation, can have a significant effect on the resulting CO 2 exchange between the water and the air (Robertson and Watson 1992). High rates of evaporation also increase vertical mixing (Woolf et al. 2016), which can have a strong effect on the photosynthetic fluxes in coastal waters due to changes in nutrient concentration, light availability, and water temperature (Genin, Lazar, and Brenner 1995). This, in turn, determines the turbidity of the water column, which is critically important to the existence and abundance of coral reefs, fisheries, and other coastal marine life (Genin, Lazar, and Brenner 1995).
Evaporation can be modeled using either the surface-energy balance or the water mass balance approach. The energy balance is often used as a reference for computations of evaporation (Assouline and Mahrer 1993). However, this method presents challenges because not all the components of the energy balance equation are widely measured in coastal areas. This is particularly true in the case of the heat storage in the water, which requires measurements of temperature and flow rates at different depths throughout the water column. Alternatively, the mass balance approach includes a lower number of variables that are more commonly observed and thus provides a more practical application. The mass balance approach, broadly used in different oceanographic models, typically consist of either bulk parameterizations based on characteristics of the air-sea interface (e.g. Kara, Rochford, and Hurlburt 2000) or formulations based on empirical parameterization of turbulence (e.g. Fairall et al. 1996;Kondo 1975). Neither of these methods, however, have been designed to account for the contributions of horizontal advection of water vapor in the atmosphere toward the nearby dry land.
The most direct method to measure evaporation from open water surfaces is the eddy-covariance (EC) technique Tanny et al. 2008).
EC towers in open waters have been traditionally installed such that the fetch of the observation will allow measuring fluxes coming exclusively from the water and not from the adjacent land (Assouline and Mahrer 1993;Bouin et al. 2012;Tanny et al. 2011). In these cases, single-tower observations provide the turbulent flux and an additional change-of-storage term for the air volume under the tower top, which can be observed through a vertical profile of water vapor concentration. This method assumes horizontal homogeneity within and around the observation footprint, and therefore, horizontal advection becomes negligible. However, when measuring evaporation near a sharp horizontal transition, the contributions of horizontal advection to the total evaporation budget need to be taken into account, even when the observed flux footprint is mostly over the water (Kenny et al. 2017). The relative importance of advection depends on the distance from the shoreline and the height of the measurement point. The relative portion of advection out of the total observed evaporation budget is generally larger when the instruments are higher (further from the surface) and closer to the coast. Although some authors have suggested the estimation of advection from multiple-elevation measurements at a single tower (Lee 1998;Leuning et al. 2008;Novick et al. 2014), a more accurate approach to calculate advection in heterogeneous surfaces is to use a multi-tower setup where two or more towers measure fluxes and vertical concentration profiles of water vapor, or any other scalar whose flux we want to measure (Aubinet et al. 2010;Feigenwinter et al. 2008;Staebler and Fitzjarrald 2004;Sun et al. 2007).
CO 2 fluxes in coastal waters are linked to changes in the solubility of aqueous CO 2 (Woolf et al. 2016), which are the result of changes in net ecosystem exchange (NEE) and changes in the carbonate system in the water (Longhini, Souza, and Silva 2015). NEE depends on the balance between heterotrophic respiration and photosynthesis, mostly attributed to coral/ zooxanthellae symbiosis and free-living algae (Small and Adey 2001). Changes in the carbonate system depend on the chemistry of the water, wave energy, and reef depth (Falter et al. 2013) and on the biogenic precipitation of calcium by corals (Jones, Ridgwell, and Hendy 2015;Small and Adey 2001). Although the biophysical characteristics of the water are important, CO 2 fluxes are also the result of atmospheric conditions that affect the water-air CO 2 exchange, in particular wind speed but also turbulence at the air-water interface and boundary layer stability (Wanninkhof 1992). Few studies have studied CO 2 fluxes under highly contrasting air-sea temperature conditions, such as those found in the Red Sea and any other desert-surrounded seas and coastal areas.
The Gulf of Aqaba has received a lot of attention due to its ecological importance and high evaporation rates. Estimates of evaporation in the Gulf of Aqaba range from 10 [mm/day] (Assaf and Kessler 1976) to more recent conservative estimates of 4.6 [mm/day] (Biton and Gildor 2011) and 5.1 [mm/day] (Ben-Sasson, Brenner, and Paldor 2009). Among the reasons for the high level of uncertainty are the seasonal variability of evaporation as well as differences in the methods used to provide these estimates. Another reason is that all previous studies did not measure evaporation directly but instead estimated it indirectly from atmospheric or water conditions. In terms of CO 2 flux, the observations are also limited. Barnes and Lazar (1993) used a metabolic model to predict the carbon uptake rates based on observations of the water conditions. Larkum, Koch, and Kühl (2003) performed laboratory measurements of carbon fluxes from corals collected in Eilat. However, we do not have any knowledge of the previous studies that directly measured CO 2 fluxes in the Gulf of Aqaba.
In this study, we used the EC approach to measure evaporation rates over land and sea and advection between them in the coastal waters of the Gulf of Aqaba (Figure 1(a)), as well as carbon fluxes between the coastal waters and the air. We used a flux tower over the water, and an additional tower over the coast with vertical profile measurements of humidity, temperature, and wind speed and direction. We measured additional environmental conditions, such as water temperature, tide height, precipitation, sensible heat flux, net radiation, turbulence, and atmospheric stability. We used this data set to determine the sensitivity of carbon dioxide to different environmental and meteorological drivers.

Study site
The study was conducted at the Interuniversity Institute for Marine Research (IUI), in Eilat, Israel, at the northern west shore of the Gulf of Aqaba (29°3 0.12ʹ N, 34°55.04ʹ E, 0 m elevation; Figure 1(a)). Using IUI data from 2007 to 2014, the site is characterized by a mean, minimum, and maximum annual air temperatures of 25.48, 6.6, and 45.0°C. The minimum occurs around the end of December, whereas the maximum occurs around late July. The water temperature has annual mean, minimum, and maximum values of 24.1, 20.6, and 32.0°C, respectively. With respect to the surface water dynamics, the site is characterized by three seasons throughout the year: sea surface cooling (October-February, corresponding to winter on land), sea surface warming and spring bloom (March-May), and the season of stable stratification and oligotrophic conditions (June-September, corresponding to summer on land). From the 2007-2014 data set, the April-May period has mean, minimum, and maximum air temperatures values of 25.8, 13.2, and 41.9°C, respectively, whereas the same values for the June-August period were 32.1, 22.1, and 45.0°C, respectively. There was no precipitation at the site during the observation period, and throughout the 5 years prior to it. The main orientation of the gulf is about 17°N, which more or less coincides with the predominantly north-north-east winds in the gulf (Afargan and Gildor 2015). The angle of the coastline at the study site (50°N), however, is shallower than the main angle of the gulf (Figure 1), which results in a large exchange of land-to-sea breezes when the wind blows predominantly from the north. During about 6% of the year, the Gulf experiences southerly winds, occurring mainly during storms in the winter and spring, whereas westerly and easterly winds are rare and recorded for about 2.5% and 6% of the year, respectively (Afargan and Gildor 2015). During our campaign, 3.3% of the time (roughly 8 days) was characterized by southerly winds and stormy conditions. Because conditions during these southerly wind events are very different from the regular season, we separated these events from the analysis of April-May and June-August data. The spatial variability of wind speed across the gulf is high, whereas the spatial variability of air temperature, water temperature, , and the advective coordinate system perpendicular to the mean angle of the coast (azimuth = 320°). (c) Scheme of the instruments setup: to the left, the land station with instruments at 4.2 and 6.4 m above the ground. To the right, the instruments located in a tower on top of the sea about 30 m from the shore. The instruments were located at 5.3 and 7.6 m above sea level at low tide. The permanent IUI station was located at 10 m above sea level at low tide. Ground level is 1.6 m higher than sea level at low tide. The watermark line delimits the western side of the control volume for calculation of fluxes. The magnitude of the variables at this point was calculated from an interpolation between the two towers. and relative humidity across the gulf is relatively small (Afargan and Gildor 2015).

Data collection
Data was collected for a period of 5 months, from 1 April to 27 August of 2009. Our campaign took place during the sea-surface warming, and the stable stratified periods, and out observations were sub-divided into two periods: April-May and June-August. 1 June was the first day where minimum nighttime air temperature surpassed the temperature of the water surface, and thus was chosen to separate the two seasons. Three flux/meteorological stations provided the data for this project: (1) one EC system and meteorological station located on a pier, protruding 32.8 m into the sea from the mid-tide shoreline (the "sea station"); (2) an EC system located on the beach, 25.3 m from the mid-tide shoreline (the "land station"); and (3) a meteorological station run by IUI located on the pier next to the sea station ("IUI meteorological station"). The sea and land stations collected meteorological data at two elevations of approximately 4 and 6 m above ground/high-water (Figure 1(c)). The IUI meteorological station provided observations at an elevation of approximately 10 m above the water level. At both the 6 and 4 m levels of the sea and land stations, we measured air temperature and relative humidity once per minute. We measured three-dimensional (3D) wind velocities at 10 Hz at the same locations. The 10 Hz CO 2 and water vapor concentration observations were available only at the sea station, except for 10 days when the land station also included these observations (see Table 1 for details of the instrumentation used). A net radiometer was installed at the 6 m level of the sea station. The IUI meteorological station provided observations every 10 min of wind speed and direction, air temperature, humidity, atmospheric pressure, as well as near-surface water temperature and tide height. The data from the IUI meteorological station are available through the Interuniversity Institute of Eilat (IUI). This data set (with observations since 2006) can be accessed through the Israel National Monitoring Program at the Gulf of Eilat (http://www.meteo-tech.co.il/eilat-yam/eilat_about_ en.asp). The instrumentation and observed variables at each station are listed in Table 1.

Despiking
We used typical quality control procedures for EC data analysis. A despiking process eliminated outliers using rules-based processing. We first filtered values according to acceptable lower and upper boundaries for each variable. We then removed all the data points whose distance to a half-hourly mean was higher than six standard deviations from the mean of a 2-minute window. A detailed description of the despiking and quality control approaches used in this study can be found in Morin et al. (2014b).

Flux processing
The calculation of water vapor flux and carbon flux through the EC technique is well documented in the literature. In this study, the flux data were processed using the approach described in Morin et al. (2014b). In summary, measurements of wind were subjected to a 3D rotation assuming the half-hourly mean vertical and crosswind components of the wind vector are zero (Finnigan et al. 2003). Given that the wind in the bottom few meters near the surface flows parallel to the surface, this rotation, in effect, sets our coordinate system parallel to the advective flow along the sea-coast surface. We assigned an advective direction of 40°as dictated by the mean angle of the coastline in our study zone (see Figure 1(b)). We also applied a correction for the temperature measurements from the sonic anemometer due to changes in the air density (Kaimal and Gaynor 1991). CO 2 flux and water vapor flux were subjected to a WPL (Webb, Pearman, and Leuning) correction (Webb, Pearman, and Leuning 1980) to account for changes in the density of dry air and water vapor and its effects on the final flux calculation. All flux and meteorological data were averaged to half hourly time steps.
Eddy-covariance data filtering and footprint model We averaged fluxes into half-hour intervals. A halfhour was rejected if more than 50% of the 10 Hz data within the 30-minute period were missing or eliminated by quality control. The half-hourly data were then subjected to a frictional velocity (u*) filter following the protocol of Reichstein et al. (2005). The minimum value allowed for a u* threshold was 0.15 m/s. Overall, from the initial data pool, 19.8% and 6.9% of the observations for the sea and the land, respectively, did not exceed this threshold and were excluded due to insufficient turbulent mixing. As our observations were near the shore, the flux observation footprints in the sea station typically included some degree of fluxes originating from the land. We used the two-dimensional expansion of the Hsieh, Katul, and Chi (2000) footprint model (Detto et al. 2006) to estimate the source area of the flux detected by the tower at each half hour. Observations with less than 70% of the sea footprint ( Figure 2) were removed from the data analysis. In the footprint calculations, we considered the change in height of the tower with respect to the sea due to changes in tide height.

Vertical interpolation of mean wind
We created vertical profiles of wind speed along the advective axis at both the sea and the land stations. Table 2 shows important markers along the vertical profiles of the sea and the land. The height of the markers of the sea was dependent on tide height, whereas the heights in the land were fixed. For simplicity, the horizontal vertical wind speed rotated along the advective axis will be referred as wind speed.
To calculate mean wind speed between upper and lower stations (see Table 2), mean wind speed was linearly interpolated between lower and upper stations, and extrapolated toward the water/land surface using a logarithmic velocity profile. The land extrapolation followed the equation: where U land z and U land low are the horizontal wind speeds in the advective axis at height z land above land and at the height of the lower station, respectively; z land is the height at which the wind speed is calculated; z land low is the height of the lower land station and is equal to 4.2 m (Table 2), and Z 0 land is the roughness length over land. In the case of the land Z 0 land was assumed constant, and calculated following Nakai et al. (2008) and Maurer et al. (2015), under stable conditions and land-only footprint.
where z land upper is the height of the upper station with respect to the land (6.4 m from Table 2) where the measurements of wind were taken. u s and u s * are the non-rotated horizontal wind and turbulence values under stable conditions, respectively.
For calculations of the vertical profile of wind speed in the sea, we had to incorporate the effect of changes in water level (WL). We expressed the wind speed at a height z sea above water level as Figure 2. Examples of two contrasting footprints on (a) 4 April and (b) 17 April. The colored areas represent an integrated percentage of the probability of a signal coming from a given source area. If the area of the footprint intercepting the sea was higher than 70% of the total area of the footprint, that half hour is accepted (case a). On the other hand, if the footprint area coming from the sea is lower than 70%, the half-hour is rejected (case b).
where U sea z is the wind speed at height z sea above the water level, z sea low is the height of the sea low station. WL ranges from 0 when the tide height is lowest to 1.31 [m] when the tide is highest. Z 0 sea is the roughness length for the sea, which was calculated according to the relation of Zilitinkevich, Perov, and King (2002): where c 1 = 0.1 and c 2 = 32, specific coefficients for coastal waters, v is the kinematic viscosity of water in [m 2 /s], and g is gravitational acceleration in [m/s 2 ].

Flux components using the 3D mass balance approach
We calculated the total evaporative flux using a flux divergence approach, discretized along a control volume whose axes are determined by the orientation of the coast, the water section from the sea station to the shoreline and the height of the towers (Figure 1(c)). Wind speed and water vapor concentration at the shoreline vertical were linearly interpolated from the land and sea stations. Within our control volume, q is the water vapor concentration and u, v, and w are the wind velocity components in the x (across shore), y (parallel to shore), and z (vertical) axes, respectively (Figure 1(c)). Δx is the distance between the sea station and the shoreline along the advective axis and is equal to 30.2 m. The height of the control volume, Δz, was equal to the height of the upper tower at the sea. The width of the volume, Δy, was equal to 1 m as we assume homogeneity along the coast. Following Sun et al. (2007), we formulated the mass balance equation for water vapor concentration to calculate total evaporation from the water surface: In Equation (5), U land z is the source or total mass of water vapor entering or leaving our control volume system. Overbars denote averages and the primes denote fluctuating quantities. Using Einstein notation, Equation (5) where the first term is storage, the second term is the advection and the third term is the flux divergence along the coordinate axis j = 1, 2, or 3 (or at the x, y, and z directions, across the shore, along shore, and up, respectively). Applying the inverse of the product rule of derivatives, we can transform the water vapor advection term in Equation (6) from its familiar form to an equivalent term composed of the difference of two surface flux terms: Discretization of source components of evaporation To discretize the continuous function in Equation (7) and calculate the different source components of evaporation, we used Gauss divergence theorem and solved for the source term S: where the operator _ Δ j ϕ ð Þ can be expressed as where subscript j is the index of the coordinate axis, superscript (i) is the downwind spatial coordinate along the j axis, and superscript (i − 1) the upwind spatial coordinate ΔA j is the aperture area of the control volume perpendicular to the j axis. The term ΔxΔyΔz is the volume of the control volume and composed of the products of the vertex lengths at the x, y, and z (1, 2, and 3) axes directions, u j is the wind velocity component in the j direction, corresponding to the u, v, and w wind components at the x, y, and z directions, respectively.
Storage. Taking the change in storage flux term from Equation (7) and defining q avg as the average of water vapor concentrations at the sea station and at the The markers for the sea change according to water level (WL). WL ranges from zero, when the tide is at its highest recorded level, to 1.31 m when it is at its lowest recorded level.
mean vertical line on the shoreline (which was interpolated from the land station), we used the central difference method to define storage as 2Δt (10) Advection. The advection component of Equation Because we assume that w = 0, there is no advection in the w (j = 3) component so that _ Δ 3 ϕ ð Þ ¼ 0. We also assume homogeneity along the shore so that _ Δ 2 φ ð Þ ¼ 0. The advection term is therefore composed of only the u component (across the shore): Substituting the aperture area in the x axis in Equation (12), and discretizing to the limits of the control volume across the x axis: sea station and the shoreline: mean water mark (wm): By canceling Δy, and discretizing the height, z, to an arbitrary number of vertical sections with a vertical width that varies with tide height (in this study we used 32 sections of~0.2 m from the water level to the level of the upper station in the sea), we obtain where Δz is variable through time and equal to the height of the higher level of the sea tower above water level divided by the number of vertical sections (32).
Flux divergence. From Equation (8) we obtain that the flux divergence in the three dimensions is equal to Because there is homogeneity in the y axis, we assume the flux divergence along this axis to be zero. In the horizontal axis, the horizontal flux divergence (HFD) is given by Due to technical problems that limited the time during which the gas analyzer at the land station was operational, the HFD was calculated only for 2 weeks in April-May and 1 week in June-August. During these times, we were able to determine that the contributions of HFD to the mass balance were small. We gap-filled the missing values of HFD using an artificial neural network (ANN), which satisfactorily predicted HFD based on continuous measurements of water temperature, air temperature, turbulence, and wind speed (R 2 = 0.38, p < 0.01).
The vertical flux divergence was equal to the turbulent flux measured by the EC tower at the tower height (E ec ). Because q'w' at the water surface is zero, E ec can be derived from the vertical flux divergence and expressed as Calculation of evaporation We define the total source of water within the control volume, S, as the product of the evaporation rate per unit water surface area, E v , and the area of the water surface at the bottom of the control volume, A w : As described earlier, E v can be defined as the sum of Storage, Adv x , HFD, and E ec . Multiplying by the volume of the control box, which is dependent on tide height, and dividing by the wet area of the ECOSYSTEM HEALTH AND SUSTAINABILITY control box, returns the evaporation in units of mass of water per unit area per unit time:

Modeling evaporation with conventional bulk formulations
We compared our observations of evaporation for the Gulf of Aqaba to predictions of evaporation from two commonly used empirical models. In the first case, we used a variation of the bulk aerodynamic model developed by Kondo (1975), which has been broadly used for sea surfaces and applied (Berman, Paldor, and Brenner 2000) and recommended for the Gulf of Aqaba (Ben-Sasson, Brenner, and Paldor 2009). In this formulation, evaporation is calculated as where ρ is the density of water vapor, L is the latent heat of vaporization, Q s and Q upper are the specific humidity at the sea surface and at the upper sea station, respectively, and u is the wind speed at the upper measurement level over the sea. C E is a bulk transfer coefficient for water vapor, which in the presence of buoyancy effects on transfers (as is the case for the Gulf of Aqaba), can be modeled according to the formulation of Kondo (1975): The higher the roughness sublayer resistance, log (z/z 0 ), the lower the bulk coefficient, and the lower the flux would be. The last term in Equation (21) represents the air-water interface resistance for the transfer of water vapor. The parameterization we used was determined by Kondo (1975) for rough sea surfaces with frequent high values of the critical roughness Reynolds number, u × h p /v, where v = 1.98 × 10 −5 is the kinematic viscosity of water and h p = 15z 0 is a roughness parameter. The conditions in our measurement location were suitable for this parameterization (data not shown). The shape parameters a and b were optimized for our study site using Markov-Chain Monte Carlo Metropolis-Hastings (MCMC-MH) algorithm (Haario et al. 2006;Haario, Saksman, and Tamminen 2001).
Additionally, we used the formulation of Kara, Rochford, and Hurlburt (2000) for modeling evaporation. The Kara formula is frequently used in global models such as the Hybrid Coordinate Ocean Model. This formulation differs from Kondo (1975) in its approach for calculating the bulk transfer coefficient, C E , to be used in Equation (20): where T water is the water temperature and T sea air is the air temperature at the upper station over the sea.
Parameters a through f were optimized for our study site using an MCMC-MH algorithm.
We evaluated the accuracy of each model using the Nasch-Sutcliffe efficiency coefficient, C: where Emo n indicates simulated values for the model of interest. The coefficient C returns values between negative infinity and 1. C = 1 represents a perfect match between measured and modeled data. More negative values indicate less accuracy whereas more positive and closer to one values indicate higher accuracy. C = 0 indicates that the model predictions are as accurate as the mean of the data.

Gap-filling and quality control
There were two types of gap-filling methods used in the study: (1) linear interpolations in time and space for environmental drivers and (2) ANN models for gap filling of fluxes. Linear interpolations in space were used to relate observations from one or two instruments in the same vertical profile to the value of missing observations of another instrument. This was applied to wind speed, air temperature (T land air ), and/or humidity (RH). Gap filling through this method was not applied to any other variable involving high-frequency fluctuations, such as u* and fluxes. Gap filling of missing wind direction was done through substitution from the nearest station. The percentage of data that was gap-filled in this way (either substituted or interpolated) was less than 20%. When no other point in the vertical profile was available, we used a bi-linear periodic function to model gaps in the variables showing a diurnal cycle, i.e. temperature, humidity, and net radiation. This substitution consisted of two linear interpolations: The first one is an interpolation between the previous and next available data values of the variable in consideration. The second one is an interpolation between values from contiguous days at the same time of day as the point needed to be gap-filled. A full description of the bi-linear function can be found in Morin et al. (2014b). Less than 11% of the data was gap-filled using this approach. The original and gapfilled data can be accessed through the Supplemental data accompanying this article.

Artificial neural network for H 2 O and CO 2 fluxes
We created ANN models for E ec , Adv x , and CO 2 fluxes. The ANN created empirical predictive models of carbon and water flux using environmental drivers causally related to the fluxes (Moffat et al. 2007;Papale and Valentini 2003). The results of the ANN had two purposes: (1) gap-fill E ec , Adv x , and CO 2 flux observations and (2) determine the most important environmental variables in the prediction of CO 2 flux. The latter approach was not applied to E ec and Adv x , since we wanted to compare the observations of H 2 O flux to common empirical formulae as described previously. As such, the environmental drivers for E ec and Adv x were chosen based on the significance of the pairwise relationship with the environmental variables that are assumed important according to the physical models. The drivers for CO 2 flux were chosen using a stepwise linear model, combined with the Akaike Information Criteria (AIC) to find a statistically justified parsimonious number of variables for the model (Morin et al. 2014a).
The architecture of all the ANNs consisted of the input layer containing predictors scaled from −1 to 1, and 2 hidden layers of 12 and 5 nodes connected through a hyperbolic tangent sigmoid (tansig) function. The 5-node layer passed through one more tansig transfer function to produce the output layer. From the input data set, a random 50% was used as a training set for the initial parameterization, another non-overlapping 25% was used as a testing set and the remaining was used as a validation set. The ANN created up to 100 models with a particular training set. The validation set was used internally to choose an appropriate model from different competing approaches. If the model produced predictions with R 2 > 0.75 on the testing data set, the model was selected, unless no satisfactory R 2 was reached after 5000 attempts, in which case we would accept a model with an R 2 higher than the mean of all R 2 + 0.1. After reaching a model with a sufficiently high score, one iteration of the model was completed. We created 100 iterations and averaged the best 10% iterations based on their R 2 value to produce the final model for a given set of environmental conditions. Environmental drivers of CO 2 flux CO 2 fluxes were calculated from EC measurements whose footprint came mostly (>70%) from the sea. Unlike water vapor, we did not have multiple probes to measure CO 2 concentrations on the land station, and therefore, calculation of CO 2 advection was not possible. The CO 2 flux we report here is only the vertical turbulent flux. We used the observed fluxes and environmental conditions to study the most important environmental drivers of CO 2 fluxes from the coastal waters, and construct an empirical model for these fluxes. To select the environmental variables to be considered in the model, we used the stepwise linear model, combined with the AIC to find a statistically justified multivariate regression model for CO 2 flux (Morin et al. 2014a). The selected variables were then used to feed the ANN model, which provided a better prediction of carbon fluxes based on the environmental drivers selected by the stepwise linear model. We used this modeled carbon flux time series to gap-fill the observations in periods when carbon flux observations were not available, due to instrument malfunction, or because the footprint was not over the water, or u* was low.

Uncertainty of flux estimates
We evaluated the uncertainty of our flux estimations coming from two sources: the gap filling of environmental drivers and (2) the uncertainty of the ANN model. For the environmental drivers, a bootstrapping technique was performed, creating 1000 synthetic data sets of all the drivers that were gap filled. The synthetic data sets were created by adding 1000 random Gaussiandistributed errors to the modeled data set. The standard deviation of these random errors was based on the statistics of model vs. observation using the same gapfilling linear model that was used to gap fill to simulated existing observations. When observations were available, we did not create synthetic data points. We ran an ANN model with 100 iterations for each of the 1000 synthetic data sets. From each ANN, we selected the best 10 models based on the R 2 for a total of 10,000 flux data sets. From this data set, we obtained the 2.5 and the 97.5 quantiles to obtain significant confidence bounds for our predictions. For daily averaged data, we provided two different types of bounds, the rough bounds, as described above, and the scored bounds. The scored bounds were created with the intention of penalizing days with more frequent gap filling of environmental drivers. To obtain scored bounds, the rough bounds were multiplied by a score from 0.5 to 1 that represented the quality of the data for that half hour. If the flux observation was available, no score was applied, as the observation does not carry gap-filling uncertainty. It was 0.5 if all the environmental drivers used to gap fill the flux estimates were product of observations, representing uncertainty only due to ANN but not due to other sources. It was 1 if all the environmental drivers were previously gap filled by one of the linear interpolations approaches in space or time. It varied between 0.5 and 1 as a proportion of the number and variables and half hours during the day for which gap filling was applied to environmental drivers.

Statistical analysis
We used paired two-tailed t-tests to compare averages of meteorological variables between the April-May and June-Aug seasons. For comparisons among the three seasons, we used analysis of variance (ANOVA). We used linear regression to compare observed fluxes to our models, to select variables as drivers for ANN, and to estimate the error in the gap-filled environmental drivers. All data processing and statistical tests were performed in Matlab R2017a. All data used in this study is available in supplementary data.

Seasonal meteorological conditions
The mean air temperature (T air ) for the year of study as measured by the IUI station was 25.33°C, which was only 0.15°C below the historic mean for the years 2007-2014. There was no precipitation during the measurement period. The mean T air and mean water temperature (T water ) were lower in April-May than in June-August by 6.3°C and 3.3°C, respectively ( Figure 3). Relative humidity (RH) was significantly higher in April-May than in June-August and significantly lower during southerly winds (Table 3). The difference in water vapor concentration in the air between the sea and the land (Q diff ) was significantly higher in June-August daytime than in the April-May daytime. On the other hand, nighttime Q diff was significantly lower in June-Aug than in April-May (Table 3). Southerly winds significantly decreased Q diff , making the air above the land sometimes more humid than the air above the sea (Table 3). Due to the wind direction during storms, only 28 usable half hours were available during storms when EC fluxes could be calculated with an appropriate footprint. We therefore excluded the storm periods from further analysis.
The mean wind direction for the period of study was 17°. Both periods, April-May and June-August, had very similar wind directions, with nighttime wind coming mostly from the north and daytime wind coming from the northeast. Sensible heat flux from the air to the sea was significantly higher during June-Aug than during April-May (Table 3). Frictional velocity (u*) was stronger during June-August than during April-May. The boundary layer stability coefficient (ζ) showed that the atmospheric surface layer above the sea was mildly stable during the entire period with stronger stability during June-August days (Table 3).

Water vapor fluxes
For the E ec ANN, we used air temperature (T air ), net radiation (R net ), vapor pressure deficit (VPD), water temperature (WT), turbulence (u*), and wind speed (u) as ANN predictor variables. We used 100 iterations for each ANN model. For advection, we used the same variables as above except WT and VPD and added wind direction. During four short periods (23 April-27 April, 3May-7 May, 21 June-25 June, 11 August−18 August), both land and sea flux stations were not operational due to power supply problems. These periods were excluded from the analysis.
Vertical turbulent flux of water vapor (E ec ) and horizontal advection (Adv x ) were modeled satisfactorily using the ANN (R 2 = 0.77 and 0.66, respectively). An average 50% of the observations in these correlations were used to train the ANN model. The neural network results were used to gap-fill missing Figure 3. Daily mean air temperature, water temperature, and water vapor concentration during the days of the study. The blue transparency around the air temperature shows the minimum and the maximum air temperature values for each day. The black marks along the lower bar demarcate the days with prevalence of southerly winds. Data were taken from the 10 m IUI station.
observations for calculation of the seasonal totals. Although the R 2 values were relatively high, the slopes of these relationships were always lower than 1, indicating that the neural network model tended to underestimate water vapor fluxes in the higher ranges. Approximately 73% and 61% of the data for April-May and June-August, respectively, were gapfilled using the ANN. Quality control statistics are listed in Table 4 for details about the effects of the filters on the data.
About 40% of the total evaporation was advected to shore and was not accounted by the tower observations of vertical turbulent flux, E ec (Figure 4). Change of storage within the observations' control volume (Storage) and horizontal turbulence flux divergence (HFD) were negligible in comparison to the high values of E ec and Adv x . Total evaporation rates (E v ) were similar during April-May and June-August but the partitioning between Adv x and E ec differed between the two seasons. The contributions of E ec relative to the total evaporation were significantly higher in April-May (7.2/11.4 [mm/day/mm/ day]) than in June-August (5.7/10.5 [mm/day/mm/ day]) (t-test, p < 0.001), while the contributions of Adv x to the total evaporation were significantly higher in June-August (5.1 [mm/day]) than in April-May (4.6 [mm/day]) (t-test, p < 0.001) ( Table 3). As seen in Figure 4, Adv x and E ec followed a similar temporal pattern (regression of daily means: R 2 = 0.44, p < 0.001, regression of half-hourly values: R 2 = 0.14, p < 0.001).
The lines depicted in Figure 4 are daily averages that include some level of gap filling. The uncertainty of our flux estimates associated with our gap filling technique, as well as the uncertainty of the ANN models is depicted in Figure 5. The areas around the lines represent 2.5% and 97.5% confidence bounds calculated based on the methods described in the heading 'Uncertainty of flux estimates'. For Ev, the 2.5% and 97.5% bounds were on average 4.0 [mm/day] below and 5.6 [mm/day] above our estimates, respectively. These rough bounds (outer area around each of the three lines in Figure 5) were calculated taking into account the average uncertainty of all the un-scored half-hours that were gap filled in a day. The compound uncertainty (inner area around the lines) was calculated using the scored uncertainty and considering the weight of the half-hours that were not gap filled into each daily average. Such bounds were on average 2.4 [mm/day] for the 2.5% quantile and 3.3 [mm/day] for the 97.5% quantile. The uncertainty in the E ec measurements was more constant due to regular u* and footprint filters, whereas the uncertainty of the advection estimates was much lower and sporadic than that of E ec , due to the much lower number of gaps in the advection data set. As a result, the uncertainty of our total evaporation estimates was highest when the advection measurements were gap filled, as in days 194-200 ( Figure 5). E ec and Adv x differed significantly between day and night. On average, E ec was 1.1 [mm/day] higher during the day than during the night (one way ANOVA, p < 0.001). As with E ec , Adv x was significantly higher during the day than during the night, with a difference of 0.6 [mm/day] (one-way ANOVA, p < 0.001) (Table 3). However, daytime advection was significantly higher in June-August than in April-May (one-way ANOVA, p < 0.001), whereas nighttime advection was significantly higher in April-May than in June-August (one-way ANOVA, p < 0.001) ( Table 3).
About 6.3% of the time, and particularly during southerly winds, half-hourly advection was negative, indicating transport of water vapor from the air above land into the sea. Negative advection was the result of negative gradients of water vapor concentration (negative Q diff ) accompanied by crosswind coming from the southwestern land ( Figure 6). When the wind blew from the north-northeast, Q diff was maintained positive, and when the wind blew from the south-southwest, Q diff was maintained negative. The build-up of this gradient was not instantaneous as it took about 1 [h] for Q diff to respond to the changes in wind direction brought by southerly winds (Figure 6). This means that southerly winds initially homogenize the boundary layer across the sea-land interface for negative advection to later take place when the southerly wind brings slightly higher humidity.

Modeling evaporation with standard formulae
Both Kondo and Kara models provided good estimates of the mean evaporation (E v ), calculated as the sum of E ec , Adv x , HFD, and Storage (see ratios of observed/predicted mean evaporation in Table 5). However, the models provided a poor representation of the half-hourly variation of evaporation as shown by the low Nash-Stutcliffe coefficient ( Table 5). The ratios between mean observed and mean calculated evaporation were closer to 1 in the Kondo model.  * Statistical difference between the April-May and June-August periods based on a t-test 2, p < 0.005. **The values in italics are ratios of observed/modeled means accompanied by the value of the Nasch-Sutcliffe coefficient (C).
The mean values are calculated with non-gap-filled observations as opposed to the total evaporation means shown in Table 3, which includes the ANN gap-filled points. However, the Nasch-Sutcliffe coefficient was closer to 1 when using the Kara model, which indicates that the Kara model provided a slightly better representation of half-hourly variation of evaporation. The Kondo model predicted evaporation as being significantly higher in June-August than in April-May, contrary to the observations. The Kara model, on the other hand, was accurate in predicting higher evaporation in April-May. Both models performed better in June-August than in April-May. They also performed better in the daytime than in the nighttime. The ratio of observed/predicted mean nighttime evaporation was close to 2 in both models, indicating a significant underestimation of the evaporation at , and their uncertainty associated with the gap-filling methods as well as the neural network model. The lower and upper limits are 2.5% and 97.5% confidence intervals. The outer area is the average uncertainty for all gap-filled half-hours on a day, whereas the inner area is the compound uncertainty that considers the weight of the half hours that were not gap-filled for that day. Figure 6. Half-hourly behavior of wind direction, total radiation, Q diff , and advection. Notice the negative values of advection occurring shortly after the wind starts to blow from the south during the southerly winds. night ( Table 5). The neural network produced the most accurate estimates of evaporation, reaching an overall R 2 of 0.77. The predictions of the ANN were more accurate in June-Aug (R 2 = 0.83) than in April-May (R 2 = 0.63).
Carbon fluxes and its environmental drivers CO 2 flux was modeled satisfactorily using the ANN (R 2 = 0.65). The fit was better in April-May (R 2 = 0.75) than in June-August (R 2 = 0.54) while the performance of the model in the daytime and nighttime was similar (R 2 = 0.64, 0.66, respectively). An average 50% of the observations in these correlations were used to train the ANN model. Figure 7 shows diurnal averages of CO 2 flux in the time of study. The areas around the line in Figure 7 represent the same confidence bounds as calculated for water fluxes in Figure 5. The 2.5% and 97.5% un-scored bounds were on average 2.0 [mm/day] below and 1.8 [mm/day] above our estimates, respectively, whereas the scored compound bounds were 1.0 [mm/day] and 0.9 [mm/day], respectively.
There was higher CO 2 uptake by the water in June-August than in April-May (Figure 7). The diurnal cycle of carbon flux had a peak assimilation rate during the night whereas during the day the net assimilation was closer to zero ( Figure 8). In the April-May period, the mean carbon flux reached positive mean values during the day (carbon coming out of the water), whereas during June-August, the carbon assimilation was always below zero, indicating continuous carbon uptake (Figure 8).
The stepwise linear model indicated that the most important drivers of CO 2 flux included humidity in the air (RH, VPD), and sea-air temperature difference (T air , T water , T diff , H). In the daytime, CO 2 flux had the strongest relationship with RH, E ec , and WT during April-May and with E ec , T diff , and wind speed during June-August. In the nighttime, VPD had the strongest relationship with CO 2 flux in April-May, followed by H, R net , and wind speed, whereas during June-August nighttime, the strongest correlation occurred with H, followed by VPD and E ec (Figure 9). The rest of the variables that had a significant relationship with CO 2 flux in the stepwise linear model and that were not rejected by the AIC analysis are shown in Figure 9. The daytime stepwise linear models reached coefficients of determination of 0.51 during April-May and 0.27 during June-August.

On the observations of evaporation
The observed evaporation in this study was higher than most estimates at the Gulf of Aqaba and in other places with similar conditions. We observed overall average evaporation rates of 10.6 [mm/day], over a five-month period which did not include the winter, whereas Biton andGildor (2011) andBen-Sasson, Brenner, andPaldor (2009)    (2011) estimated a rate for the entire Gulf using an oceanic model, whereas Ben-Sasson, Brenner, and Paldor (2009) calculated evaporation rates at our study site using an energy balance approach. Ben-Sasson, Brenner, and Paldor (2009) calculated a summer mean of 2.7 [mm/day] and a winter mean of 6.3 [mm/day]. This summer estimate is lower than our observations of 10.5 [mm/day for our short summer, the June-August season. Al-Subhi (2012) calculated evaporation rates of 5.41 [mm/day] for the period of June-August and 5.43 [mm/day] for the period of April through May. Other studies that measured evaporation using EC measurements in a similar environment also provided lower daily averages of evaporation. Assouline and Mahrer (1993) reported observing an annual rate of 4.1 [mm/day] in a small lake in Northern Israel, whereas for summer days in a small reservoir in Northern Israel, Tanny et al. (2008) found evaporation rates of 5.5 [mm/day]. These two EC observations, however, did not account for advection. The only previous estimate that is in good agreement with our estimate is the one provided by Assaf and Kessler (1976), who estimated the annual evaporation at the Gulf of Aqaba to be 3.65 [m/year], or 10 [mm/day] considering an average day throughout the year. The large differences in estimated evaporation rates demonstrate the complexities in measuring evaporation in desert-surrounded semi-enclosed basins. Johns et al. (2003) found that similar discrepancies exist for estimates of evaporation in the Persian Gulf and that these are in part caused by the uncertainties on the measurement of advection of water vapor.

Model estimates versus observations
Although the models' predictions were somewhat biased and with a moderate goodness of fit, the positive, non-zero Nasch-Sutcliffe coefficients indicated that the models provided good predictions of the mean of the observed data ( Table 5). Part of this result possibly occurred because both Kondo and Kara models were not able to capture the shortterm effects of advection in coastal waters, which are dependent on wind direction and sea-land water vapor gradients ( Figure 6). Assouline et al. (2008) found that advective conditions result in a large deviation of the ratio between latent and sensible heat from expected values and thus perturb the scaling of fluxes by models that use the Monin-Obukhov similarity theory, such as the Kondo model. Similarly, Assouline and Mahrer (1993) found that the highest differences between calculated and measured values occurred also during abnormal periods, during which the wind brings hot and dry air toward the water surface, which is generally the standard condition in our site.

Requirements for observations of coastal water evaporation
As our observations show, horizontal advection of water vapor plays a role in total evaporation rate from the coastal waters, particularly during shortterm maxima of evaporation rates. Higgins et al. (2013) found that to minimize the effect of advection in measurements of turbulent fluxes over sharp transitions, the ratio between the height of the tower (z sea ) and its downwind distance from the edge of the transition (Δx) must be less than 0.02. In our study, mean z sea = 7.03 [m] and mean Δx = 32.8 [m] (both distances are not fixed, as they depend on the wind direction and tide height, the mean horizontal distance was calculated in the direction perpendicular to the mid-tide coast line), and thus mean (z sea /Δx) = 0.21. Our measurements were therefore within the range at which the effect of advection must not be neglected, and indeed, our two-tower system observed a high contribution of advection as a component to the total water vapor flux. Following the recommendation by Higgins et al. (2013), it is possible to observe evaporation using a single tower that does not need to take into account the effect of advection. Nonetheless, such a tower would have to be installed about 150 [m] from the coast if it was at a height of 3 [m] above the mean tide. About 3 [m] is roughly the lowest position where it is technically feasible to install an EC system over the sea due to tide height variation, waves, and sea spray that will quickly destroy the instruments if in direct contact with sea water. Even at this low height, locating a system at such a distance from land is challenging because of restricted access to power, limited serviceability, and the challenges of accounting for the effects of waves on EC flux measurements from a floating platform. As a cost-effective approach to observe evaporation from coastal waters, we therefore suggest using a two-tower system that includes an EC-tower over sea and two vertical profiles of water vapor concentrations and wind speeds, and thus, can directly measure advection.
On the energy balance and contributions of advection to total evaporation Considering the energy balance above the water, we obtain that where R net is net radiation, LE is latent heat flux, H is sensible heat flux, and G is the water heat flux into the water (negative values) or released to the atmosphere (positive values). We assume that the change of heat storage in the air under the tower is negligible. Although we did not measure heat storage in the water, we did measure all the other three components of Equation (24), which allowed us to calculate heat storage as shown in Table 3. Average positive R net values indicated that the radiation reaching the surface was higher than the long-wave radiation emitted back to the atmosphere, except during the night, when the shortwave radiation from the sun becomes zero. H was negative in all seasons during both the day and the night, indicating heat flux from the air into the water for most of the time. The values of H were more negative during the summer, and during the day, when the air is much warmer than the water. LE was on average 313 [W/m 2 ], from which approximately a third (131 [W/m 2 ]) were a product of advection. The average magnitude of LE was much higher than the magnitude H, of −32.9 [W/m 2 ], producing Bowen ratios of ca. 0.1, typical of tropical/ Mediterranean oceans. Because of this high evaporation rates, and the high incoming radiation, the water became a net source of heat to the air, indicated by an average G of 98.3 [W/m 2 ]. The calculated heat flux from the water was much higher in the spring than in the summer. During the day, when the air is warmed by the adjacent desert, the water removed heat from the air. During the night, however, the water became a large source of heat toward the air (Table 3), which provided the energy and turbulent mixing necessary for the unusually high values of evaporation observed at night. Horizontal advection of water vapor across a wetto-dry transition has been studied extensively in irrigated agricultural fields surrounded by semi-arid zones, where it creates a particular problem since water is scarce and needs to be used efficiently. Several studies have demonstrated that the energy brought by advection of sensible heat in such agricultural fields makes up to 30%, 50%, and 28% of the total energy available for evapotranspiration (Brakke, Verma, and Rosenberg 1978;Hanks, Allen, and Gardner 1971;Prueger, Hipps, and Cooper 1996), respectively. Alfieri et al. (2012) found a contribution of advection of between 6-100 [W/m 2 ] in an irrigated cotton field surrounded by drier bare ground and grass. The differences were higher when the cotton was more developed. Higgins et al. (2013) found advection of 0-60 [W/m 2 ] in a lake surrounded by agricultural fields. These values are lower than the mean values we observed in this study (131 [W/m 2 ]) because in the montane environment where their measurements where conducted, the lake-land contrast of temperature and vapor fluxes were much smaller than in the Gulf of Aqaba.
For the case of water-land transitions in desertengulfed and G is the water heat fluxreas, Johns et al. (2003) found that advection of water vapor in the Straits of Hormuz was higher than the computed evaporation for the gulf, a result they attributed to poor sampling of evaporation in coastal areas and uncertainties in the advective flux. In that case, evaporation from the coast was assumed to be equal to the evaporation calculated in open waters far from the coast, which ignores the large advective flux caused by hot and dry sea breezes from the Arabian Peninsula into the gulf.
For the case of the Red Sea, Eshel and Heavens (2007) predicted that advection of water vapor should not contribute a significant input to the balance of water vapor in the boundary layer in the Northern Red Sea. They calculated that in the Northern Red Sea the lateral advective moisture flux is much smaller than the vertical flux and that atmospheric subsidence is the main mechanism through which water vapor in the boundary layer is balanced. Although they acknowledge the importance of the disparities of water vapor concentration between the sea and the land, they consider this discrepancy only a means by which subsidence occurs. We hypothesize that the main reason for the observed disagreement with Eshel and Heavens (2007) is that they calculated advection across the entire channel of the Red Sea (ca. 200 km of length), whereas the present study focused on observing the evaporative water balance in the coastal waters above the coral reef, roughly 30 m away from land. In addition, the mean values of cross-channel sea breezes determined by Eshel and Heavens (2007) (0.12 [m/s]) were much smaller than the values found in the present study (1.2 [m/s]). This difference, in turn, can be attributable to the orientation of the Gulf of Aqaba, which is almost perpendicular to the main axis of the Northern Red Sea that Eshel and Heavens (2007) studied (see Figure 1(a)).
The coral reef and its surrounding lagoon act as a sink of atmospheric CO 2 Our measurements indicate that the coral reefs in the Gulf of Aqaba are acting as a sink of atmospheric CO 2 . These results are different from those found by Gattuso et al. (1996), who observed that a reef flat in Moorea was a source of CO 2 to the atmosphere during the austral summer. However, they are consistent with studies on other reefs that have found that coral ecosystems act as a CO 2 sink (Cyronak et al. 2014). Barnes and Lazar (1993) studied CO 2 fluxes in the Gulf of Aqaba and found gross primary productivity of 11.8 [µmol/m 2 /s] with a production/respiration ratio of 1.6:1. The resulting net CO 2 fluxes are negative (indicating net uptake) and close to our average CO 2 fluxes of −1.3 [µmol/m 2 /s]. Barnes and Lazar (1993) suggested that the high productivity in Eilat is due to an extensive epilithic algal community rather than high calcification rates.
We found that evaporation as measured by the EC tower was negatively correlated with CO 2 flux. This meant that upward flux of water vapor was correlated with downward flux of CO 2 . This implies that the same environmental drivers that facilitate evaporation also determine the rates of CO 2 uptake. Using a stepwise model approach, we determined that the effective variables that predict daytime CO 2 fluxes during the daytime are the relative humidity, water temperature, the temperature difference between the sea and the air, and wind speed, whereas during the night sensible heat flux and vapor pressure deficit were more important. The fact that temperature and sensible heat were important drivers indicates that in our system the effects of buoyancy on CO 2 flux are important. The variables that drive evaporation such as vapor pressure deficit or relative humidity were also important in determining CO 2 flux. Wind speed was an important driver of CO 2 flux in both seasons, both during the daytime and the nighttime. Studies of CO 2 flux over the open ocean have shown that there is a positive relationship between wind speed and bulk fluxes of CO 2 (Wanninkhof 1992;Wanninkhof 2014). However, we found that this relationship was negative in both seasons ( Figure 9); higher wind speeds resulted in more negative CO 2 fluxes (more carbon entering the sea).
It is possible that additional environmental drivers are affecting the diffusivity of CO 2 in the water and therefore its rates of emission. Such drivers could be the dynamics of the exchange of warm water between the Red Sea and the Gulf (Cyronak et al. 2014) or changes in the salinity and temperature in the "skin temperature" of the sea, which have been considered by others as main drivers of CO 2 flux (Woolf et al. 2016). Other authors have reported that daily flux of CO 2 into the reef waters increases with daily irradiance (Gattuso et al. 1996) but we found that radiation is not an important driver of CO 2 flux in our observations. In addition, we found that the flux of CO 2 was higher during the night than during the day. Combination of these two facts suggests that photosynthesis is never light limited during the day in our observation site. Furthermore, at the intra-daily scale, other processes related to changes in the carbonate system of the water (Falter et al. 2013), and vertical mixing drive the instantaneous flux rates. At night when the water is colder, flux rates are higher due to higher diffusivity of CO 2 into the colder water (Longhini, Souza, and Silva 2015). However, more data on the physical, chemical, and biological characteristics of the water are necessary to better understand these processes and to offer better high-frequency estimates of CO 2 flux that can be coupled with the atmospheric conditions above the water.

Conclusions
We have directly measured evaporation in the coastal waters of the Gulf of Aqaba and calculated the contributions of advection of water vapor to total evaporation, showing that calculations of evaporation rates in coastal waters with commonly used formulas can approximate mean values of evaporation but have poor representativeness of the temporal variation of fluxes. This can result in overestimation or underestimation of seasonal evaporation in coastal waters, which can have implications to the modeling of the ecology of coastal marine organisms. We have also shown that wind direction can be a key driver of advection and thus must be taken into account in the calculations of evaporation. We show that the coral reef is a net sink of carbon during the period of our measurements. Typically, photosynthetic activity in the Gulf of Aqaba is higher in the winter and therefore our observations represent the minimal seasonal uptake rates. Our seasonal carbon uptake estimates of this coastal ecosystem are probably further underestimated by up to 40% due to lack of measurements of carbon advection, which should affect carbon flux estimates similarly to water vapor. The main environmental drivers of CO 2 flux were humidity, water temperature, sensible heat flux, and wind speed. Carbon flux was correlated to evaporation due to parallel environmental drivers affecting the vertical movement of water vapor and CO 2 . We conclude that considering an estimate of advection in coastal waters can improve evaporation and carbon flux models for coastal waters or land-water transition areas. This two-tower study approach can be applied to other desert-enclosed water bodies where land-sea gradients are strong.