Submarine Groundwater Discharge in a Coastal Bay: Evidence from Radon Investigations

: Jiaozhou Bay, an urbanized coastal bay located in the southern part of Shandong Peninsula, China, has been deeply a ﬀ ected by anthropogenic activities. Here, the naturally occurring 222 Rn isotope was used as a tracer to assess the submarine groundwater discharge (SGD) in this bay. The time series of 222 Rn concentrations in nearshore seawater were monitored continuously over several tidal cycles at two ﬁxed sites (Tuandao (TD) and Hongdao (HD)) during the dry season in spring and the wet season in autumn of 2016. 222 Rn concentrations in seawater were negatively related to the water depth, indicating the inﬂuence of tidal pumping. A 222 Rn mass balance model revealed that the mean SGD rates were 21.9 cm / d at TD and 17.8 cm / d at HD in the dry season, and were 19.5 cm / d at TD and 26.9 cm / d at HD in the wet season. These rates were about 8–14 times the discharge rates of the local rivers. Enhanced groundwater inputs occurred at HD in the wet season, likely due to the large tidal amplitudes and the rapid response to local precipitation. Large inputs of SGD may have important inﬂuences on nutrients levels and structure, as well as the water eutrophication occurring in coastal waters.


Introduction
Submarine groundwater discharge (SGD) is defined as the flow of all continental margin waters from the seabed to the coastal ocean, which includes submarine fresh groundwater discharge (SFGD) and recirculated saline groundwater discharge (RSGD) [1,2]. As an important process of land-ocean interactions, SGD has received increasing attention in recent decades in the field of hydrogeology. Accumulating evidence has shown that SGD flux can be comparable to or even higher than riverine flux [2][3][4][5]. On a global scale, SGD flux is 3 to 4 times greater than the freshwater fluxes into the oceans by rivers [6]. Moreover, SGD contains high levels of nutrients [7][8][9], carbon [10,11], and metals [12,13]. As a result, a small SGD flux may carry a large quantity of dissolved materials from land to ocean [14,15]. These large fluxes of dissolved materials can lead to water eutrophication and acidification in coastal areas [16][17][18], which seriously influence marine ecosystems [19,20]. Therefore, it is of great significance to correctly understand and evaluate the SGD flux for coastal environment protection.
A heavily urbanized coastal bay (Jiaozhou Bay, China) has been greatly influenced by anthropogenic activities such as rapid population growth and industrial expansion. From 1966 to 2008, an area of 180 km 2 of land has been reclaimed from the coastal regions by backfilling along the northern and western coastline of Jiaozhou Bay [35]. The nutrient concentrations (dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP)) in Jiaozhou Bay were increasing from the 1960s to the 1990s [36,37]. The bay waters underwent seasonal eutrophication and red tides occurred frequently because of large discharges of pollutants into the bay [38,39]. However, previous studies mainly have focused on regional SGD and associated nutrient fluxes in Jiaozhou Bay [40,41]. Little information is available on the temporal variations of SGD in this region. Therefore, it is of great significance to study the variability of SGD with different seasons and sites in Jiaozhou Bay. In this study, the variations of SGD were assessed for Jiaozhou Bay over several tidal cycles at two fixed sites (Tuandao and Hongdao) in spring and autumn, based on continuous 222 Rn observations. The difference in SGD between different seasons and sites was analyzed, and the relationship between the SGD and multiple influence factors, including water depth, season, location, and local river discharge, were explored.

Study Site
Jiaozhou Bay, a typical semiclosed bay (35°58′ N~36°18′ N, 120°04′ E~120°23′ E), is located in the southern part of Shandong Peninsula, China (Figure 1), linking to the Yellow Sea through a narrow channel 3.1 km wide. The average water depth is 7 m, and the bay mouth has a maximum water column depth of 60~70 m [42]. The average tidal range is 2.71 m and the largest tidal range is 6.87 m [43]. Tidal prism in 2008 decreased by 28% compared to the value in 1935 and the water residence time increased from 36 to 41 days due to land reclamation [44]. The water area had been reduced to 300 km 2 by 2015, due to the influence of human activities, and tidal flats occupy approximately one-third of the total water area [45][46][47]. There are several different types of terrains surrounding the bay: Jimo Basin in the north, an alluvial plain in the northwest, Laoshan Mountain in the east, and Xiaozhu Mountain in the southwest [48]. The west and north of the bay have Quaternary unconsolidated sediments, which There are several different types of terrains surrounding the bay: Jimo Basin in the north, an alluvial plain in the northwest, Laoshan Mountain in the east, and Xiaozhu Mountain in the southwest [48]. The west and north of the bay have Quaternary unconsolidated sediments, which can easily store groundwater and be recharged by rainfall. Figure 2 shows the daily precipitation and evaporation in 2016 with an annual mean rainfall of~484 mm. The mean precipitation in spring 2016 and autumn 2016 (during the sampling period for this study) was 0.95 mm/d and 1.33 mm/d, and the mean evaporation was 3.37 mm/d and 2.62 mm/d, respectively. Many seasonal rivers were flowing into Jiaozhou Bay, including Yanghe River, Moshui River, Licun River, Dagu River, and Baisha River [42,49]. The Dagu River contributes at least 80% of the freshwater from the local rivers each year, with an average discharge rate of 1.98 × 10 6 m 3 /d [50,51].
Water 2020, 12, x FOR PEER REVIEW  3 of 16 can easily store groundwater and be recharged by rainfall. Figure 2 shows the daily precipitation and evaporation in 2016 with an annual mean rainfall of ~484 mm. The mean precipitation in spring 2016 and autumn 2016 (during the sampling period for this study) was 0.95 mm/d and 1.33 mm/d, and the mean evaporation was 3.37 mm/d and 2.62 mm/d, respectively. Many seasonal rivers were flowing into Jiaozhou Bay, including Yanghe River, Moshui River, Licun River, Dagu River, and Baisha River [42,49]. The Dagu River contributes at least 80% of the freshwater from the local rivers each year, with an average discharge rate of 1.98 × 10 6 m 3 /d [50,51].

Time-Series Deployments
Continuous 222 Rn monitoring was conducted at two fixed sites during the year's contrasting hydrological conditions (dry season in spring and wet season in fall). The first test site in Tuandao (TD) was located in the south, near the mouth of the bay. The second site in Hongdao (HD) was located in the north and surrounded by a more enclosed section of the bay, compared to TD ( Figure  1). The two stations were separated by approximately 14.4 km. The time-series measurements of 222 Rn concentrations were performed using a radon-in-air monitor and RAD AQUA [52]. Radon counts were integrated over 1 h cycles, giving measurement uncertainties of 10 ~ 30% for each data point. In the dry season, continuous measurements of 222 Rn concentrations in the nearshore seawater were conducted from 9 to 11 April 2016 (49 h) at TD, and from 11 to 14 April 2016 (68 h) at HD. In the wet season, 222 Rn concentrations were monitored from 10 to 12 November 2016 (50 h) at TD, and from 13 to 15 November 2016 (47 h) at HD.
During each time-series deployment, seawater was continuously pumped from ~0.5 m above the sea bottom and filtered through a cartridge filter to screen out any particulates before entering the AQUA system. The air-water exchanger was connected to the RAD7 [52] detector by a closed air loop. To ensure the accuracy of measurement, no bubble was introduced into the system. The 222 Rn concentration in seawater was determined by multiplying the 222 Rn concentration in the water column and the partition coefficient (α ). Schubert et al. introduced an easily applied equation to derive the partition coefficient [53], which depends on both water temperature and salinity, as shown in Equations (1-3): ) 100 100 ( ) 100 ln( ) 100 ( ln where A( 222 Rn) and A( 222 Rnm) are the 222 Rn concentrations (Bq/m 3 ) before and after calibration, respectively; α is the partition coefficient; β is the Bunsen coefficient; T is the temperature (°C); S is salinity (ppt); and a1-a3 and b1-b3 are six variable parameters.

Time-Series Deployments
Continuous 222 Rn monitoring was conducted at two fixed sites during the year's contrasting hydrological conditions (dry season in spring and wet season in fall). The first test site in Tuandao (TD) was located in the south, near the mouth of the bay. The second site in Hongdao (HD) was located in the north and surrounded by a more enclosed section of the bay, compared to TD (Figure 1). The two stations were separated by approximately 14.4 km. The time-series measurements of 222 Rn concentrations were performed using a radon-in-air monitor and RAD AQUA [52]. Radon counts were integrated over 1 h cycles, giving measurement uncertainties of 10~30% for each data point. In the dry season, continuous measurements of 222 Rn concentrations in the nearshore seawater were conducted from 9 to 11 April 2016 (49 h) at TD, and from 11 to 14 April 2016 (68 h) at HD. In the wet season, 222 Rn concentrations were monitored from 10 to 12 November 2016 (50 h) at TD, and from 13 to 15 November 2016 (47 h) at HD.
During each time-series deployment, seawater was continuously pumped from~0.5 m above the sea bottom and filtered through a cartridge filter to screen out any particulates before entering the AQUA system. The air-water exchanger was connected to the RAD7 [52] detector by a closed air loop. To ensure the accuracy of measurement, no bubble was introduced into the system. The 222 Rn concentration in seawater was determined by multiplying the 222 Rn concentration in the water column and the partition coefficient (α). Schubert et al. introduced an easily applied equation to derive the partition coefficient [53], which depends on both water temperature and salinity, as shown in where A( 222 Rn) and A( 222 Rn m ) are the 222 Rn concentrations (Bq/m 3 ) before and after calibration, respectively; α is the partition coefficient; β is the Bunsen coefficient; T is the temperature ( • C); S is salinity (ppt); and a 1 -a 3 and b 1 -b 3 are six variable parameters.
In parallel with continuous monitoring of 222 Rn, the electrical conductivity, water temperature, and water levels were measured automatically in 30 min intervals by a CTD-Diver. The barometric pressure and temperature were recorded by a Baro-Diver placed in the air near the continuous 222 Rn monitoring sites. Continuous 222 Rn concentration in the atmosphere was also measured by an additional RAD7 detector. Both RAD7 detectors were programmed to integrate the 222 Rn counts every hour, whether in the air or the water. Wind speed was monitored by a portable hand-held anemometer over a 30 min period.

226 Ra Enrichment and Analysis
To estimate the 222 Rn supported by 226 Ra in seawater, 60 L of seawater was collected from the 222 Rn monitoring sites to extract 226 Ra. After filtering, the seawater was slowly passed through a polyamide fiber, containing about 25 g of Mn fibers at a flow rate of less than 1 L/min. The fibers were thoroughly washed with distilled water to remove all particles and salts, and then the fibers were taken to the laboratory for measurement. In the laboratory, the fiber samples were sealed for 22 days (7 times the half-life of 222 Rn). After 226 Ra reached its decay equilibrium with 222 Rn and its daughter, the 226 Ra concentrations were analyzed with a RAD7 detector, following the method of Kim et al. [54]. To improve the accuracy, the measurement time of each sample was increased to 12 h. The measurement uncertainty of 226 Ra was 10-20%.

Determination of 222 Rn in Pore Water
Six sediment samples from Jiaozhou Bay were collected to evaluate the 222 Rn concentration in pore water based on sediment equilibration experiments [55] (Figure 1). In the laboratory, 100 g sediment samples were mixed with 500 mL of seawater in radon-free Erlenmeyer flasks, and then sealed and oscillated continuously for about 30 days. Given the half-life of 222 Rn, the 226 Ra and 222 Rn in the pore water of sediments and the overlying water reach decay equilibrium after~4 weeks. The overlying water was then transferred to a 250 mL sampling bottle and measured using a RAD7 detector with a RAD-H 2 O accessory [52].

Radon Mass Balance Model
Burnett et al. developed a radon mass balance model based on continuous 222 Rn observations [31]. In seawater systems, the 222 Rn sources included SGD input, sediment diffusion, 226 Ra decay, and tidal input. Sinks of 222 Rn included tidal output, atmospheric evasion, decay loss, and mixing loss. The radon mass balance is shown in Equation (4): where dI/dt is the difference in flux of 222  Note that the radon mass balance model alone did not allow a separation of these two SGD components of SFGD and RSGD [56]. Therefore, a simple two end-member mixing model on salinity was used to estimate the fresh groundwater fraction for the two seasons [57,58]: where S GW is the groundwater salinity in the mixing zone; S fw and S bay are the average salinities of fresh groundwater and bay water, respectively; and f fw is the fraction of SFGD.

222 Rn Inventory and 226 Ra Supply
In the dry season, the 222 Rn concentrations during the monitoring period ranged from 17.0 to 136.7 Bq/m 3 (average ± standard deviation: 57.8 ± 30.0 Bq/m 3 ) at TD and from 9.7 to 175.5 Bq/m 3 (77.4 ± 30.7 Bq/m 3 ) at HD. In the wet season, the 222 Rn concentration ranged from 13.1 to 68.4 Bq/m 3 (31.0 ± 11.4 Bq/m 3 ) at TD, and from 41.9 to 293.8 Bq/m 3 (111.9 ± 48.4 Bq/m 3 ) at HD (Figure 3). The average tidal range was larger in the dry season (4.28 m) than that in the wet season (2.86 m) at TD. On the contrary, there was a lower tidal range in the dry season (3.31 m) than in the wet season (4.18 m) at HD. As shown in Figure 3, the concentration of 222 Rn changed periodically with the tides and had a negative correlation with water depth. This indicated that low and high concentration of 222 Rn in seawater was occurred during flood tide and ebb tide, respectively. The general fluctuations of 222 Rn concentration during the monitoring period were known to respond to tidal pumping and the hydraulic gradient. Most likely, the hydraulic gradient between seawater and groundwater increased at low tides, which caused large amounts of seepage and higher 222 Rn concentrations. Conversely, the hydraulic gradient decreased at high tides. The 222 Rn concentration in the nearshore water was diluted by mixing with offshore water at high tides, which contributed to less seepage and lower 222 Rn concentrations. Interestingly, there was a lag time between low tide and the peak of 222 Rn concentration in this study, and the lag times were different at different low tides ( Figure 3). This phenomenon was consistent with previous studies [9,59]. Wu et al. showed that the phase delay may be associated with the distance from the shoreline, water depth, and topographic conditions [59].
Water 2020, 12, x FOR PEER REVIEW 5 of 16

222 Rn Inventory and 226 Ra Supply
In the dry season, the 222 Rn concentrations during the monitoring period ranged from 17.0 to 136.7 Bq/m 3 (average ± standard deviation: 57.8 ± 30.0 Bq/m 3 ) at TD and from 9.7 to 175.5 Bq/m 3 (77.4 ± 30.7 Bq/m 3 ) at HD. In the wet season, the 222 Rn concentration ranged from 13.1 to 68.4 Bq/m 3 (31.0 ± 11.4 Bq/m 3 ) at TD, and from 41.9 to 293.8 Bq/m 3 (111.9 ± 48.4 Bq/m 3 ) at HD (Figure 3). The average tidal range was larger in the dry season (4.28 m) than that in the wet season (2.86 m) at TD. On the contrary, there was a lower tidal range in the dry season (3.31 m) than in the wet season (4.18 m) at HD. As shown in Figure 3, the concentration of 222 Rn changed periodically with the tides and had a negative correlation with water depth. This indicated that low and high concentration of 222 Rn in seawater was occurred during flood tide and ebb tide, respectively. The general fluctuations of 222 Rn concentration during the monitoring period were known to respond to tidal pumping and the hydraulic gradient. Most likely, the hydraulic gradient between seawater and groundwater increased at low tides, which caused large amounts of seepage and higher 222 Rn concentrations. Conversely, the hydraulic gradient decreased at high tides. The 222 Rn concentration in the nearshore water was diluted by mixing with offshore water at high tides, which contributed to less seepage and lower 222 Rn concentrations. Interestingly, there was a lag time between low tide and the peak of 222 Rn concentration in this study, and the lag times were different at different low tides ( Figure 3). This phenomenon was consistent with previous studies [9,59]. Wu et al. showed that the phase delay may be associated with the distance from the shoreline, water depth, and topographic conditions [59]. As shown in Figures 3a,d, the salinity has a significantly positive correlation with water depth, indicating that the salinity peaked at high tides, and the variation of salinity tended to increase with time, which may be due to an increased inflow of high salinity seawater from the open ocean, accompanied by a gradual increase in tidal range. The salinity fluctuation (29.5-36.3 ppt) at HD in the dry season (Figure 3b) was the largest (Figure 3). Figure 4 shows the temporal variation of 226 Ra in seawater at both sites in two seasons. In the dry season, the 226   As shown in Figure 3a,d, the salinity has a significantly positive correlation with water depth, indicating that the salinity peaked at high tides, and the variation of salinity tended to increase with time, which may be due to an increased inflow of high salinity seawater from the open ocean, accompanied by a gradual increase in tidal range. The salinity fluctuation (29.5-36.3 ppt) at HD in the dry season (Figure 3b) was the largest (Figure 3). Figure 4 shows the temporal variation of 226 Ra in seawater at both sites in two seasons. In the dry season, the 226   To evaluate the 222 Rn flux attributed to SGD, it is necessary to calculate the variations in storage flux of 222 Rn (dI/dt) in seawater between two adjacent time intervals. The excess 222 Rn inventory was defined as the product of the excess 222 Rn concentration (subtracting 226 Ra from total 222 Rn) and water depth [31]: where I is the 222

Tidal Effects
Radon inventories in seawater can be strongly controlled by tides. 222 Rn from the open sea can flow into the seawater system with the incoming water during the flood tide and would flow out of the seawater system with the outgoing water during the ebb tide. Following Zhang et al. [9], the flux of 222 Rn (Ftide) subjected to tide transportation can be expressed using Equation (7) where  To evaluate the 222 Rn flux attributed to SGD, it is necessary to calculate the variations in storage flux of 222 Rn (dI/dt) in seawater between two adjacent time intervals. The excess 222 Rn inventory was defined as the product of the excess 222 Rn concentration (subtracting 226 Ra from total 222 Rn) and water depth [31]: where I is the 222

Tidal Effects
Radon inventories in seawater can be strongly controlled by tides. 222 Rn from the open sea can flow into the seawater system with the incoming water during the flood tide and would flow out of the seawater system with the outgoing water during the ebb tide. Following Zhang et al. [9], the flux of 222 Rn (F tide ) subjected to tide transportation can be expressed using Equation (7): 222 Rn o f f when f lood tide h t+∆t −h t ∆t × 222 Rn t,t+∆t when ebb tide (7) where h t+∆t and h t are the water depth (m) at time t and t + ∆t, respectively; 222

Atmospheric Loss
Radon is a slightly soluble gas in water. When the 222 Rn concentration in water is greater than that in air, 222 Rn will diffuse into the atmosphere through the water-air interface. The atmospheric loss of 222 Rn was mainly related to its migration coefficient and the concentration gradient of the water-gas interface [61]: where k is the gas transfer velocity (m/h); C w and C air are the 222 Rn concentrations in seawater and atmosphere, respectively (Bq/m 3 ). α is the partition coefficient related to temperature and salinity (see Equation (2)). When C w > αC air , 222 Rn diffuses from the water column to the air. Lambert et al. give the expression of the gas transfer velocity k of 222 Rn [62]: where µ is the wind speed (m/s) and S c is Schmidt constant for radon at a given water temperature [63].

Atmospheric Loss
Radon is a slightly soluble gas in water. When the 222 Rn concentration in water is greater than that in air, 222 Rn will diffuse into the atmosphere through the water-air interface. The atmospheric loss of 222 Rn was mainly related to its migration coefficient and the concentration gradient of the water-gas interface [61]: where k is the gas transfer velocity (m/h); Cw and Cair are the 222 Rn concentrations in seawater and atmosphere, respectively (Bq/m 3 ). α is the partition coefficient related to temperature and salinity (see Equation (2)). When 222 Rn diffuses from the water column to the air. Lambert et al. give the expression of the gas transfer velocity k of 222 Rn [62]: where μ is the wind speed (m/s) and Sc is Schmidt constant for radon at a given water temperature [63]. Figure 5 shows

Diffusive Flux From Bottom Sediments
Radon exchange between sediments and the overlying seawater was affected by the concentration gradient, as well as biological and physical disturbances. Sediment source of radon is important due to radon's short half-life. The 222 Rn flux through sediment diffusion was obtained by performing a sediment cultivation experiment of radon in this study [64]: where Fsed is the diffusion flux of 222 Rn from sediments (Bq/m 2 h); λ is the decay constant of 222 Rn (0.181d −1 ); Dm is the molecular diffusion coefficient related to temperature; and Ceq is the 222 Rn concentration (Bq/m 3 ) in the pore water of sediments. C0 is the 222 Rn concentration (Bq/m 3 ) in the overlying seawater. After correcting the wet density and porosity of the sediment, Ceq and Dm can be obtained by the following equations: where Csed is the 222 Rn concentration in wet sediments (Bq/kg); V is the volume of the overlying seawater in the sediment cultivation experiment; wet ρ is the wet density of the sediments (kg/m 3 ); and Msed is the mass of the sediments (kg). According to the laboratory particle analysis, the porosity of the sediments in Jiaozhou Bay was 0.3 and the wet density of the sediments was 1205 kg/m 3 . The average 222 Rn concentration in the pore water obtained from the sediment cultivation experiment was 4745 Bq/m 3 at TD (S1-S3 in Figure 1) and 4788 Bq/m 3 at HD (S4-S6 in Figure 1). In the dry season, the average value of the molecular diffusion coefficient (Dm) was 8.75 × 10 −6 cm 2 /h at TD and 9.86 × 10 −6 cm 2 /h at HD. Based on Equation (11), the variations in 222 Rn flux across the sediment-water interface ranged from −2.88 × 10 −2 to 3.46 × 10 −2 Bq/m 2 h at TD and from −6.78 × 10 −2 to 2.88 × 10 −2 Bq/m 2 h at HD. In the wet season, the average value of the molecular diffusion coefficient (Dm) was 1.03 × 10 −5 cm 2 /h at TD and 9.70 × 10 −6 cm 2 /h at HD. The variations in 222 Rn flux across the sediment-water interface ranged from −1.97 × 10 −2 to 2.39 × 10 −2 Bq/m 2 h at TD and from −5.50 × 10 −2 to 3.75 × 10 −2 Bq/m 2 h at HD. Compared to the flux of 222 Rn from tides, the diffusion flux

Diffusive Flux From Bottom Sediments
Radon exchange between sediments and the overlying seawater was affected by the concentration gradient, as well as biological and physical disturbances. Sediment source of radon is important due to radon's short half-life. The 222 Rn flux through sediment diffusion was obtained by performing a sediment cultivation experiment of radon in this study [64]: where F sed is the diffusion flux of 222 Rn from sediments (Bq/m 2 h); λ is the decay constant of 222 Rn (0.181d −1 ); D m is the molecular diffusion coefficient related to temperature; and C eq is the 222 Rn concentration (Bq/m 3 ) in the pore water of sediments. C 0 is the 222 Rn concentration (Bq/m 3 ) in the overlying seawater. After correcting the wet density and porosity of the sediment, C eq and D m can be obtained by the following equations: where C sed is the 222 Rn concentration in wet sediments (Bq/kg); V is the volume of the overlying seawater in the sediment cultivation experiment; ρ wet is the wet density of the sediments (kg/m 3 ); and M sed is the mass of the sediments (kg). According to the laboratory particle analysis, the porosity of the sediments in Jiaozhou Bay was 0.3 and the wet density of the sediments was 1205 kg/m 3 . The average 222 Rn concentration in the pore water obtained from the sediment cultivation experiment was 4745 Bq/m 3 at TD (S1-S3 in Figure 1) and 4788 Bq/m 3 at HD (S4-S6 in Figure 1). In the dry season, the average value of the molecular diffusion coefficient (D m ) was 8.75 × 10 −6 cm 2 /h at TD and 9.86 × 10 −6 cm 2 /h at HD. Based on Equation (

Mixing Loss and SGD Flux
After dI/dt was corrected for atmospheric loss, tidal effects, and diffusion from sediments, the net 222 Rn flux could be estimated, i.e., the algebraic sum of F SGD and F mix . The radon decay loss in Equation (4) was not considered because the decay loss of 222 Rn within 1 h was negligible. The net flux of 222 Rn (F net ) can be expressed as: To avoid an overestimation of SGD, the maximum negative value of F net was chosen as the mixing loss [31,61]. As shown in Figure 6, in the dry season, the net 222 Table 1 shows the 222 Rn sources and sinks in the dry and wet seasons. The 222 Rn fluxes of atmospheric loss represented only a small portion of the total 222 Rn sources. The average 222 Rn concentrations of groundwater at four stations (GW1-GW4 in Figure 1) were used to represent the groundwater end-member values at both sites. The 222 Rn concentrations of groundwater ranged from 2850 to 7400 Bq/m 3 (4883 ± 1974 Bq/m 3 ) in the dry season, and from 2868 to 7960 Bq/m 3 (4563 ± 2246 Bq/m 3 ) in the wet season. The 222 Rn flux attributed to the SGD should be divided by the 222 Rn concentration in the groundwater end-member value entering the system. Figure 7 shows the variations in the SGD flux during the time-series measurements at both sites in two seasons. The SGD flux was negatively correlated with tidal height at TD in two seasons, which indicated that the mechanism of tidal pumping was an important process for driving SGD on a short time scale at TD in different seasons. Whereas, the SGD flux was positively correlated with tidal height at HD in different seasons (Figure 7). In the dry season, the SGD flux ranged from 0 to 74.4 cm/d (21.9 ± 18.3 cm/d) at TD and from 0 to 72.5 cm/d (17.8 ± 18.2 cm/d) at HD. In the wet season, the SGD flux varied from 0 to 64.5 cm/d (19.5 ± 17.1 cm/d) at TD, and from 0 to 139.8 cm/d (26.9 ± 29.2 cm/d) at HD. The SGD flux in the wet season at HD was the largest and about 1.5 times that of HD in the dry season. The SGD flux in the wet season at TD was about the same as in the dry season. According to Equation (4) and the data of field measurements, the SGW, Sfw, and Sbay were 30.6, 24.8, and 31.0 ppt in the dry season, respectively, and were 32.3, 29.5, and 33.1 ppt in the wet season, respectively. Therefore, the fresh groundwater fractions in the dry and wet seasons were 7% and 22%, respectively.

Uncertainty Analysis
The estimate of SGD was influenced by multiple uncertain factors, among which the determination of groundwater end-member can make the biggest difference [9,14]. Other factors include the 222 Rn concentrations from sediments, water level, and wind speed, which contribute to smaller SGD fluxes.
In this study, the 222 Rn concentrations of the four groundwater samples were collected from different wetland types (sand beaches, mudflats, tidal marshes, and estuarine intertidal zones). These samples were used as the groundwater end-member value in the mass balance calculation. The mean values and standard deviations of 222 Rn concentration of the groundwater samples were 4883 ± 1974 Bq/m 3 in the dry season and 4563 ± 2246 Bq/m 3 in the wet season. To evaluate the uncertainty in the SGD calculation caused by the groundwater end-member value, a variation in the mean groundwater end-member value of one sigma standard deviation (±40% for the dry season and ± 49% for the wet season) was considered in this study. Therefore, if the groundwater end-member value increases by 40% in the dry season or by 49% in the wet season, the mean SGD rate will decrease by 29% or 33%, respectively. If the groundwater end-member value decreases by 40% in the dry season or by 49% in the wet season, the mean SGD flux will increase by 67% or 96%, respectively.

SGD Fluxes for Different Sites and Seasons
In this study, the SGD rate showed significant differences between the TD and HD sites. The SGD rate in the dry season was larger at TD (21.9 cm/d) than at HD (17.8 cm/d), which was According to Equation (4) and the data of field measurements, the S GW , S fw , and S bay were 30.6, 24.8, and 31.0 ppt in the dry season, respectively, and were 32.3, 29.5, and 33.1 ppt in the wet season, respectively. Therefore, the fresh groundwater fractions in the dry and wet seasons were 7% and 22%, respectively.

Uncertainty Analysis
The estimate of SGD was influenced by multiple uncertain factors, among which the determination of groundwater end-member can make the biggest difference [9,14]. Other factors include the 222 Rn concentrations from sediments, water level, and wind speed, which contribute to smaller SGD fluxes.
In this study, the 222 Rn concentrations of the four groundwater samples were collected from different wetland types (sand beaches, mudflats, tidal marshes, and estuarine intertidal zones). These samples were used as the groundwater end-member value in the mass balance calculation. The mean values and standard deviations of 222 Rn concentration of the groundwater samples were 4883 ± 1974 Bq/m 3 in the dry season and 4563 ± 2246 Bq/m 3 in the wet season. To evaluate the uncertainty in the SGD calculation caused by the groundwater end-member value, a variation in the mean groundwater end-member value of one sigma standard deviation (±40% for the dry season and ± 49% for the wet season) was considered in this study. Therefore, if the groundwater end-member value increases by 40% in the dry season or by 49% in the wet season, the mean SGD rate will decrease by 29% or 33%, respectively. If the groundwater end-member value decreases by 40% in the dry season or by 49% in the wet season, the mean SGD flux will increase by 67% or 96%, respectively.

SGD Fluxes for Different Sites and Seasons
In this study, the SGD rate showed significant differences between the TD and HD sites. The SGD rate in the dry season was larger at TD (21.9 cm/d) than at HD (17.8 cm/d), which was consistent with the tidal range at both sites. Whereas, the SGD flux in the wet season at HD was about 1.4 times that of TD, which was also consistent with the tidal amplitude at two sites. Moreover, the sediment characteristics at both sites are different. The strata around TD is the bedrock of Laoshan Mountain, which has a low permeability. However, the surface sediments in the Jimo Basin, around HD in the northern part of the bay, consist of sand and silt with high permeabilities. This may affect the SGD flux at the two sites [14].
The SGD rates at the TD site did not show large seasonal variations because of the bedrock with low permeability around TD, which is unlikely to be sensitive to rainfall. Figure 7 shows that the SGD flux at TD in two seasons had a negative correlation with water depth during the monitoring period. Therefore, the SGD rate at TD was slightly larger in the dry season than in the wet season, due to large tidal amplitude. Whereas, the SGD rate at HD was larger in the wet season (26.9 cm/d) than in the dry season (17.8 cm/d), which was still consistent with the average tidal range in the two seasons. Figure 2 shows that the total precipitation in Jiaozhou Bay during the spring (March, April, and May 2016) and autumn (September, October, and November 2016) were 87.2 mm and 120.7 mm, respectively. Moreover, the rainfall in the month before the field monitoring periods in spring and autumn were 2.8 mm and 27.9 mm, respectively, which was consistent with the SGD flux at HD in the dry and wet seasons. In addition, the HD site is surrounded by surface sediments of sand and silt with high permeabilities. Also, the northern part of the bay has Quaternary unconsolidated sediments around HD that easily store groundwater, and the sediments were recharged by rainfall. Therefore, the SGD rate at HD was larger in the wet season than in the dry season, due to higher rainfall in the wet season that led to a greater hydraulic gradient between groundwater and seawater [51]. Collectively, the highest SGD rate for both sites was observed at HD in the wet season, likely due to the sedimentary lithology of sand and silt with high permeability, large tidal amplitude, and the rapid response to local precipitation.

Comparisons with Local River Inputs
The product of the SGD rate and the seepage area of Jiaozhou Bay determined the magnitude of the SGD flux into the bay [14]. Detailed information on the water area and tidal flats of the bay, and the average discharge rate of Dagu River was well described in Section 2.1. Tidal flats of sand and silt sediments with a large-scale seepage face mainly occur in the north and northwest of the bay due to a low slope and large tidal amplitudes [41]. To ensure the accuracy of the SGD flux estimation, the mean SGD rate (17.8 cm/d in the dry season and 26.9 cm/d in the wet season) at HD near tidal flats was used for estimating the magnitude of the SGD flux into the bay. Therefore, the SGD flux was about 8-14 times the freshwater inputs from the Dagu River.

Comparisons with Previous Studies
Most previous studies have focused on regional SGD in Jiaozhou Bay [40,41,65]. However, in this study, two fixed sites were selected for estimating the temporal variations of SGD over spatial and temporal scales. The influence of different seasons on the SGD was analyzed, both qualitatively and quantitatively.
In addition, several studies have assessed the SGD flux in Jiaozhou Bay using different methods. Yuan et al. showed that the SGD flux was 9 cm/d in the wet season (September and October 2011) and 4 cm/d in the dry season (April and May 2012) using the radium mass balance models ( 224 Ra and 226 Ra) [65]. Qu et al. used a generalized Darcy's law to evaluate the SGD in four wetland types, including sandy beaches, mudflats, tidal marshes, and estuarine intertidal zones, with a range of 3.6 × 10 −3 to 7.6 cm/d [41]. Zhang et al. applied the 222 Rn mass balance model to estimate the SGD rate (12.3 cm/d in January 2016 and 17.8 cm/d in July 2015) at TD, which was 25% lower than the results found in this study, likely due to the different seasons and the determinations of groundwater end-members [51]. However, compared with the SGD rates in other studies, it is suggested that the SGD in this study is in the range of those in previous studies ( Table 2). The minor differences between the estimations may be a result of different study areas and different times.