Water temperature control on CO2 flux and evaporation over a subtropical seagrass meadow revealed by atmospheric eddy covariance

Subtropical seagrass meadows play a major role in the coastal carbon cycle, but the nature of air–water CO2 exchanges over these ecosystems is still poorly understood. The complex physical forcing of air–water exchange in coastal waters challenges our ability to quantify bulk exchanges of CO2 and water (evaporation), emphasizing the need for direct measurements. We describe the first direct measurements of evaporation and CO2 flux over a calcifying seagrass meadow near Bob Allen Keys, Florida. Over the 78‐d study, CO2 emissions were 36% greater during the day than at night, and the site was a net CO2 source to the atmosphere of 0.27 ± 0.17 μmol m−2 s−1 (x̅ ± standard deviation). A quarter (23%) of the diurnal variability in CO2 flux was caused by the effect of changing water temperature on gas solubility. Furthermore, evaporation rates were ~ 10 times greater than precipitation, causing a 14% increase in salinity, a potential precursor of seagrass die‐offs. Evaporation rates were not correlated with solar radiation, but instead with air–water temperature gradient and wind shear. We also confirm the role of convective forcing on night‐time enhancement and day‐time suppression of gas transfer. At this site, temperature trends are regulated by solar heating, combined with shallow water depth and relatively consistent air temperature. Our findings indicate that evaporation and air–water CO2 exchange over shallow, tropical, and subtropical seagrass ecosystems may be fundamentally different than in submerged vegetated environments elsewhere, in part due to the complex physical forcing of coastal air–sea gas transfer.

Shallow seagrasses of the tropics and subtropics are highly productive ecosystems, supporting economically and ecologically important food webs in surrounding regions. Seagrass meadows also cover and protect belowground organic carbon stocks which are significant to the global carbon cycle (Duarte et al. 2005;Fourqurean et al. 2012a;Fourqurean et al. 2012b). While these "blue carbon" reservoirs contain significant amounts of seagrass-derived carbon, nonseagrass autochthonous and autochthonous material often composes around 50% of this blue carbon stock (Kennedy et al. 2010;Röhr et al. 2018). During periods of high seagrass primary production and CO 2 consumption, the water column partial pressure of CO 2 (pCO 2 ) may fall below that of the atmosphere, causing a net transfer of atmospheric CO 2 into the water (Perez et al. 2018). Given the relative ease with which dissolved oxygen (DO) can be measured, ecosystem metabolic rates are often assessed with DO, which is converting into units of carbon assuming a specific O 2 : CO 2 ratio . Recent underwater O 2 eddy covariance studies have yielded well-resolved measurements of seagrass production (Long et al. 2015), showing that some seagrass systems are near metabolic balance (Attard et al. 2019) or can shift between autotrophy and heterotrophy . Benthic chamber-based DO studies have also shown net heterotrophy at some sites (Asmala et al. 2019). This use of DO as a proxy for CO 2 may be appropriate in siliciclastic and noncalcifying systems (Attard et al. 2019), where the ratio of O 2 : CO 2 consumption may be near 1: 1. However, in carbonate-dominated systems, changes in total alkalinity due to CaCO 3 formation and dissolution along with anaerobic alkalinity inputs will cause pCO 2 changes independent of O 2 . In net calcifying systems (common among tropical and subtropical seagrass), metabolic rates estimated using DO may be biased toward autotrophy, to an extent that is proportional to the ratio of net ecosystem production to calcification (Van Dam et al. 2019a). Likewise, DO-based methods may underestimate metabolism when carbonate minerals are being dissolved. Consequently, there is a clear need for direct measurements of air-water CO 2 exchanges in tropical and subtropical seagrasses. Such measurements, in conjunction with rigorous ecosystem metabolism studies, may begin to help resolving the question of whether these important ecosystems are sources or sinks of carbon to the atmosphere.
The subject of air-water CO 2 exchange over seagrasses has received much attention in recent years, owing in large part to the emerging interest in coastal vegetated habitats as carbon sinks (the concept of "blue carbon"; Nellemann et al. 2009), and the growing recognition of calcification as a potential confounding factor (Howard et al. 2018). Understanding airwater CO 2 exchange above seagrass meadows is key to their importance in global carbon budgets, yet few studies have measured CO 2 fluxes over seagrasses, and nearly all have taken place in temperate, tidally influenced systems (Polsenaere et al. 2012;Tokoro et al. 2014). The only prior measurements of CO 2 flux in Florida Bay (with floating domes) found clear regions of CO 2 emission and uptake throughout the bay (DuFore 2011), but lacked assessment of diel variability, thus presenting an incomplete picture of CO 2 exchanges and ecosystem metabolism. Because of the difficulty in directly measuring CO 2 flux, it is often estimated by combining measured pCO 2 with an empirical parameterization attempting to capture the physical forcing of turbulence at the air-water interface. However, the factors governing air-water gas exchange in these shallow, wind-exposed, and often fetch-limited systems are poorly understood and likely site-specific.
Underpinning this diverse forcing of gas transfer is the physical structure of the air-water interface, which is bounded on both sides by diffusive sublayers in the air and water. The exchange of gas across this boundary is rate-limited by diffusion through each of these sublayers, although for gasses of relatively low solubility (e.g., CO 2 , CH 4 ), the resistance is largely on the water side. While the size of the water-side diffusive sublayer is typically considered most sensitive to wind-driven mixing, other factors like convection and bottom-generated turbulence can be key when conditions are calm (Upstill-Goddard 2006;Wanninkhof et al. 2009). As a result, wellestablished wind-and current-based gas transfer parameterizations developed for the open ocean (Jiang et al. 2008;Wanninkhof 2014) and coastal waters (Ho et al. 2016) may not apply in seagrass meadows, where thermal forcing, surfactants, and other factors can affect variations in gas transfer independent of wind speed. This uncertain physical forcing of gas transfer has challenged prior efforts to constrain CO 2 fluxes over tropical and subtropical seagrasses, which make up the largest fraction of seagrass area globally (Green and Short 2003). Therefore, direct measurements of CO 2 flux in these "blue carbon" systems are highly desirable (Macreadie et al. 2019), especially when combined with a rigorous assessment of gas transfer.
In addition to investigating seagrass ecosystem productivity and carbon storage, air-water turbulent fluxes are also relevant to other factors affecting seagrass ecosystem health. For example, seagrasses in shallow and tropical waters are sensitive to extreme temperatures, with die-offs typically occurring during warm late summer months (Zieman et al. 1989;Robblee et al. 1991;Koch et al. 2007). The causes of these die-offs are varied and complex, but are made worse during conditions of high water temperature and hyper-salinity, which reduce O 2 solubility and allow H 2 S gas to accumulate to toxic levels (Borum et al. 2005). Therefore, a better understanding of the factors affecting salinity and temperature over seagrass meadows is relevant to the prospect of future seagrass die-offs. In areas with little freshwater inputs through surface or groundwater channels, like Florida Bay, seasonal and eventscale variations in salinity are driven by the relative rates of evaporation and precipitation (Nuttle et al. 2000;Swart and Price 2002;Lee et al. 2006), and by water mass advection, when salinity is spatially variable. Direct measurements of evaporation are not yet available for Florida Bay, and model estimates show methodological uncertainty (Lee et al. 2006;Price et al. 2007).
In this study, we used atmospheric Eddy Covariance to make the first direct measurements of CO 2 and H 2 O (evaporation) exchanges in Florida Bay, a well-studied seagrassdominated system. Our objective was to identify the dominant physical and biogeochemical processes governing turbulent fluxes over shallow submerged seagrasses. We describe a set of unique drivers that causes air-water CO 2 exchange in this system to differ from results in seagrass meadows elsewhere. Evaporation rates (latent heat fluxes) are also assessed and placed into the context of past and potential future seagrass die-offs. We also use heat transfer as a proxy for the gas transfer velocity, enabling the first preliminary assessment of the potential physical drivers of gas exchange in shallow subtropical seagrasses. Collectively, our findings support the role of temperature as a critical driver of air-water CO 2 and Van Dam et al.
Seagrass CO 2 flux latent heat exchange. These findings improve our understanding of air-sea CO 2 and H 2 O exchange in seagrasses, clarifying the role of tropical and subtropical seagrass ecosystems in broader regional and global carbon and water cycles.

Site description
Florida Bay is a large and shallow embayment stretching between the Florida Keys and the coastal Everglades. Fine carbonate mud throughout the bay is colonized by mixed seagrasses (mostly Thalassia testudinum, also Halodule wrighti, and Syringodium filiforme) and calcareous green macroalgae (Zieman et al. 1989). Primary producers are heavily phosphorus-limited to the northeast, but P delivered by tidal exchange with the Gulf of Mexico supports greater primary production in western Florida Bay (Fourqurean et al. 1992;Armitage and Fourqurean 2016) This study took place near the center of this P-limited productivity gradient (Fig. 1), near Bob Allen keys (5.027 N 80.681 W). Here, a shallow basin of 1-2 m depth is surrounded by mud banks which significantly restrict tidal exchange (Wang et al. 1994). Tidal amplitudes ($ 0.05 m) and water currents (< 2 cm s −1 ) are very low, and are driven by winds rather than lunar tides (Holmquist et al. 1989;Long et al. 2015).Water residence times in this region of Florida Bay are 6-12 months (Lee et al. 2006). Rates of net ecosystem metabolism (Long et al. 2015) and belowground organic carbon stocks are moderate for the Bay (Fourqurean et al. 2012b), despite phosphoruslimited seagrass primary production (Armitage and Fourqurean 2016). Underwater O 2 -based eddy covariance measurements made in 2012 (Long et al. 2015) showed that the gross primary productivity of the benthic community at Bob Allen Keys (151 ± 23 mmol O 2 m −2 d −1 ) was intermediate of sites to the east (Duck [68 ± 6 mmol O 2 m −2 d −1 ]) and west (Rabbit Key Basin [190 ± 27 mmol O 2

Eddy covariance and biometeorological measurements
Our eddy covariance system (Li-Cor, U.S.A.) was installed on an existing permanent data collection structure operated by the Everglades National Park at Bob Allen keys (water quality station BOBF1). This eddy covariance system consists of a rapidly responding (10 Hz) closed-path infrared gas analyzer (Li-7200RS), a sonic anemometer (Gill Windmaster Pro), and bio-meteorological sensors for net solar radiation (Kipp & Zonen NR Lite2), photosynthetically active radiation (PAR; Li-190R Quantum Sensor), and precipitation. Wind speed was corrected to a height of 10 m (U 10 ) assuming neutral conditions (Large and Pond 1981). Additional parameters such as water depth, temperature, and salinity were measured hourly by the Everglades National Park staff at the same location. Hourly average water temperature and salinity records were linearly interpolated to 30 min to match the measurement frequency of the other results. Eddy covariance instrumentation was operational from 22 March 2019 until 07 June 2019, for a total of 78 d. Air-water fluxes of CO 2 , latent heat, and sensible heat were determined from high-frequency (10 Hz) measurements of CO 2 , H 2 O and temperature, and were calculated at 30-min intervals using an automated routine in EddyPro (Li-Cor). Equations used for calculating latent heat, sensible heat, and CO 2 fluxes are described in Foken et al. (2012). Positive latent and sensible heat fluxes represent heat loss from the water. In order to more comprehensively assess the energy balance, we also estimated the water-column heat storage (J; W m −2 ) for each 30 min interval as: where ρ is the average water density for that time period (kg m −3 ), C P o is the specific heat of water (J kg −1 K −1 ), ΔT/Δt is the average rate of water temperature change (K s −1 ), and h is the water depth (m). While this approach necessarily ignores the heat flux into sediments and any lateral advection, information is lacking to estimate these fluxes. In order to assess the role of water-side convective forcing on CO 2 and latent heat fluxes, we calculated buoyancy flux (B) according to Czikowsky et al. (2018). First, the effective heat flux (Q eff ) was calculated as the downwelling net solar radiation (R n ) minus the latent (LE) and sensible (H) heat fluxes.
The buoyancy flux (B; m 2 s −3 ) was calculated as in Podgrajsek et al. (2015): where g (m s −2 ) is the gravitational acceleration, α (unitless) is the thermal expansion coefficient, and C pw (J kg −1 K −1 ) is the specific heat of water. When calculated in this manner, positive buoyancy fluxes indicate a loss of density from the water, which thereby becomes more buoyant. The result of this is greater stratification and an increasing stability of the sea surface microlayer, which can hamper rates of gas transfer. On the other hand, when buoyancy flux is negative, there is a gain in density at the water surface, enhancing convective overturn, and thereby gas transfer. We then converted all negative values of buoyancy flux (equation not valid for buoyancy fluxes > 0) into a water-side convective velocity scale (w *w , m s −1 ) according to the equation (Jeffery et al. 2007;Van Dam et al. 2019b): Data postprocessing and screening Postprocessing of high-frequency data in EddyPro involves a variety of steps including tilt correction, time-lag compensation, Webb-Pearman-Leuning corrections (Webb et al. 1980), as well as high/low pass filtering. In addition to these steps, the 30-min data were further screened by a variety of criteria intended to ensure that final data were representative of the appropriate flux footprint (i.e., the calculated footprint given wind speed and direction only contained the intended focus area, the submerged seagrass ecosystems), and to exclude times when conditions were nonstationary (i.e., statistics of flow are consistent through time), or when the IRGA optics appeared dirty. These conditions, the problem each criterion seeks to remedy, and the quantity of data lost to each criterium are shown in Table 1. The quality control procedure implemented by EddyPro (Mauder and Foken 2006) has become an internationally standard approach to postprocessing of eddy covariance data (see FLUXNET and AmeriFlux), and we have adopted the same approach here. Table 1. Criteria used to screen eddy covariance data, including brief justification and fraction of measurements that did not meet each step (see "Eddy covariance and biometeorological measurements" section for a description of screening variables). The QC code generated by EddyPro ranges from 0 (no suspected issues) to 2 (flux results highly suspect), and we took the conservative approach of only using flux results when the QC code was flagged as "0." This screening procedure was the most common cause of data loss, accounting for 33% of rejected flux data. Data failing the EddyPro QC threshold also often failed other criteria, which collectively accounted for the removal of 43% of records. It was also important to screen data when atmospheric conditions were excessively unstable. Lower atmospheric stability can be assessed quantitatively through Monin-Obukhov similarity theory, where the length scale L (m) is the height above the sea surface at which shear and buoyant forcing are equally important. When L is less than the measurement height (z), the turbulent processes measured there are influenced by buoyant, rather than shear forcing. Therefore, we can use the ratio of z/L as an indicator of the relative importance of shear and buoyant forcing (Smedman et al. 2007). When z/L is greater than 0, the boundary layer is stable, and vertical transfer is driven by wind-shear. Negative values indicate convective, unstable conditions. Hence, we screened data that fell below a z/L threshold of −0.55 (Czikowsky et al. 2018), as described in Table 1. The intent of this step was to remove time periods when turbulent exchanges were influenced by regions outside of the intended flux footprint, due to excessive convective forcing and minimal wind shear. Such conditions accounted for the rejection of 4% of records, were evenly distributed between day and night, and occurred exclusively when wind was below 5 m s −1 . Last, the size of the eddy covariance flux footprint increases with z/L, such that the footprint may become problematically large (extending outside of the intended focus area) when z/L exceeds 0.15 (Mørk et al. 2014). Maximum z/L during the study period was 0.146, suggesting that such nonstationary and footprint mismatch issues were unlikely. The cumulative effect of all of our data screening procedures can be seen graphically in Fig. 1, where the eddy covariance flux footprint initially intersected with the mangrove island to the west, before data screening. Following our screening procedure, the flux footprint was limited to the relatively homogeneous area just east (< 500 m) of the eddy covariance tower, with a benthic cover of patchy seagrass (Thalassia testudinum and Halodule wrightii) and calcareous green macroalgae in 1-2 m water ( Fig. 1, inset).

Physical drivers of CO 2 and latent heat fluxes
Air-water gas exchange can be expressed using a "bulk transfer" relationship where F gas is the gas flux across the sea surface (measured by eddy covariance in this study), ΔC is the concentration difference between the water and air, and k w is the gas transfer velocity: While a rigorous assessment of gas transfer requires waterside measurements of CO 2 concentration (not measured in this study), it may still be possible to estimate k w and investigate its drivers using air-water exchange of heat as a proxy for gas transfer. Prior studies have found positive relationships between the rate of air-water heat transfer (k g ) and that of airwater gas transfer (k w ), although k g often exceeds k w (Zappa et al. 2003). Given this general relationship between k g and k w , our direct measurements of air-water heat exchange (sensible + latent) can be used to estimate k w (Wanninkhof et al. 2009). As an extension of Eq. 5, heat transfer can be parameterized as: where ΔT is the water-air temperature gradient and F heat is the net heat flux (sensible + latent). Rearranging Eq. 6, we can calculate k g as follows: We can then convert k g into a gas transfer velocity (k w ) according to Schmidt (Sc) to Prandtl (Pr) scaling: Values of k w less than 0 cm h −1 were excluded from subsequent analyses due to heat fluxes being inconsistent with airwater temperature gradient. While gas transfer velocities above 50 cm h −1 have been observed under high wind speed and wave-breaking conditions in the open ocean (Smith et al. 2011) and during storms in Fjords (Mørk et al. 2016), such conditions did not occur in this study, causing us to exclude k w > 50 cm h −1 as outliers (Jonsson et al. 2008).

In situ measurements
In addition to the micrometeorological measurements described previously, a variety of hydrographic and biogeochemical sensors were deployed at the eddy covariance tower over varying intervals. First, an upward-facing Acoustic Doppler Current Profiler (ADCP) was deployed at the base of the tower, but given the shallow water depth (< 2 m), all ADCP bins were combined into a water-column average speed and direction. Additionally, a Datasonde capable of measuring temperature and salinity (EXO2, Yellow Springs International, U.S.A.) was deployed at the tower over various intervals, covering a total of 129 h. A direct comparison of our temperature and salinity measurements with those from the Everglades National Park showed a linear relationship with a slope of 1.0 ± 0.3 and an r 2 of > 0.96, and as such, we used our measurements to gap-fill when Everglades National Park data were absent or faulty. Unless otherwise noted, all error values are reported as mean ± standard deviation (SD).

Energy balance
Over the study period, this site was a net source of sensible and latent heat to the atmosphere, with an average sensible and latent heat fluxes of 15.1 ± 13.8 W m −2 (xc ± SD) and 100.5 ± 51.8 W m −2 , respectively (Fig. 2c). Latent and sensible heat fluxes were positive for 99.9% and 90% of measurements, respectively, and heat fluxes were dominated by latent heat flux, which was greater than the sensible flux 99% of the time. The Bowen ratio (sensible flux: latent flux) was therefore quite low, with an average of 0.15 ± 0.18 (xc ± SD), typical of shallow tropical waters (Polsenaere et al. 2013;Rey-sánchez et al. 2017). The energy budget exhibited very poor closure (slope < 0.2) under the simple assumption that net solar radiation is balanced only by the sum of sensible and latent heat fluxes (Fig. 2a). However, when the heat storage term (J) is included as an energy sink, the slope indicating energy closure is nearly 1: 1, indicating better energy balance closure (Fig. 2b). Water-column energy storage due to solar heating is therefore an important component of the energy budget, as clearly endorsed by the large diel signature in water temperature (Fig. 3b). The mean evaporation rate (3.6 ± 1.8 mm d −1 ; xc ± SD) was an order of magnitude greater than precipitation (0.32 ± 4.4 mm d −1 ; xc ± SD), in close agreement with prior estimates for the Florida Bay dry season (Price et al. 2007).
Further evidence for net evaporation can be seen in the salinity and temperature record, which shows a clear trend of increasing water temperature, rising from 20.34 C to 33.56 C and salinity, which rose $ 14% from 37.19 to 42.25 over the study period (Fig. 3).
Temporal trends in air-water CO 2 /latent heat fluxes and associated parameters A footprint analysis (Kljun et al. 2015) indicated that an area less than 100 m east of the eddy covariance tower provided the majority (83%) of CO 2 flux measurements (Figs. 1, 4a). While CO 2 flux was generally distributed evenly across wind directions, it was slightly greater in the window of 70-90 . In this narrow window, U 10 was also slightly higher ( Fig. 4b), offering a partial explanation for the increase in CO 2 flux. Nevertheless, the flux footprint for this data set was generally homogeneous and dominated by seagrass beds overlain by a $ 1.5 m water depth (Fig. 1). Mean U 10 was moderate, at 6.6 ± 1.8 m s −1 (xc ± SD), and the presence of the trade winds prevented any development of diel cycle in wind speed or direction. Maximum U 10 during the study period was 13.6 m s −1 . Over the study period, the majority of 30-min records of z/L (1631 out of 2588, or 63%) were within the unstable but very close to neutral (UVCN) zone (Fig. 5). These conditions of moderate atmospheric stability have been Van Dam et al. Seagrass CO 2 flux implicated in relative gas transfer enhancement (Smedman et al. 2007;Sahlée et al. 2008;Andersson et al. 2018), and we develop this concept in greater detail later in the "Discussion" section.
Water temperature was almost always greater than air temperature, and the average daily range in water temperature (2.8 C) was greater than the diel range in air temperature (2.0 C). As a result, the air-water temperature gradient (ΔT = T water − T air ) was almost always positive and increasing throughout the day, from a dawn minima to a dusk maxima, mirroring the trend in water temperature (Fig. 6d). The likely cause of this was the dominance of easterly "trade" winds. Air parcels reaching the site were therefore of a marine origin, with stable temperatures reflecting thermal equilibrium with ocean water, rather than a typical signature of daytime solar heating over land.
Over the study period, this site was a net source of CO 2 to the atmosphere, with an average CO 2 flux of 0.27 ± 0.17 μmol m −2 s −1 (xc ± SD). Daytime CO 2 emissions (0.29 ± 0.17 μmol m −2 s −1 ; xc ± SD) was 17% greater than the nighttime CO 2 release (0.25 ± 0.17 μmol m −2 s −1 ; xc ± SD), and the only observed period of net CO 2 uptake occurred between 02 and 05 May, following a storm that generated high winds, relatively low water temperature, and no measurable rain. This period of CO 2 uptake was limited only to nighttime hours and was offset by CO 2 emissions during the day. A short period of CO 2 uptake was also observed during the nighttime hours of 24-25 March, when water temperature was near the minimum for the study period. While the study site was affected by many other storm events and generated at times significant rainfall rates and/or high wind speeds, they had no clear impact on trends in CO 2 flux.
Diel variability in both CO 2 and latent heat flux were large, typically around 0.15 μmol m −2 s −1 and > 100 W m −2 , respectively (Fig. 5). Hourly variations in CO 2 flux were also out of phase with PAR, which peaked around 12:00 h (Fig. 6b), compared with the CO 2 flux maxima near 15:00 h. Maximum CO 2 flux was typically just before sunset, while minimum CO 2 flux was often near sunrise. The large diel variability in ΔT coincided with a strong diel excursion in latent heat flux (Fig. 6b), such that peak latent heat flux (18:00 h) occurred at  with points colored blue for u * < 0.2, and orange for u * > 0.2 (a). The difference between the patch at a footprint of $ 75 m, and the footprint patch between 100 and 400 m is caused by the footprint calculation in EddyPro, which switches between two different algorithms at a u * threshold of 0.2 m s −1 . FCO 2 as a function of wind direction, with points colored by U 10 (b).
the same time as maximum water temperature, after PAR (and net solar radiation, not shown) has already peaked and begun to decrease. Buoyancy flux exhibited a diel trend similar to that in PAR and net solar radiation (not shown) with a maximum gain in buoyancy during the middle of the day, and buoyancy loss (convection) at night (Fig. 6c).

Correlations between measured fluxes and physical drivers
While prior studies have found net solar radiation to be a good predictor of evaporation in Florida Bay (Price et al. 2007) and in the large and shallow lake Taihu, China (Lee et al. 2014) and Lake Erie, U.S.A. (Shao et al. 2015), we found weak correlations between latent heat flux and net solar radiation (R 2 < 0.1; root mean squared error = 50.1). Instead, latent heat flux was best predicted by an exponential function of waterside convection (w *w ) explaining 33% of the variability in latent heat flux. Likewise, significant positive correlations were found between latent heat flux and vapor pressure deficit, ΔT (Liu et al. 2009), and U 10 . The strength of the correlation between latent heat flux and w *w may be inflated due to the fact that latent heat flux itself is included in the calculation of w *w . The correlation between latent heat flux and air temperature was very poor (R 2 = 0.01), while that with water temperature was slightly better (R 2 = 0.07), further indicating that variations in ΔT were dominated by water temperature variations, rather than air temperature. The interactive effects of wind and temperature on latent heat fluxes can be clearly seen in Fig. 7c, where at a given wind speed, latent heat flux increases with increasing temperature. It was necessary to present absolute values of CO 2 flux (jCO 2 fluxj) in Fig. 7e-h because the occasional negative CO 2 flux measurements (CO 2 uptake) would have obscured the relationship between the various physical drivers and CO 2 flux. The strongest correlations for jCO 2 fluxj were with w *w (R 2 = 0.13) and U 10 (R 2 = 0.10), while less predictive power was provided by ΔT and vapor pressure deficit (R 2 < 0.10). A linear correlation between CO 2 and latent heat flux (R 2 = 0.37) was also evident, reinforcing the concept that latent heat and CO 2 fluxes shared common physical drivers. Last, ADCP water velocity was not correlated with fluxes of CO 2 or latent heat, indicating that bottom-driven turbulence is at most a secondary driver of turbulent fluxes (p < 0.05; R 2 = 0.013 for CO 2 flux and R 2 < 0.01 for latent heat flux). While there was a significant negative linear relationship between water depth and CO 2 flux (increasing CO 2 flux with decreasing water depth), this correlation was weak (r 2 < 0.1) perhaps due to the very narrow tidal range at this site. Together, these relationships suggest that wind-shear and convective forcing were the dominant drivers of turbulent fluxes, exceeding the impact of bottom-generated turbulence.

Variations in estimated gas transfer velocity (k w )
In this study, we used heat transfer as a proxy for gas transfer to make the first estimates of k w over a seagrass meadow. This k w compares well with three commonly used gas transfer models: Wanninkhof 2014 (k W14 ), Jiang et al. 2008 (k J08 ), and Ho et al. 2016 (k Ho16 ) when considered as an average over the study period (Fig. 8a), and also captures much of the same short-term variability (Fig. 8b). The inclusion of current speed and water depth into the model (only k Ho16 ) also yielded gas transfer velocities similar to our calculated k w , but only for the $ 2.5 d for which ADCP data were available (Fig. 8b). Together with the poor relationship between CO 2 flux and water velocity, this finding suggests that bottom-generated turbulence was not a significant driver of air-water gas exchange.
In fact, the average difference between our calculated k w and the models: k W14 , k J08 , or k Ho16 was not significantly different from zero, at −0.78 ± 9.0 cm h −1 (k w − k W14 ; xc ± SD), 0.55 ± 8.5 cm h −1 (k w − k J08 ), and −3.6 ± 9.3 cm h −1 (k w − k Ho16 ). This agreement between parameterized and calculated k w further supported wind as a driver of gas transfer in shallow waters. However, it was also apparent that these parameterizations tended to over-predict k w when wind-shear was high. This effect is evident in Fig. 7c,d, where the difference between calculated and modeled k (k w − k W14 and k w − k J08 ) becomes more negative with increasing u * . This overprediction was greatest when ΔT exceeds 2-3 C, indicating that these common models of gas transfer may overpredict k w , especially during the day as water heats relative to the air.

Discussion
Temperature-dependence of summer CO 2 flux in Florida Bay We documented persistently high temperature differences between air and water (ΔT) throughout the study period. The likely cause of this was intense solar heating of shallow water, combined with easterly "trade winds" delivering maritime air masses of consistent temperature. Coinciding with this diel variability in ΔT, we also observed a large diel range in latent heat flux and CO 2 flux, which were both greatest in the late afternoon and at a minimum near sunrise (Fig. 6b). Not only were diel trends in CO 2 flux large, they were of opposite phase of what has been observed in other, noncarbonate, seagrass systems (Polsenaere et al. 2012;Tokoro et al. 2018). We also observe a significant positive correlation between CO 2 flux and ΔT (Fig. 7e), although the strength of this correlation was weak relative to U 10 and w *w ( Table 2). The reason for this correlation between CO 2 flux and ΔT is unclear, but could be related to factors like the temperature-dependence of pCO 2 (Takahashi et al. 1993), or the temperature sensitivity of ecosystem respiration and calcification/dissolution. While robust Fig. 8. Violin plot (a) of mean k w (estimated from heat fluxes, at in situ Sc), compared with gas transfer models (k W14 , k J08 , k Ho16 ), showing all records from the study period. ADCP current data were only available for a short period, explaining the small sample size for the Ho16 parameterization which considers water velocity. Selected time-series for 2 weeks in early April 2019 (b). Figures (c-e) show the difference between estimated and parameterized k w as a function of the friction velocity (u * ), with points colored by ΔT. Table 2. Best fit models (linear or exponential) for latent heat flux (LE) and absolute value CO 2 flux (jFCO 2 j), for the set of drivers shown in Fig. 7. Vapor pressure deficit is abbreviated here as VPD. All equations presented have slopes significantly different from zero (α = 0.05). carbonate chemistry measurements are required to fully address this question, we can easily estimate a general thermal impact on pCO 2 variability. First, if we assume air-water CO 2 equilibrium, we can use the average atmospheric CO 2 concentration from this study (401.5 μatm) as an initial pCO 2 to which we can apply an isochemical temperature dependence of ∂ln(pCO 2 )/∂T = 0.0423 C −1 (Takahashi et al. 1993). This temperature sensitivity has an empirical relationship with salinity, ranging from 0.0332 at a salinity of 0-5 (Joesoef et al. 2015) to 0.0423 at a salinity of 36 (Takahashi et al. 1993).
Combining this with the average daily range in water temperature over the study period (2.31 C), we can estimate a temperature-driven pCO 2 increase of 10.28% [exp(0.0423 × (2.31 C)) = 1.1028], or an increase of 41.3 μatm from the initial pCO 2 of 401.5 μatm. Using an average gas transfer velocity (k w ) from this study (11.7 cm h −1 , or 0.117 m h −1 ) and CO 2 solubility (26.1 μmol m −3 μatm −1 ), this pCO 2 increase should enhance CO 2 flux by 26.1 × 41.3 × 0.117 = 126.1 μmol m −2 h −1 , or 0.035 μmol m −2 s −1 . This thermal effect can explain a similar magnitude decrease in CO 2 flux as water cools over the evening. The average diel range in CO 2 flux during the study was $ 0.15 μmol m −2 s −1 , meaning that this thermal pCO 2 effect can explain 0.035/0.15, or $ 23% of the typical range in CO 2 flux over the study period. While this thermal effect is clearly significant, a separate source of CO 2 for the remaining CO 2 emissions (during both day and night) must exist. Plausible sources are net ecosystem respiration and calcification, which should both increase during daytime warming. Florida Bay seagrasses are highly productive, with underwater eddy covariance studies (Long et al. 2015) showing large benthic O 2 production at Bob Allen (25 mmol O 2 m −2 d −1 , or 0.29 μmol O 2 m 2 s −1 ). However, seagrass productivity is only one component of net ecosystem metabolism, which is in turn just one factor contributing to net CO 2 production or consumption and ultimate exchange with the atmosphere. As a case in point, seagrass beds elsewhere in Florida Bay have been found to be net heterotrophic (Van Dam et al. 2019a), offering a partial explanation for prior observations of pCO 2 above (Millero et al. 2001) or near (Yates et al. 2007) equilibrium with the atmosphere in this area. Likewise, Berg et al. (2019) attributed a shift from net O 2 production during the spring to O 2 consumption in the summer (in a temperate Zostera marina meadow) to temperaturesensitive increases in ecosystem respiration during the warm summer months.
We also suspect a variety of other benthic processes as CO 2 sources/sinks affecting CO 2 flux, including calcification/carbonate dissolution, denitrification, and net sulfur biogeochemistry. First, while Florida Bay overall should be a net producer of carbonate minerals (Bosence 1989), only some regions are indeed net calcifying (Yates and Halley 2006;Turk et al. 2015), with other regions apparently net dissolving (Van Dam et al. 2019a). Assuming all CO 2 produced or consumed by calcification exchanges with the atmosphere, recent estimates of carbonate dissolution (Van Dam et al. 2019a) and precipitation (Turk et al. 2015) in this region can be used to explain air-water CO 2 fluxes of −0.14 to 0.57 μmol CO 2 m −2 s −1 , respectively. These estimates are well within range of our CO 2 flux observations, indicating that carbonate precipitation may be responsible for the net CO 2 emission at this seagrass meadow. Next, organic carbon respiration via nitrate reduction (denitrification, dissimilatory nitrate reduction to ammonium) is an important heterotrophic pathway in seagrass sediments. While the CO 2 produced via denitrification is likely orders of magnitude below measured CO 2 flux (Eyre and Ferguson 2002;Welsh et al. 2000), the alkalinity generated by net denitrification may be important to the carbonate equilibria of overlying water. Likewise, while rates of sulfate reduction are high in Florida Bay (Walter et al. 2007), CO 2 produced by this mechanism appears to be at rates far below our CO 2 flux measurements. Nevertheless, the reoxidation of sulfide in the root zone enhances carbonate dissolution (Ku et al. 1999), with uncertain impacts on watercolumn carbonate chemistry.
Given the recent interest in the carbon sequestration potential of seagrasses (Duarte et al. 2005Fourqurean et al. 2012a), and the growing recognition that CaCO 3 reactions must be included in "blue carbon" accounting (Mazarrasa et al. 2015;Macreadie et al. 2017;Howard et al. 2018), there is a clear need for studies combining direct measurements of CO 2 flux with determinations of carbonate precipitation/dissolution, as well as sulfate and nitrate reduction rates. The bulk of measured CO 2 emissions in this study cannot yet be attributed to a direct biogeochemical source, providing further motivation for such studies.

Physical drivers of the gas transfer velocity (k w )
In the present study, we estimated k w using heat as a proxy for gas transfer and found that it agreed reasonably well with k w parameterized by wind speed (k W14 , k J08 ) or a combination of wind, current speed, and water depth (k Ho16 ). On average, the difference between k W and k W14 was not significantly different from zero (i.e., k w −k W14 = −0.98 ± 8.6 cm h −1 [xc ± SD]). However, these parameterizations all tended to overpredict calculated k w especially at higher friction velocity (u * ) and ΔT (Fig. 8b,c). Prior estimates of k w using similar heat-transfer approaches have largely found this k w to instead exceed the directly measured gas transfer velocity, indicating that the apparent gas transfer suppression in this study is not an artifact. Other studies in stratified lakes (MacIntyre et al. 2010), rivers (Berg and Pace 2017), Everglades freshwater wetlands (Ho et al. 2018b), and Amazon floodplains (Polsenaere et al. 2013) have observed similar gas transfer suppression during periods of heating under light to moderate winds. It is likely that this gas transfer suppression was caused by thermal stratification, generated during the day by large buoyancy fluxes into the water (Fig. 6c). This thermal stratification will dissipate when buoyancy fluxes reverses sign after dusk. However, some time is required for stratification to be broken down by net heat fluxes as well as wind-and current-driven shear, eventually causing air and water temperature to approach equilibrium. Therefore, while we did observe a clear diel pattern in ΔT when plotted as a climatological hourly mean (Fig. 6d), there were many occasions where positive ΔT anomalies persisted through the evening and night hours (not shown). Accordingly, we observed that these periods of gas transfer suppression (k w < k W14 ) occurred both during the day and at night, when ΔT has not yet decreased. This effect is illustrated in Fig. 9c-f, where significant relationships exist between both u * , w *w , and k w . While the nighttime correlation between w *w , and k w is significant, there is no relationship during the day, when large buoyancy fluxes into surface water (Fig. 6c) suppressed convective mixing. Furthermore, we observed decreases in k w with increasing ΔT (at a given u * or w *w ), such that k w was rarely above 10 cm h −1 when ΔT exceeded $ 4 C.
A closer examination of atmospheric stability conditions provides further evidence for the inferred suppression of gas transfer. As discussed above, atmospheric conditions were frequently in the UVCN zone, where −0.1 < z/L < 0. Such conditions may enhance gas transfer via periodic energetic downdrafts adding turbulence to the air-water interface (Smedman et al. 2007;Sahlée et al. 2008;Andersson et al. 2018). In line with previous findings, we indeed found the greatest departure of k w from k W14 when z/L was in intermediate, UVCN conditions (Fig. 9b). This deviation from prediction occurred in the positive direction (k w enhancement) when ΔT was low, and in the negative direction (k w suppression) when ΔT exceeded 2-3 C. Wind speed was also relatively high during UVCN conditions (Fig. 9a), but could not explain the deviation of k w from prediction, because k w −k W14 was not clearly positive or negative in the UVCN zone under elevated U 10 . While convection and wind-shear appeared to be dominant drivers of CO 2 and latent heat fluxes, the greatest convective enhancement occurred under UVCN conditions. Future studies should more directly address the combined effects of water-side convection and UVCN conditions on turbulent exchanges in shallow water.
Collectively, our findings suggest buoyancy-driven suppression of gas transfer during the day (and convective enhancement at night). This buoyant suppression was broken down in the late afternoon by relatively large latent heat fluxes out of the water (Fig. 6b), which acted to reverse thermal stratification via convective overturn (negative buoyancy flux). These day-time latent heat fluxes imply similarly large evaporative water losses, which should increase the salinity at the sea surface microlayer (Asher et al. 2014). Subsequent density-driven overturn should then act to restore a stable density structure. However, these large latent heat fluxes out of the water also coincided with the greatest under-prediction of k w relative to k W14 (between 15:00 h and 20:00 h, not shown), suggesting continued gas transfer suppression despite latent heat fluxdriven salinization of the sea surface microlayer. These unstable, evaporation-driven, density anomalies can be maintained near the sea surface due to intralayer tension (Wurl et al. Fig. 9. Interaction between UVCN conditions and thermal stratification, including scatter plot of z/L vs. k w −k W14 (b), and U 10 (a). Points are colored by k w −k W14 in (a), and by ΔT in (b). Dotted vertical lines represent UVCN conditions (−0.1 < z/L < 0). Scatter plots for day (right) and night (left) w *w (c, d), u * (e, f), and k w , with the points colored by ΔT. Regression lines are not shown for (e) because the slope was not significantly different from zero (p > 0.05). 2019), perhaps allowing the larger buoyant structure to continue to dampen gas transfer. This is despite the expected effect of density-driven overturn, which should enhance vertical mixing and turbulent exchange.
Together, our results indicate that turbulent fluxes from these shallow waters cannot be explained by a simple stagnant boundary layer model, whereby air-water gas exchange is ratelimited by diffusion across a thin boundary layer. Instead, a more apt conceptual model may be one of eddy diffusion (Lamont and Scott 1970;Wang et al. 2015), whereby gas transfer is driven by the cascade of energy transfer toward molecular dissipation (Wanninkhof et al. 2009;Garbe et al. 2014). The findings of this study, based on direct eddy covariance measurements of CO 2 flux, contribute to our understanding of the physical process of gas transfer, specifically the balance between shear and convective forcing in shallow, coastal waters. Future studies can build on this work by making water-side pCO 2 measurements, facilitating the direct calculations of gas transfer velocity.

Evaporation-precipitation balance in Florida Bay
The clear trend of increasing salinity over the study period, with the transition into the wet season (Fig. 3b), may seem counterintuitive but fits well with the current understanding of Florida Bay's water balance. Net surface water inputs to Florida Bay are small, and have little effect on salinity trends, especially in central and western regions of the bay (Nuttle et al. 2000). Here, the water budget can be considered as a simple evaporation-precipitation balance (Nuttle et al. 2000;Swart and Price 2002;Price et al. 2007), where any small changes in water level are compensated for by slow tidal exchange across sub-basins. Over our study period, cumulative evaporation and precipitation were 518.6 mm and 172 mm, respectively, creating a water deficit of 346.6 mm. Water lost to evaporation is replaced through lateral exchange with adjacent basins and eventually the Gulf of Mexico (Lee et al. 2006(Lee et al. , 2008. Morphologically, Florida Bay is a collection of broad and shallow basins surrounded by mud banks (Fig. 1). These banks are cut by narrow channels, which serve as the only path for cross-basin water exchange (Wang et al. 1994), causing long water residence times, of $ 6-12 months (Lee et al. 2006). Any change in the size of these channels will affect water residence times, causing variable impacts of the evaporation-precipitation balance on salinity trends.
An additional factor affecting the local water balance is clearly the rate of evaporation (LE), which if not balanced by lateral mixing would have caused salinity to rise by 13.5, rather than the observed change of 4.9. While net solar radiation can be a good predictor of evaporation over longer, monthly and seasonal time scales in Florida Bay (Price et al. 2007) and in large and shallow lakes like Taihu (Lee et al. 2014) and Lake Erie (Shao et al. 2015), we found no such relationship between latent heat flux and net solar radiation. Instead, we found that peak latent heat flux (at $ 16:00 h) occurred nearly 4 h after the mid-day peak in net solar radiation (Fig. 6). Such time lags between latent heat flux and net solar radiation have been documented previously in eastern Florida Bay (Price et al. 2007) and in shallow lakes (Czikowsky et al. 2018). Instead, we found that latent heat fluxes were driven by convective (ΔT, w *w ) and shear (U 10 ) forcing in conjunction with vapor pressure deficit. This has been observed previously in a subtropical reservoir (Liu et al. 2009), and related to frontal patterns which can deliver either cold or dry air (enhancing latent heat flux) or warm and moist air (suppressing latent heat flux). We also found that diel trends in latent heat fluxes were opposite of those observed at temperate (Shao et al. 2015;Rey-sánchez et al. 2017) and subtropical (Liu et al. 2009) sites where latent heat flux is typically greater at night. Instead, the high daytime ΔT caused latent heat fluxes to be similar to other tropical sites (Polsenaere et al. 2013). Last, we show that large buoyancy fluxes into the water during the day (11:00 h peak of 4.6 × 10 −7 m 2 s −3 ) effectively suppress daytime latent heat fluxes. At night, buoyancy loss from the water (nighttime average of −1.3 × 10 −7 m 2 s −3 ) drives convective enhancement of latent heat flux (and CO 2 flux).
In summary, observed latent heat flux trends were in stark contrast with other temperate (Shao et al. 2015;Rey-sánchez et al. 2017) and subtropical (Liu et al. 2009) sites where latent heat flux is greatest at night, driven by a cooling of the air relative to the water. While we were able to quantitatively explore the drivers of latent heat flux, and link it with seasonal variations in salinity, more data are required to understand the complete, annual, evaporation-precipitation balance. Additional studies using modeling approaches or rainfall isotopic measurements could explore the ultimate fate of evaporative water losses from Florida Bay, and whether it is of any significance to the freshwater budget of the Florida Everglades.

Regional significance
As discussed in "Temporal trends in air-water CO 2 /latent heat fluxes and associated parameters" section, the unique diel excursion in ΔT at this site (greater ΔT during the day than at night) appears to be related to the marine origin of air reaching central Florida Bay. Due to the presence of the easterly "trade" winds at this latitude, air parcels arriving at this site are near thermal equilibrium with the coastal ocean. This, in conjunction with the shallow water depth (< 2 m) and intense solar heating, causes diel trends in ΔT to be dictated by water temperature rather than air temperature. We document the importance of this diel excursion in ΔT on: (1) increasing pCO 2 (lowering CO 2 solubility), (2) buoyancydriven daytime suppression of turbulent air-water exchange, and (3) nighttime convective enhancement of turbulent fluxes. These factors are likely important at other shallow, subtropical sites, with relatively stable air temperature. In contrast, at shallow coastal sites adjacent to large land masses, continental weather patterns may cause air temperature to be variable, relative to water temperature, causing an opposite diel trend in ΔT (increasing during the night and decreasing during the day). This may be one mechanism by which turbulent fluxes at coastal sites fringing land may be fundamentally different from those at ocean-dominated sites. We expect that air mass origin (either maritime or continental) will have a significant effect on air-water transfer of gas and energy in shallow coastal waters. This topic could be addressed with simple analyses of existing eddy covariance data sets.

Conclusion
We report on the first direct measurements of air-water exchange of CO 2 and heat over a tropical seagrass ecosystem. Against convention, CO 2 flux was largely positive (out of the water), with peak emissions during warm afternoon hours. A more complete explanation of these CO 2 emissions will require water column measurements along with a full annual cycle of direct CO 2 flux determinations by eddy covariance and water pCO 2 measurements. Nevertheless, it appears that $ 23% of the daytime increases in CO 2 flux (likewise, the night-time decrease in CO 2 flux) can be ascribed to the direct, thermodynamic effect of rising water temperature on decreasing CO 2 solubility. The remaining CO 2 emissions may be the combined product of net ecosystem heterotrophy, carbonate precipitation/dissolution, sediment denitrification, and net sulfur biogeochemistry. The relative importance of each process on the carbonate equilibria and ultimate air-water CO 2 exchange is presently unknown. However, it is clear that temperature played a key role in enhancing latent heat fluxes, which generated a water deficit over the short study period of 346.6 mm, causing dramatic increases in salinity. Our correlation analysis indicates that evaporation rates cannot be predicted by a simple relationship with net radiation (as previously assumed), but instead respond to a complex interaction between vapor pressure deficit, air-water temperature differences (ΔT), and convective (w *w ) as well as wind (U 10 ) forcing. As elevated salinity and temperature have been associated with historic seagrass die-offs in Florida Bay and elsewhere, our findings should help ecosystem managers predict when and where future die-off events may occur. In other words, Florida Bay's water budget cannot be understood through solar radiation alone. We describe a suite of factors that should instead be used to forecast changes in temperature and salinity over daily to weekly time scales. Last, we showed that it is possible to estimate the gas transfer velocity (k w ) using heat exchange as a proxy for gas exchange. Using this new approach, we found that k w was less than predicted by commonly used empirical parameterizations based on wind alone or wind and water velocity. This dampening of k w became more prominent with increased daytime ΔT, suggesting that this gas transfer suppression is related to thermal stratification at the sea surface. We point to the apparent contradiction of thermal gas transfer suppression during periods of rapid evaporation, when salinization at the sea surface should instead enhance gas transfer via density-driven vertical mixing. Our results suggest the competing roles of buoyancy fluxes and latent heat fluxes as modulators of gas transfer in shallow, tropical waters. Future studies could build on this work, with water-side measurements of pCO 2 , temperature, and salinity to further address this question.
Together, these findings show that air-water CO 2 exchange in shallow, tropical seagrass ecosystems does not function as other submerged vegetated environments in deeper or more temperate waters. Our results point to temperature as a critical factor governing: (1) evaporative water loss, (2) diel cycles in CO 2 flux (as well as pCO 2 ), and (3) the physical process of airwater gas transfer. By investigating the factors affecting air-sea CO 2 and H 2 O exchange in tropical seagrass meadows, this study will improve our understanding of these important ecosystems in broader regional and global cycles of carbon and water.

Data availability statement
All data sets generated during this project are published on the data sharing repository Figshare (https://doi.org/10.6084/ m9.figshare.13118396.v1). Further requests for data or methods sharing can be directed toward the corresponding author.