In the eddy covariance technique, lateral heat fluxes in the atmosphere, surfaces, and subsurfaces are often ignored under the assumption of a homogeneous surface. Among lateral heat fluxes, the surface and subsurface fluxes, which might affect the surface energy balance closure over a heterogeneous surface, are less studied. Many wetlands are heterogeneous, with a mix of vegetated areas and shallow open water. This study examined the subsurface lateral heat fluxes between a reed bed and the adjacent water of a subtropical wetland in Hong Kong, China. An array of temperature and water-level sensors was installed in the soil of the reed bed and in the adjacent water. An eddy covariance system was also set up on the reed bed. The subsurface lateral heat fluxes were largest close to the interface of the reed bed and water and decreased as the distance from the interface increased, as expected. However, the subsurface lateral heat fluxes could not account for the energy imbalance because their magnitudes were relatively small and fluctuated in phase with the energy budget residuals during the winter months. The uncertainties of the turbulent fluxes and the lateral heat fluxes were estimated to be 10%–20% and 10%–30%, respectively. This study explored another potential reason behind the energy imbalance of the eddy covariance technique. The results enhance the understanding of water and energy exchanges between a terrestrial biotope and the surrounding water, which might further generate insights into the biochemical processes in wetlands.
Observations from single or multiple eddy covariance systems, such as those in the global FLUXNET database (https://fluxnet.ornl.gov/), always lack energy balance closure when the data are analyzed for surface–atmosphere exchanges of energy, momentum, water, and trace gases. There have been various efforts to examine the contribution of heterogeneity to the energy imbalance. Some studies found atmospheric exchange and turbulent motion to be the primary causes of the energy imbalance when the eddy covariance technique was applied over heterogeneous landscapes (Foken et al. 2004; Foken 2008; Stoy et al. 2013). For instance, Kochendorfer and Paw U (2011) identified the role of atmospheric advection in the energy balance using a large array of eddy covariance systems. Higgins et al. (2013) introduced an analytical approach to explore the advection transport mechanism over a land–lake transition and calculated the lateral flux of water vapor in the atmosphere.
The heterogeneous effects are site specific (Kidston et al. 2010; Moncrieff et al. 1996); it is therefore necessary to examine the different land covers and their specific effects. The FLUXNET research sites are of different land covers such as forest, shrubland, cropland, and urban areas. However, there have been very few related studies on wetlands, limiting our understanding about the energy budget over a wetland surface. Wetlands constitute approximately 6% of the Earth’s land surface (OECD 1996), supporting a wide range of wildlife habitats. Understanding each component of the energy budget over the wetland surface can improve our understanding of the biophysical processes, and thus our preservation and restoration of wetland ecosystems (Palanisamy and Chui 2013). For example, evapotranspiration is a key ecohydrologic process in wetlands that links the vegetation, energy budget, and hydrologic system. However, many measurement techniques such as the eddy covariance system (Kochendorfer and Paw U 2011; Stoy et al. 2013) and the use of scintillometer (McGloin et al. 2014) and estimation approaches (Drexler et al. 2004; Abtew and Melesse 2013) are influenced by the heterogeneous effects. Therefore, it is of great significance to investigate the potential reasons behind the energy imbalance over the heterogeneous surfaces of wetlands.
However, even with a comprehensive experiment design such as those of Kochendorfer and Paw U (2011) and Higgins et al. (2013), atmospheric advection induced by different air temperatures and water vapor amounts still cannot fully explain the energy imbalance. Surface heterogeneity can not only change the wind field and moisture distribution in the atmosphere (Froidevaux et al. 2013; Higgins et al. 2013), but also results in different water levels, moisture contents, and temperatures in the subsurface, possibly inducing lateral water and energy exchanges in the subsurface, as shown in Fig. 1. Such subsurface lateral energy exchanges could be a form of energy gain or loss and thus might influence the measurement of energy fluxes and trace gases. This, however, has always been overlooked in previous studies, in which energy components such as net radiation, ground heat flux, and turbulent fluxes were examined without considering the lateral energy exchange with the surrounding environment. Therefore, the lateral heat flux between different landscapes could be a potential reason behind the energy imbalance.
To gain a more integral understanding of the impact of heterogeneous surfaces and to supplement existing studies that mostly examine the atmospheric lateral exchange in nonwetland environments, the conductive and advective lateral heat fluxes in the subsurface of a subtropical wetland were estimated, and their potential relationships with energy imbalance were examined. In this study, it was hypothesized that neglecting lateral heat fluxes in the subsurface could contribute to the energy imbalance over the heterogeneous surface and that the level of energy balance closure could be improved if the lateral heat fluxes were considered.
Field measurements were performed using an eddy covariance system and an array of temperature and water-level sensors. The eddy covariance system was used to obtain the energy budget of a vegetated area, and the sensors were used to capture the temperature and water-level gradients across the interface between the vegetated area and the adjacent water. Footprint analysis was performed to confirm that the source area of the eddy covariance measurements was only the vegetated area, and data of low quality were filtered out. Lateral heat fluxes were estimated, and their spatial variations as well as their temporal fluctuations from December 2015 to September 2016 were investigated. Moreover, the uncertainty of each energy term, including that of the lateral heat fluxes, was quantified. Finally, sensitivity analysis for the thermal and hydraulic conductivity was implemented to examine the potential variations of conductive and advective lateral heat fluxes in different environments.
a. Site description and experiment setup
The field site of this study was within the Mai Po Nature Reserve (MPNR), which is located in northwest Hong Kong, China. Hong Kong has a typical wet, subtropical climate, and MPNR is composed of several different biotopes of vegetated areas and shallow open water, as shown in Fig. 2. The entire experiment was set up in a reed bed within a closed rain-fed pond dominated by Phragmites australis and surrounded by several water channels. The reed bed is an irregular trapezoid (Fig. 2a) with a width of around 250 m and a length of 350 m, and the water channels are about 3–5 m wide and 1 m deep. The site is managed by the World Wildlife Fund for Nature–Hong Kong (WWF-HK), but there was no artificial pumping or release during the monitoring period. An eddy covariance system was set up in the reed bed to obtain basic micrometeorological characteristics and surface energy components, such as the sensible and latent heat fluxes. Considering the potential impact of the internal boundary layer and upscale fetch, the eddy covariance tower was set to be 5.5 m high and 12 m from the water (Fig. 2b). A sonic anemometer (Gill WindMaster) and infrared gas analyzer (model LI-7500) were installed perpendicular to the prevailing wind direction, measuring 10-Hz meteorological variables such as wind velocity. Net radiation was recorded by the net radiometer at 30-min intervals using the model of single-component Kipp and Zonen NR Lite2. The ground heat flux was calculated as the average of the measurements from two soil heat flux plates (model Hukseflux HFP01) that were installed 0.1 m below the ground surface near the eddy covariance tower and recorded using the same time intervals as the other energy components. Other instruments in the system, such as a relative humidity probe, also sampled at 30-min intervals. Moreover, an array of temperature sensors (HOBO U22-001 Water Temp Pro v2 data logger) were spaced horizontally at a depth of 0.2 m in the soil of the reed bed and in the adjacent water channel to record the temperatures at 30-min intervals, marked as Ti in Fig. 2b. Three water-level sensors (In-Situ Aqua TROLL 200 and Level TROLL 300) were installed within a 1-m-deep piezometer at 6 and 12 m from the soil–water interface, as well as in the adjacent water. The piezometers had openings from a depth of 0.5 m to the piezometer bottom and were wrapped with a geotextile fabric to prevent soil from entering the piezometers. Although the piezometers were not screened across the water table, the water levels within the piezometers could be regarded as water table levels because the vertical distances between the water table and the openings were short. The two arrays of sensors together captured the horizontal temperature and water-level gradients, which were used to deduce the subsurface lateral heat flux and water exchange between the reed bed and the adjacent water. The entire monitoring lasted from October 2015 to September 2016; however, monitoring was occasionally interrupted because of equipment failures (i.e., February and March 2016 for the temperature sensors and August and September 2016 for the water-level sensors).
For the energy balance calculation, the red trapezoid that surrounds the reed beds in Fig. 2a was assumed to be the control volume. There were several narrow water channels within the control volume, but their influence on the eddy covariance measurements was ignored because the width of the water channels was negligible compared to that of the entire area. Moreover, wetlands are mostly flat or have small slopes, with water that is shallow and slow moving (Orme 1990; U.S. EPA 2008). The water channels of this study were shallow, with a depth of around 1 m. They were also within an enclosed environment, without artificial water pumping or release, and were subjected to only natural processes such as precipitation, evapotranspiration, and wind. Therefore, the water in the channels was assumed to be stagnant, and the depth of the water channel was assumed to be the same (1 m) throughout the year. All energy components of the control volume, not only the energy over the surface but also the lateral heat flux over each interface of the reed bed and adjacent water, are shown in Fig. 1. There are three water channels embedded within the reed bed (i.e., six interfaces) and one water channel on each of the four sides of the reed bed area (i.e., another four interfaces). Phragmites australis was distributed evenly, and there was almost no other vegetation within the reed bed. Therefore, the lateral heat fluxes to the water channels in different directions were assumed to be equal, and 10 soil–water interface areas were accounted for during the calculation of the subsurface lateral heat fluxes. Furthermore, García et al. (2003) investigated the variations of temperature over the width, length, and depth of reed bed systems that treated wastewater and found that the temperature did not vary much in the direction perpendicular to the water flow. Kadlec (2006) also found that there was no major temperature gradient in the direction perpendicular to the surface flow of a wetland ecosystem. Hence, one transect of the temperature and water-level sensors should be sufficient to identify the major temperature and water-level gradients along the reed bed and open water. Moreover, the residual of the energy budget was in the vertical direction, whereas the lateral heat fluxes were horizontal. To facilitate the comparison between the vertical and horizontal energy fluxes, a ratio between the vertical interface area and the horizontal reed bed surface area (i.e., =/A) was defined, where is the area of the soil–water interface through which the lateral subsurface heat flux was transported (m2), and A is the reed bed surface area, which contributed to the measurement of vertical energy fluxes (m2). Both and A were estimated based on the measurements of width and length using satellite images of Google Earth (Mai Po Nature Reserve, 22.483933°N, 114.037633°E, accessed 25 October 2016).
b. Data preprocessing
Condensation occurs when there is heavy rainfall and can cause unsteady and invalid eddy covariance measurements (Heusinkveld et al. 2008; Burba and Anderson 2010; Mauder and Foken 2011). Hence, data with negative latent heat fluxes, suggesting the occurrence of condensation, were removed. The common assumptions behind the eddy covariance technique include a horizontal homogeneous terrain and stationary airflow with trivial vertical wind velocity (McGloin et al. 2014). Data quality was therefore assessed using the steady-state test and integral turbulence characteristics recommended by Foken et al. (2004), and data potentially violating these assumptions were removed. Göckede et al. (2004) suggested that data of the nonrotated mean vertical wind component exceeding a threshold absolute value of 0.35 m s−1, which corresponds to either strong turbulent effects induced by topography or incorrect sensor orientation, should be excluded. After data quality checking and filtering, a footprint analysis was conducted to confirm that the design of the experiment (e.g., the height of the eddy covariance tower) could eliminate the influence of the lateral heat flux in the atmosphere. The detailed atmospheric characteristics and footprint analysis performed to exclude the atmospheric advection between the reed bed and water are described in appendix A. After all the above filtering, around 7000 data points remained, which were distributed rather evenly over the monitoring period, except for February, when the eddy covariance system did not work because of a power failure.
c. Energy budget of the reed bed
The energy balance budget over the reed bed is shown as Eq. (1), and the methods for calculating each component are illustrated here:
where Rn is the net radiation absorbed by the reed bed; LE is the latent heat flux, which is used by evapotranspiration; H is the sensible heat flux, which is the turbulent exchange with air induced by the temperature difference; G is the heat flux transmitted into the ground and is the sum of the soil heat flux and soil heat storage; S is the energy storage in the air and canopy below the gas analyzer sensor; and and are the surface and subsurface lateral energy exchanges with the surrounding water, respectively. (The units of all the terms are expressed as W m−2.)
Parameters and were measured by the sensors mentioned in section 2a. Based on the eddy covariance data recorded at 10 Hz, LE and H were calculated as the mean covariance of water and sensible heat fluctuations with the vertical wind velocity fluctuations, respectively (Baldocchi et al. 1988):
where is the air density at a given air temperature (kg m−3); is the latent heat of vaporization of water (kJ kg−1); , and denote the fluctuations of vertical wind velocity w (m s−1), air temperature T (°C), and specific humidity q, respectively; is the specific heat capacity of air at constant pressure (J g−1 K−1); and the overbars indicate the average of the product over the monitoring interval.
The storage fluxes of LE and H were estimated from the concentrations of gases and were based on a one-point profile (EddyPro Software 2017). In short-vegetation ecosystems (e.g., crops and grasses), the storage term is expected to be small and is frequently computed on the basis of one single concentration measurement at the height of the eddy covariance flux measurement (Moureaux et al. 2006; Béziat et al. 2009; Aubinet et al. 2012).
The surface lateral heat flux, which is mainly induced by precipitation and runoff, was estimated using (Gosnell et al. 1995):
where is the temperature of the runoff water (°C), is the water temperature in the water channel (°C), is the flow rate of the runoff (m3 s−1), is the heat capacity of water (J g−1 K−1), is the density of water (kg m−3), and A is the surface area of the reed bed (m2). In this study, the runoff temperature was approximated as the wet bulb temperature, and was estimated by the rational method (Kuichling 1889; Thompson 2006), in which the flow rate of the runoff is regarded as proportional to the precipitation:
where the runoff coefficient C is a function of soil type, land use, drainage area slope, etc. (Meshgi et al. 2015) and I is the rainfall intensity (mm h−1).
The subsurface lateral heat flux consisted of two parts: conductive transfer induced by the temperature difference between the reed bed and the adjacent water and the advective heat transfer induced by the water exchange between the reed bed and the adjacent water at various water levels. The conductive heat transfer was calculated by the temperature gradient approach:
where λ is the bulk thermal conductivity (W m−1 K−1), which is dependent primarily upon the soil bulk density and the soil water content (Lu et al. 2007); is the temperature difference between the reed bed and open water (); and and are, respectively, the soil and water temperatures obtained by the temperature sensors (°C). The bulk thermal conductivity was calculated using the weighted mean method based on the distances from the soil–water interface:
where and are, respectively, the horizontal distances (m) of the measurement in the soil and that in the water from the soil–water interface, and and are the thermal conductivities of soil and water, respectively.
The advective heat transfer induced by the different water levels was calculated with
where (°C) is the temperature difference between the soil water temperature and water temperature , and q is the hydraulic loading rate (m day−1), which was estimated by Darcy’s law:
where and are, respectively, the water-level difference of the reed bed and the adjacent water and the distance between them (m), and K is the soil hydraulic conductivity (m day−1), which was estimated using the auger-hole method (Hooghoudt 1934). The rate of water rise within each piezometer after the removal of soil and water was measured and analyzed. Combining Eqs. (9) and (10), the advective energy was calculated as
d. Estimation of vertical temperature profiles
The slope of the reed bed was very small; therefore, the temperature sensors installed at the same depth were assumed to be at the same elevation. The temperature at a depth of 0.2 m was adopted to calculate the lateral heat fluxes. However, temperature varies with depth both in soil and water, which might result in different lateral heat fluxes. To examine the possible impact of such variations with depth, the Van Wijk and De Vries model (Hasfarther et al. 1976; Holmes et al. 2008) was applied to estimate the temperatures at other depths (i.e., 0.1, 0.3, 0.5, and 0.8 m).
with the 24-h moving average temperature (°C) in a temporal resolution of 30 min and a time lag between the two depths:
The damping depth can be calculated from
where W is the thermal diffusivity (m2 s−1), which was assumed to be constant with depth and time (Holmes et al. 2008). The thermal diffusivities were 6.22 × 10−7 m2 s−1 and 1.43 × 10−7 m2 s−1 for soil and water, respectively (Blumm and Lindemann 2006; Oladunjoye and Sanuade 2012). The term is the angular frequency and is a function of the period (s). Because the diurnal variation of the soil and water temperatures corresponds with solar radiation, the period here was assumed to be one day (86 400 s). Because of the different temperatures in summer and winter at the site, two weeks of data, one in summer and one in winter, were analyzed.
e. Uncertainty analysis
Sources of uncertainty in this study were mainly the random uncertainties of the eddy covariance system and the uncertainty in the lateral heat fluxes, which propagated from the uncertainty in the measurements of temperature and water levels. In the surface energy budget [Eq. (1)], the net radiation and soil heat flux were measured by a radiometer and soil heat flux plates, respectively, and their uncertainties were similar at various sites (Twine et al. 2000). The net radiation is probably the most accurate and stable measurement among the major surface energy components, with an accuracy of 5%–7% with reasonable care and knowledge (Halldin and Lindroth 1992; Twine et al. 2000; Foken 2008). Soil heat flux measurements have uncertainties associated with the accuracy of measurements and the spatial variability, giving a random uncertainty of 15%–20% (Twine et al. 2000). Moreover, Twine et al. (2000) suggested that the total random uncertainty in net radiation and soil heat flux is only about half of the uncertainties in turbulent fluxes (e.g., sensible and latent heat fluxes). Therefore, this study focused on the estimation of random uncertainties of the turbulent and lateral heat fluxes.
1) Uncertainty in turbulent flux measurements
Turbulent fluxes averaged over a limited time period have random errors because of the stochastic nature of turbulence (Lenschow et al. 1994; Rannik et al. 2006, 2016). According to Eqs. (2) and (3), the random uncertainty in eddy covariance fluxes is mainly from the calculated covariance between the vertical wind velocity w and the scalar of interest c, which can be evaluated as one standard deviation of the random uncertainty of turbulent flux observed over an averaging period T (T = 30 min in this study; Rannik et al. 2016). Three statistical analysis methods are commonly applied: the “paired tower,” “24-h differencing,” and “model residual” approaches. Given that this study had only one tower and there was no highly tuned empirical model (Aubinet et al. 2012) for turbulent flux estimations, the 24-h differencing approach was adopted to quantify the uncertainty induced by the eddy covariance system (Hollinger and Richardson 2005; Richardson et al. 2006, 2008). Two flux measurements and made at a single tower exactly 24 h apart (to minimize diurnal effects) were regarded as similar and paired measurements. Similar environmental condition criteria were used to include only data in which the difference between and could largely be attributed to random error rather than to different environmental forcing. This study used the criteria suggested by Hollinger and Richardson (2005) (i.e., photosynthetically active photo flux density within 75 μmol m−2 s−1, air temperature within 3°C, wind speed within 1 m s−1, and a vapor pressure deficit within 0.2 kPa) as they yield an acceptable balance between having “similar” environmental conditions and a sufficiently large sample size of measurement pairs so that the uncertainty could be adequately estimated (Richardson et al. 2006, 2008).
2) Uncertainty in lateral heat fluxes
The lateral heat fluxes were calculated from Eqs. (7) and (11), respectively, for conductive and advective heat fluxes between the reed bed and the adjacent water. Random errors could have been induced by the measurements of temperature and water level, as well as by the estimation of thermal and hydraulic conductivities. The Taylor series method (TSM) was applied to calculate the total uncertainty of the lateral heat flux propagated from the uncertainties of all the above variables and coefficients (Taylor and Kuyatt. 1994; Coleman and Steele 2009).
For an experimental result r that can be expressed as a function of J measured variables :
The correlated random error terms have traditionally been assumed to be zero (Coleman and Steele 2009), and so the propagation equation actually used was
where represents the uncertainties in the measured variables . To estimate the total uncertainty propagated from the uncertainty of each variable, is expressed as the relative uncertainty , where is the uncertainty of variable (see appendix B for the detailed transformation and simplification). The uncertainty of the conductive lateral heat fluxes can then be derived as
Similarly, the uncertainty of advective lateral heat fluxes can be obtained by
f. Sensitivity analysis for thermal and hydraulic conductivities
It was difficult to estimate accurately the thermal conductivity of the soil because it not only depends on soil constituents, such as quartz and organic matter contents, but also on the size, shape, and spatial arrangement of soil particles (Hillel 1988; van der Velde et al. 2009). It was also difficult to accurately deduce hydraulic conductivity due to spatial variation (i.e., heterogeneity). To more clearly understand the influence of thermal and hydraulic conductivities, the sensitivity of the lateral heat fluxes to the assumed values of the thermal and hydraulic conductivities was examined. The lateral heat fluxes closest to the soil–water interface were expected to be the largest; hence, the sensitivity analysis was conducted using temperature and water-level measurements closest to the interface. Sensitivity analysis was carried out for the thermal conductivity of the saturated soil at 25°C ranging from 0.6 to 4 W m−1 K−1 and for several different orders of magnitude of hydraulic conductivity: 0.01 m day−1 (e.g., for peat), 1 m day−1 (e.g., for sandy loam), 100 m day−1 (e.g., for well-sorted sand), and 1000 m day−1 (e.g., for well-sorted gravel; Bear 1972; Mesri et al. 1997). The relative variations of conductive and advective lateral heat fluxes as compared to the original values in this study were deduced according to Eqs. (7) and (11), respectively.
a. Surface lateral heat flux
Based on some previous studies (e.g., Chow et al. 1988; Chin 2000), the rational runoff coefficient of the reed bed was estimated to be 0.35. The specific heat capacity of water is about 4.184 J g−1 °C−1 (Perlman 2016), and the surface lateral heat flux was estimated to be based on Eq. (6), with a magnitude of 0.1–1 W m−2, which indicated that the surface lateral heat flux was directly related to the rainfall intensity and the temperature difference between runoff and channel water. However, the rainy period accounted for only around 10% of the monitoring period. Moreover, more than 75% of the data gathered during the rainy period was filtered out based on the filtering criteria explained in section 2b. Given the small amount of data gathered during the rainy period and its seemingly insignificant influence, this study focused on the potential effect of the subsurface lateral heat flux on the energy imbalance; the surface lateral heat flux was excluded from here onward.
b. Horizontal variation of subsurface lateral heat flux
Although the heat capacity and thermal conductivity vary with temperature, the changes should not be significant. For example, the heat capacities of water at 10°, 20°, and 30°C are 4.192, 4.182, and 4.180 J g−1 K−1, respectively (Harr et al. 1984; Marsh 1987; Sengers and Watson 1986), and the corresponding thermal conductivities are 0.579, 0.597, and 0.613 W m−1 K−1, respectively (Farouki 1981). Therefore, the heat capacity and thermal conductivity were assumed to be constant in this study. Applying the respective values of 0.60 and 1.40 W m−1 K−1 for the thermal conductivities of water and of soil with a very high moisture content at 20°C (Johansen 1975; Peters-Lidard et al. 1998), the bulk thermal conductivities at different locations of the reed bed–water channel system were estimated to be 1.0, 1.37, and 1.39 W m−1 K−1 along the direction from the water to the reed bed. Figure 3 displays the temperatures at the different measurement locations and the conductive lateral heat fluxes between those different locations in the reed bed and the adjacent water channel. The temperatures fluctuated in a similar pattern at the different locations, and the soil temperature fluctuations decreased with increasing distance from the interface of the reed bed and water. The water temperature (i.e., T4) fluctuated more dramatically than the soil temperature (i.e., from T1 to T3). For example, the largest temperature difference during the monitoring period was as large as 25°C for the water but about 13°C for the soil. Because of the various temperature gradients, the conductive lateral heat fluxes as well as their fluctuations also decreased with increasing distance from the interface. The conductive lateral heat flux closest to the interface fluctuated from −10 to 10 W m−2. However, at 12 m from the interface, it was close to zero and could be regarded as negligible.
Furthermore, the conductive lateral heat fluxes exhibited obvious linear relationships, as shown in Fig. 4. The differences in , , and (i.e., the conductive lateral fluxes between temperature sensors 1 and 4, 2 and 4, and 3 and 4, respectively) were induced by the temperature difference, the distance, and the bulk thermal conductivity. With increasing distance from the interface, the difference in the bulk thermal conductivities decreased because they were all dominated by the soil thermal conductivity. Therefore, the linearity between and was more apparent, with a better fitness of R2 of 0.992 (p value < 0.001 at the 95% confidence level). The linearities between and and between and were less obvious, with R2 values of 0.872 ( p value < 0.001 at the 95% confidence level) and 0.859 ( p value < 0.001 at the 95% confidence level), respectively.
As shown in Fig. 5a, the water level in the channel was higher than those in the reed bed during most of the monitoring period, and the fluctuation of the water level in the channel was also larger. The water-level fluctuations were about 0.4, 0.25, and 0.15 m in the water channel and 6 and 12 m in the reed bed. Applying the auger-hole method, the hydraulic conductivities (i.e., K6 and K12) were estimated to be 0.008 and 0.003 m day−1 at 6 and 12 m from the interface, respectively. The water exchange induced by the water-level difference was from the water channel to the reed bed during most of the monitoring period (76.6%), which may be due to the rapid transpiration of Phragmites australis (Zhou and Zhou 2009; Headley et al. 2012). The corresponding advective lateral heat flux was small, on the order of 0.01 W m−2, because of the small temperature difference and hydraulic conductivity (Fig. 5b). Similar to the conductive lateral heat flux, the advective heat fluxes 6 and 12 m from the interface also displayed linear relationships, with an R2 value of 0.765 (p value < 0.001 at the 95% confidence level; Fig. 5c).
c. Vertical variation of subsurface lateral heat flux
The temperatures of and the temperature differences between the reed bed and the water at different depths for one week in winter and another week in summer are shown in Fig. 6. In general, the diurnal temperature fluctuations close to the land surface were larger in both soil and water. Moreover, the variation induced by insolation was smaller for soil than for water. For depths larger than 0.3 m, the fluctuation was not obvious, and the temperature difference between the reed bed and the water was insignificant. Therefore, the variations in the shallower depths (i.e., 0.1, 0.2, and 0.3 m) were further examined. In winter, the temperature difference at 0.1 m fluctuated between −4.2° and 5.6°C. The magnitude of the temperature difference was about 2–3 times larger than that at 0.2 and 0.3 m, both of which fluctuated between −1.2° and 2.8°C. In summer, because of stronger radiation, the temperature difference at 0.1 m was about 4–5 times that at 0.2 and 0.3 m. However, even with the largest temperature difference at a depth of 0.1 m in summer, the magnitude of the conductive lateral heat flux was found to be less than 2 W m−2.
d. Comparison between subsurface lateral heat flux and energy budget residual
Before investigating the potential contribution of the subsurface lateral heat flux to the energy imbalance, the energy balance closure (EBC) was first examined, and the results are displayed in Fig. 7. The EBC was evaluated by the slope of the linear regression between incoming and outgoing energy fluxes, excluding the lateral heat fluxes. A value of 1 indicates a perfect energy closure. A value of less than 1 implies that the energy input is greater than the combination of turbulent fluxes, ground heat flux, and storage terms, and vice versa. The EBC values for each month (p values less than 0.001 at the 95% confidence level) are shown in Fig. 7. The straight line in each figure has a slope (i.e., EBC) of 1, covering cases with perfect energy balance closures. The EBC values in August–November were much better than those in other months, especially in April and May. During February, little data were collected because of equipment power failure. The net radiation, being the only energy input, was larger in May–September than in other months, with a magnitude from −100 to 800 W m−2. The sum of turbulent fluxes, ground heat flux, and energy storage displayed similar magnitudes as that of net radiation, and the values fluctuated more in warmer months.
To examine whether the subsurface lateral heat flux could be the potential reason behind the energy imbalance, the conductive and advective lateral heat fluxes were summed up and compared with the residual of the surface energy budget, which is the difference between incoming energy (net radiation) and outgoing energy (ground heat flux, turbulence fluxes, etc.). The comparisons between the lateral heat flux (i.e., the sum of conductive and advective fluxes) and the energy budget residual are shown in Fig. 8. This comparison was performed for the 8 months during which temperature and eddy covariance measurements were both available. The lateral heat flux was rather insignificant as compared to the energy budget residual, with a magnitude of less than 1 W m−2. Even with the largest temperature difference at a depth of 0.1 m in summer, reported in section 3c, the magnitude of the conductive lateral heat flux was less than 2 W m−2; it could therefore be regarded as negligible. The positive value of the lateral heat flux indicated that the energy flow direction was from the reed bed to the adjacent water, which was in the same phase as the energy loss in the reed bed. However, among the 8 months investigated, only 2 months (December 2015 and January 2016) were in the same phase; the other six months were all in the opposite phase. In addition, the residual of the energy budget in August and September 2016 was apparently smaller than that in the other months, which was consistent with the good EBC shown in Fig. 7.
e. Random uncertainty of energy fluxes
1) Random uncertainty of turbulent fluxes
Using the 24-h differencing method, the probability distributions of random uncertainties of sensible and latent heat fluxes were obtained and are displayed in Fig. 9. The environmental condition criteria filtered out around 50% of the observations. A double-exponential (Laplace) distribution provided a good fit and could capture the strong central peak and the outliers. Such a distribution of errors was also found in some previous studies (e.g., Hollinger and Richardson 2005; Rannik et al. 2016). The latent heat flux, being larger in magnitude than the sensible heat flux, also had larger uncertainty (Fig. 9). For the latent heat flux, the mean difference was 1.86 W m−2, with a standard deviation of 61.75 W m−2; for the sensible heat flux, the two statistical parameters were −0.45 and 26.35 W m−2, respectively. This suggests that the eddy flux data are heteroscedastic, with the error increasing with the absolute magnitude of the flux. The uncertainties of the sensible and latent heat fluxes were about 7.2% and 6.5%, respectively, with fluctuations of the sensible heat flux from −162.1 to 205.2 W m−2 and that of latent heat flux from −334.2 to 611.3 W m−2.
2) Random uncertainty of lateral heat fluxes
The uncertainties of lateral heat fluxes as a function of distance x and temperature T were estimated. Assuming reasonable errors of = 0.01 m and = 0.2°C for measurements involving a measurement tape and a temperature sensor (HOBO Pro v2), respectively, then , , and were derived to be 5.2%, 3.5%, and 28.3%, respectively. Then, according to Eq. (18), the uncertainty of the conductive lateral heat fluxes was found to be 29.0%. Similarly, the uncertainty of advective lateral heat fluxes was estimated as 9.3% using Eq. (19).
f. Sensitivity analysis of thermal and hydraulic conductivity
The relative variation of conductive lateral heat fluxes at different soil thermal conductivity values and the relative variation of advective lateral heat fluxes at various hydraulic conductivity values are listed in Table 1. The results suggested that the variation of the conductive lateral heat flux was relatively small at different thermal conductivities of the saturated soil, but the advective lateral heat flux varied dramatically with the different hydraulic conductivities. Various environments with different soil properties can have different orders of magnitude of hydraulic conductivities and water fluxes (Chui et al. 2016). These imply that although the advective lateral heat flux in this study was not significant, it might not be negligible in other cases and might be one reason behind the energy imbalance.
a. Relationship between the lateral heat flux and energy imbalance
Based on the results of the footprint analysis, the source area of the eddy covariance tower was only the reed bed (see appendix A for details). Therefore, there should not be atmospheric advection induced by a heterogeneous surface to account for the energy imbalance. Regarding the subsurface lateral heat flux and taking December 2015 and May 2016 in the first row of Fig. 8 as examples, the lateral heat flux was from the reed bed to the water in December 2015 but in the opposite direction in May 2016. If the small magnitude is neglected, the lateral heat flux from the reed bed to the water might have partly contributed to the “unaccountable” energy loss in December 2015. Then, one could speculate that the lateral heat flux from the water to the reed bed in May 2016 should result in better EBC. However, this speculation was contradicted by the lower EBC in May and the higher EBC in December. Moreover, even for months in which the subsurface lateral heat flux was from the water to the reed bed (e.g., May and August), there were substantial differences among the EBC values (i.e., 0.790 and 0.916, respectively, for May and August; Fig. 7). Furthermore, there was no substantial difference between EBC values considering rainy days and clear days (i.e., 0.81 vs 0.83), which indicates that the surface lateral heat flux did not significantly affect the EBC value. Given the relatively small magnitude of the lateral heat flux compared to the energy budget residual and the flux being an energy gain rather than an “unaccountable” loss in certain months, it was concluded that the subsurface lateral heat flux was not the main cause of the energy imbalance.
b. Factors influencing the subsurface lateral heat flux
1) Influence of the reed bed areal extent
Only with an assumed ratio between the vertical interface area and the horizontal reed bed surface area (i.e., ) could the lateral heat flux be compared to the vertical energy balance residual. However, the larger the interface area or the smaller the reed bed area, the more important the lateral heat flux is when compared with the vertical energy balance residual. Therefore, the lateral heat flux might not always be negligible. For example, Gong et al. (2016) found that in a semiarid ecosystem with low water content, the root uptake and lateral flows from the interspaces of bare soils accounted for over 40% of the annual water loss from the shrub cover. Therefore, ignoring the exchanges between the shrub cover and bare soils in such an environment would result in substantial errors in heat flux estimations. However, when the area with vertical energy exchange is substantially larger than that with lateral energy exchange, the influence of lateral heat flux would become minimal. It is therefore reasonable and common to ignore the influence of lateral heat flux when considering vertical energy budget over a large scale.
2) Influence of other site-specific characteristics
It is also challenging to estimate other influential system characteristics that are site specific. First, both the conductive and advective lateral heat fluxes are induced by temperature differences, and both soil and water temperature fluctuations are driven by air temperature and net radiation variations. The relationship between air and water temperature is relatively easy to obtain. However, that between air and soil is rather site specific, being influenced by various factors such as land surface cover (e.g., vegetation), soil water content, etc. For example, Chudinova et al. (2006) suggested that the soil temperature of frozen ground was influenced by air temperature but also by other conditions, especially snow depth and duration of snow cover. In addition, the conductive and advective lateral heat fluxes are closely related to the thermal and hydraulic conductivities, respectively. The thermal conductivity of nonmetals is approximately constant at high temperatures but it decreases at low temperatures (Debye 1912; Shulman 2004). Lu et al. (2007) proposed a model for predicting soil thermal conductivity for a wide range of water contents at a specific temperature for different soil types and found that the higher the water content, the higher the thermal conductivity, especially for soil with high porosity. The hydraulic conductivity is a proportionality coefficient that expresses the relationship of the rate of water movement to the hydraulic gradient in porous media (Soil Survey Staff 2012). Therefore, it is correlated to soil properties such as pore size, grain size, and soil texture (van Genuchten 1980), and its value ranges from 10−7 m day−1 for unweathered clay to 105 m day−1 for well-sorted gravel. The results of the sensitivity analysis showed that the advective lateral heat fluxes might contribute more to the energy imbalance and should not be neglected in the energy budget of an environment with high hydraulic conductivity.
For the site of this study, given that the soil was saturated and the temperature was relatively high most of the time, the effective thermal conductivity was considered high. On the other hand, the temperature difference between the soil of the reed bed and the adjacent water was not very significant, which might be due to the high water content of the reed bed. Therefore, the conductive lateral heat flux was not very significant. The soil in the reed bed was loam and its hydraulic loading rate was low, which resulted in the insignificant advective lateral heat flux. Moreover, the thermal properties of soil are closely correlated to the soil water content. For this site, the soil of the reed bed was almost saturated throughout the monitoring period, which might also result in the small temperature difference and the heat flux between the reed bed and the adjacent water.
3) Potential influence of gas translocation
Plants in waterlogged soils have internal aeration from aerial parts, and they leak oxygen from their roots to the rhizosphere (Drew and Lynch 1980; Armstrong and Beckett 1987). However, with the formation of aerenchyma and the radial oxygen loss from roots, oxygen permeability was found to decline in the rhizospheres of typical wetland plants (Armstrong et al. 2000). As a defense mechanism, the permeability to phytotoxins and water also declines with decreasing oxygen permeability. Although the rapid transpiration of Phragmites australis might increase the water uptake, the aerenchyma formation and phytotoxin defense might reduce the permeability of oxygen, as well as water permeability. This could be one potential reason behind the small water exchange as well as the insignificant advective lateral heat flux between the reed bed and its adjacent water in this study.
c. Understanding subsurface water and energy exchanges for biochemical processes
The young roots of Phragmites australis and particularly the basal laterals are likely to be the most important oxygen source in the saturated soil (Armstrong and Armstrong 1988). At the soil and water interface, the oxygen could be transported through the interstitial water and affect the biochemical processes in this zone (Gambrell and Patrick 1978). In addition, the energy fluxes between the water and the biotope could drive oxidation–reduction reactions (Kadlec and Wallace 2009), which could shift the aeration status and further influence the microbial and chemical processes, as well as the biological availability of major nutrients in the soil (Patrick et al. 1985; Gambrell et al. 1976; Dušek et al. 2008). Moreover, for constructed wetlands, pollutant removal is sensitive to internal flow patterns. Pollutants often dissolve in water and are transported to biofilms for biochemical transformation (Kadlec and Wallace 2009). Therefore, understanding the subsurface lateral water and energy exchanges in the wetland system might benefit the study of these related biochemical processes.
d. Study limitations
There were several limitations in this study. First, even though the reeds were distributed rather evenly over the reed bed and the temperature gradient should be small in the direction perpendicular to the flow of water and energy (Garcia et al. 2003; Kadlec 2006), more transects of temperature and water-level measurements at different locations might give more representative results. Moreover, the water within the channels was assumed to be stagnant and to have a constant depth of 1 m throughout the monitoring period. However, the water level fluctuated by around 0.2 m, mainly due to precipitation and evaporation. Such water-level fluctuations would not change the order of magnitude of the lateral heat fluxes and would not significantly affect the conclusions of this study, but they would lead to some uncertainty. Last but not least, the reed bed was assumed to be homogenous, and the hydraulic conductivity was measured in each of the piezometers once. The possible spatial variability of the hydraulic conductivity could therefore have led to some errors during the estimation of lateral advective heat fluxes.
This study quantified the subsurface conductive and advective lateral heat fluxes between a reed bed and the adjacent water in the heterogeneous environment of a subtropical wetland and examined their potential contributions to the surface energy imbalance. The surface energy balance closure varied in different months, ranging from 0.790 to 0.916. The subsurface lateral heat fluxes were largest close to the interface of the reed bed and the water and decreased as the distance from the interface increased, as expected. However, the magnitudes of both the conductive and advective lateral heat fluxes were not substantial when compared to the other energy terms (e.g., sensible and turbulent heat fluxes), particularly the advective lateral heat flux, which was on the order of just 0.01 W m−2. Moreover, the subsurface lateral heat flux was in phase with the energy budget residual only in some months but out of phase for most of the time. Therefore, the subsurface lateral heat fluxes could not account for the energy imbalance because of their small magnitudes and phase differences. The uncertainties of turbulent fluxes (i.e., sensible and latent heat fluxes) and lateral heat fluxes were respectively estimated to be about 10%–20% and 10%–30%. However, given that both conductive and advective lateral heat fluxes were not substantial, the uncertainties would not affect the above conclusions. Furthermore, the sensitivity analysis of thermal and hydraulic conductivity indicated that the conductive lateral heat flux might not vary much among various landscapes. However, the advective lateral heat flux could vary dramatically because of the wide range of hydraulic conductivities for different soil types. Although the subsurface lateral heat fluxes were concluded to not be the main reason behind the energy closure problem of the study site, this study enhanced our understanding of the water and energy exchanges between a terrestrial biotope and its surrounding water, which might further generate some insights into the biochemical processes in the heterogeneous environment of wetlands.
This work was funded by the General Research Fund of the Research Grants Council, Hong Kong (project code 17202114). The authors thank the Agriculture, Fisheries and Conservation Department of Hong Kong and WWF-Hong Kong for facilitating the field work at Mai Po Nature Reserve. They would also like to thank their colleague, Chan Tsz On, for his assistance during the experimental setup and field data collection.
Footprint Analysis and Seasonal Variation of Wind
a. Footprint analysis
The upwind-contributing or source areas (i.e., footprints) of the eddy covariance system vary with measurement height, atmospheric stability, surface roughness length, etc. Two footprint models were applied to determine and confirm that the source area was within the reed bed, so that there was minimal atmospheric advection induced by heterogeneity.
The first model is the “simple footprint parameterization” model described by Kljun et al. (2004), which assumes a homogeneous surface. It is a three-dimensional backward Lagrangian footprint model that can determine the source area for a wide range of atmospheric stabilities (Kljun et al. 2004). This model estimates the distance of the location that contributes the most to the measurements and the percentages of the measurements from locations at various distances from the eddy covariance system. The other model is SCADIS, which is the footprint model developed by Sogachev and Lloyd (2004) and Sogachev and Sedletski (2006). It describes the scalar distribution and can be utilized over a heterogeneous surface. SCADIS is essentially a canopy–atmospheric boundary layer and a scalar transport one-and-a-half-order closure model (Sogachev and Lloyd 2004), which is based on the Reynolds equations for turbulent flow transformed from the Navier–Stokes equations (Monin and Yaglom 1971) and considers the temporal variation and horizontally averaged momentum. The first model estimates the distance of the location that contributes the most to the measurements and the percentages of the measurements from locations at various distances from the eddy covariance system. The second model gives the crosswind integrated footprint function , which presents the relative contribution of each upwind surface element to the measured scalar flux, as well as the cumulative footprint function along the flow direction.
Although several water channels were embedded within the contributing area, their widths were only 3–5 m, and their contributions were neglected given that the entire area was 250 m × 350 m. The distribution of locations that contributed the most to the eddy covariance measurements (Fig. A1a) and the cumulative distribution function of distances from the eddy covariance system within which 90% of the collected data originated from (Fig. A1b) were obtained from Kljun’s model. The spatial variation of contributions (Fig. A1c) was obtained from the SCADIS footprint calculation. The areas that contributed the most to the eddy covariance measurements were mostly within 40 m. Also, the source area that contributed 90% of the measurements was mostly within 250 m from the eddy covariance system. Moreover, the results from SCADIS indicated that the contribution from locations beyond 250 m was very minimal. Therefore, the results of both footprint models suggested that the effective source area contributing the turbulent fluxes was within 250 m.
b. Seasonal variation of wind
The wind roses derived using the wind velocity from the eddy covariance system are shown in Fig. A2 for different seasons. The predominant wind directions in winter, spring, and autumn were similar, that is, northeast. However, because of the tropical monsoon climate, the primary wind direction in summer was southeast. The average wind speed at the study site in winter and spring was relatively larger than that in summer and autumn. Combining the results of the wind roses and the site surface characteristics, the predominant wind was considered to be from the reed bed for all seasons except for a small percentage coming from the open water in the west in summer and autumn (Fig. A2).
Combining the analysis from both parts, the effect of atmospheric advection between the reed bed and the adjacent water was minimized. The data during which the wind was coming from the water or out of the reed bed were also filtered out.
Taylor Series Method for Uncertainty Analysis
For an experimental result r that could be expressed as a function of J measured variables :
An analogous TSM equation for the random standard uncertainty of the calculated result is (Coleman and Steele 2009)
where represents the uncertainties in the measured variables . There is a correlated random error term containing a covariance factor for each pair of measured variables whose “random” variations are not independent from each other. It is assumed that the relationship given by Eq. (B2) is continuous and has continuous derivatives in the domain of interest. The correlated random error terms have traditionally been assumed to be zero (Coleman and Steele 2009), so that the propagation equation actually used was
To obtain the relationship between the total uncertainty and the partial uncertainty contributed by each independent variable, each term of Eq. (B3) was first divided by . Each term on the right-hand side of Eq. (B3) was also multiplied by :
where is the relative uncertainty of the result and represents the relative uncertainties for each variable. The uncertainty magnification factors (UMFs) are defined as
The UMF for a given indicates the influence of the uncertainty in that variable on the uncertainty in the final result. The influence of the uncertainty in the variable is magnified as it propagates through to the final result, when the UMF value is greater than 1 and vice versa.
In this study, the conductive and advective lateral heat fluxes were respectively calculated using Eqs. (7) and (11). Taking Eq. (7) as an example and applying Eq. (B4) to estimate the uncertainty of conductive lateral heat fluxes,
The UMFs are