Submarine fresh groundwater discharge into Laizhou Bay comparable to the Yellow River flux

Near- and off-shore fresh groundwater resources become increasingly important with the social and economic development in coastal areas. Although large scale (hundreds of km) submarine groundwater discharge (SGD) to the ocean has been shown to be of the same magnitude order as river discharge, submarine fresh groundwater discharge (SFGD) with magnitude comparable to large river discharge is never reported. Here, we proposed a method coupling mass-balance models of water, salt and radium isotopes based on field data of 223Ra, 226Ra and salinity to estimate the SFGD, SGD. By applying the method in Laizhou Bay (a water area of ~6000 km2), we showed that the SFGD and SGD are 0.57 ~ 0.88 times and 7.35 ~ 8.57 times the annual Yellow River flux in August 2012, respectively. The estimate of SFGD ranges from 4.12 × 107 m3/d to 6.36 × 107 m3/d, while SGD ranges from 5.32 × 108 m3/d to 6.20 × 108 m3/d. The proportion of the Yellow River input into Laizhou Bay was less than 14% of the total in August 2012. Our method can be used to estimate SFGD in various coastal waters.

S ubmarine groundwater discharge (SGD) is any and all flows of water on continental margins from the seabed to the coastal ocean, regardless of fluid composition or driving force 1 . SGD, which is driven by both terrestrial and marine forcing components, comprises terrestrial fresh groundwater (SFGD) and re-circulated seawater (RSGD) 1,2 . As an essential part of SGD, SFGD is an important source of freshwater, nutrients, contaminants, and other chemicals to the coastal waters, and has significant impacts and implications on coastal environment and ecology 3,4 . With near half of the global population residing within 100 km of the coastline 5 , it is important to find more freshwater resources including near-and off-shore fresh groundwater for the relief of water scarcity in densely populated coastal megacities 6 . Many previous studies have found that regional scale SGD is several times greater than river discharge [7][8][9] . For example, Kim et al. 7 estimated the magnitude of SGD into the Yellow Sea to be at most 300% of the river water input. Moore et al. 8 showed that the SGD flux is probably between 0.8 and 1.6 times the river flux to the Atlantic. Moore 9 reevaluated the SGD to a large section of the South Atlantic Bight and found the annual average SGD flux is three times greater than the river fluxes. However, there are few studies on large scale SFGD as compared with those on SGD, and SFGD with magnitude comparable to large river discharge is never reported. In this paper, we proposed a method for combining mass-balance models of water, salt and radium isotopes to estimate the SFGD, SGD. By applying the method in Laizhou Bay, the SFGD, SGD, and the proportion of the Yellow River input into Laizhou Bay were estimated and discussed.
Laizhou Bay, located between 37.05uN , 37.80uN and 118.9uE , 120.35uE (WGS84 reference system), is one of the three major bays in the Bohai Sea, China (Fig. 1). The natural coastline of the Bay is relatively straight and extends from the Qimu Cape to the Yellow River Estuary. Along the coastline, there are at least four different types of depositional environments: the Yellow River delta in the northwest, the alluvial plain in the southwest, the marine deposit plain in the south, and the hilly area in the east [10][11][12] . The coast areas west of Hutouya ( Fig. 1) are alluvial or marine deposit plains with aquifers mainly composed of permeable coarse material. The coast area east of Hutouya is hilly with a coastal plain. Groundwater occurs mainly in the Quaternary aquifers of the coastal plains. The average annual rainfall is ,640 mm and occurs mostly from June to September, while the potential evaporation is approximately 2050 mm 13 . The seawater depth of the Bay generally increases as the offshore distance increases, and has an average depth of ,8 m. Based on salinity measurements at surface, middle and bottom at the stations S2-S6 shown in Fig. 1, where the seawater is deepest in whole Laizhou Bay, it was found that the Bay is vertically well mixed and can be treated as a single layer in terms of the large-scale water motion.
There are many rivers flowing into the Bay, including the Yellow River, Jiaolai River, Xiaoqing River and Wei River. Particularly, the Yellow River, as the second and sixth largest river in China and the world, respectively, is the largest that discharges into the Bohai Sea. However, the Yellow River discharge has significantly decreased since the 1950s due to both climate change and human activities 14,15 . Seasonal variations of the Yellow River discharge are also significant, with the minimum occurring usually in April and May, and the maximum occurring usually in July and August (Fig. 2). The annual average discharge of the Yellow River in 2012 was 7.23 3 10 7 m 3 /d. In July and August 2012, its discharge rate was at least 50 times greater than the sum of the discharge rates of all other rivers flowing into Laizhou Bay.
The Yellow River deposits large amounts of sediments, creating a fast-growing delta area. Its course into the Bohai Sea changed a number of times over the past several decades 16 . The present channel through the delta was formed artificially in 1996. The mouth of the Yellow River is located between Bohai Bay and Laizhou Bay, and the shifting of the mouth to the north may significantly reduce its direct discharge into Laizhou Bay.
As one of the three bays in the Bohai Sea, China, and with a shoreline of ,320 km and area of ,6000 km 2 , Laizhou Bay is an important coastal environment. As well, it is subject to a variety of environmental stresses, and hence provides an archetype of a semienclosed bay for which ecological functioning is a sensitive issue. Most existing studies on SGD in the Bohai Sea are restricted to the Yellow River delta [17][18][19][20][21] , the largest estuary in the Bohai Sea. In order to estimate SFGD, SGD, and the Yellow River input into Laizhou Bay, activities of radium isotopes ( 223, 226 Ra) and salinities were measured (Fig. 1). Natural radium isotopes are ideal tracers for effective and efficient assessment of SGD since they are conservative chemically and widely enriched in groundwater relative to surface waters 7,8,21,22 . In August 2012, we collected eight groundwater samples, six river water samples, and 44 seawater samples for measurement of radium isotopes ( 223, 226 Ra) and salinity (Supplementary Table S1). In May 2014, we measured the salinities at the same locations.

Results
Spatial distribution of radium isotopes and salinity. Figure 3 shows the spatial distributions of the two radium isotopes ( 223, 226 Ra) and salinity within Laizhou Bay. The activity distributions of these two isotopes have the following common features: (1) the activities were significantly higher in the west and south than in the east and north of the Bay; and (2) the activities were very high in the estuary and nearshore areas and they generally decreased with the offshore distance (Figs. 3a and 3b). The seawater salinities in August 2012 were higher in the east than in the west of the Bay (Fig. 3c), where the activities of the radium isotopes were comparatively high. The seawater salinities in May 2014 were similar to those in August 2012 overall, but with a local low-salinity area near the middle of the eastern coast (Fig. 3d), indicating considerable SFGD since there are no river inputs in that region.
Determination of flushing time. To explore the dynamics of coastal water, one should first estimate the flushing time T f [T] for a bay, i.e., the ratio of the mass or volume of a constituent (V) to its renewal rate (Q) 23 , or T f 5 V/Q. Combining the radium isotope method 24 and the tidal prism method 25 , we obtained the following equation to estimate the flushing time T f based on measurements of radium isotopes, where F( 223 Ra/ 226 Ra) is the 223 Ra/ 226 Ra activity ratio of the flux into the bay, I( 223 Ra/ 226 Ra) is the 223 Ra/ 226 Ra activity ratio in the bay, and l 223 5 0.061d 21 is the decay constant of 223 Ra. Using equation (1), we estimated T f 5 36.      Assessment of SFGD and SGD. To quantify the freshwater fluxes into Laizhou Bay, we developed water-mass and salt-mass balance models that include seawater from outside the Bay, river input, SFGD, precipitation (P T ) and evaporation (E T ), all estimated over the period of T f d immediately before the observation time (taken as the mid-point of the 8-d field sampling period).
Coupling the salt-mass balance model with the water-mass balance model, we derived the SFGD flux (Q SFGD ) as (see Water and Salt Mass Balance Model in Supplementary Information): where r YL is the ratio of the Yellow River input into Laizhou Bay to the total input into the Bohai Sea (Q Y ) during the flushing time immediately before the observation time; Q r,i is runoff of the i th river other than the Yellow River into the Bay (Supplementary where 226 Ra TYL is the total input of 226 Ra from the Yellow River discharge and suspended particle desorption; 226 Ra Tr is the total input of 226 Ra from discharges and suspended particle desorption of all rivers other than the Yellow River; 226 Ra gw and 226 Ra Bay are the radium activity in groundwater and Bay water, respectively; and 226 Ra bg is the background activity. Detailed calculations of these parameter values are described and summarized in Supplementary  Table S4. Figure 4a shows how the SFGD and SGD fluxes change with the proportion of the Yellow River input into Laizhou Bay (r YL ). The annual average flux Q Y2012 57.23 3 10 7 m 3 /d of the Yellow River in 2012 was used as the reference, which is more than 20 times the fluxes of all the other rivers flowing into Laizhou Bay. One can see that both SFGD and SGD decrease as r YL increases. Note that SGD is much greater than SFGD since SGD includes RSGD. From Fig. 4a, since SFGD $ 0, we have SGD $ 5.13Q Y2012 and r YL # 0.396. When r YL 5 0, we obtain the upper-bound estimate of SGD as 8.57Q Y2012 and the corresponding value of SFGD as 0.88Q Y2012 . In general, we have SGD 5 (8.57 2 8.69r YL )Q Y2012 , and SFGD 5 (0.88 2 2.22r YL )Q Y2012 .
Seasonal variations of the Yellow River flux were significant. The minimum flux occurred in April-May and was less than one-tenth of the maximum flux in July-August (Fig. 2) Figure 4b shows how the SFGD predicted by equation (2) using the salinity data observed in May 2014 changes with r YL , the proportion of the Yellow River input into Laizhou Bay, dur-ing April and May in 2014. As r YL increases from 0 to 1, the SFGD decreases from 0.79Q Y2012 to 0.34Q Y2012 . As expected, the impact of the Yellow River flux on the SFGD was significantly reduced by comparison with that of August 2012 (Fig. 4a).
The precipitation in Laizhou Bay is much larger in June-August than in March-May, and so the groundwater table in the coastal area should be higher (i.e., higher hydraulic head) in August than in May owing to rainfall infiltration into ground surface. Due to the increased landward gradient in hydraulic head, it follows that SFGD (August 2012) $ SFGD (May 2014). Using the lower-bound estimate for SFGD in May 2014 (0.34Q Y2012 ), from Fig. 4a one can see that if SFGD in August 2012 $ 0.34Q Y2012 , then r YL (August 2012) was less than 0.24 and SGD (August 2012) $6.48Q Y2012 (as is indicated by the thick red-dashed vertical line in Fig. 4a).
Although we do not have direct evidence to conclude that r YL (May 2014) #0.24 since r YL varies with time, it should be considerably less than unity. In addition, the Yellow River mouth, which is located between Bohai Bay and Laizhou Bay, is northerly-oriented, as shown in photographs taken by NASA Landsat in May 2014 or Fig. 2 of Xu et al. 21 . This fact may significantly limit the input of the Yellow River into Laizhou Bay. Thus, it is at least not unreasonable to assume r YL (May 2014) #0.50. From Fig. 4b one can see that r YL (May 2014) #0.50 implies SFGD (May 2014) $0.57Q Y2012 . Since SFGD (August  Fig. 4a one can see that SFGD in August 2012 was in the range (0.57 , 0.88)Q Y2012 , and r YL (August 2012) was less than 0.14 and SGD in August 2012 was between 7.35Q Y2012 (as is indicated by the thin red-dashed vertical line in Fig. 4a) and 8.57Q Y2012 .
The seawater circulation in the Bohai Sea shows apparent seasonal variations, which in turn affect the path of the Yellow River discharge 28 . During the summer months, the monsoonal winds blow from the south in this region, thus creating a cyclonic gyre within the Bohai Sea 26,28 . These may limit the input of the Yellow River into Laizhou Bay and support the small estimated value of r YL (August 2012).

Discussion
Our tracers-orientated estimation of SGD is approximately ten times as large as SFGD. SGD is much greater than SFGD since SGD5 SFGD 1 RSGD, and many previous studies show that tides, waves, fluid density gradient, storms, geothermal gradient, and seabed topography can result in RSGD [29][30][31][32][33][34][35][36] . As a result, the SGD values estimated by the tracers-orientated method and mechanism-orientated method such as traditional hydrogeological approach often do not match. Li and Jiao 29 reviewed the studies of tidal contribution to SGD and found that the RSGD induced by tides in the intertidal zone is at least 1 , 2 orders of magnitude smaller than SGD estimated by radium isotope tracers. The order of magnitude of RSGD contributed by waves is the same as that of tides 30,31 . The density-driven RSGD is usually much less than SFGD [32][33][34] . The RSGD induced by geothermal gradient and seabed topography is at least 4 orders of magnitude smaller than the SGD estimated by radium isotope tracers [35][36][37] . Thus, SGD estimated by radium isotope tracers is much larger than that of the total sum of SGD induced by all the above factors/ mechanisms. Solving this challenging problem needs long-term efforts involving close combination of tracers-orientated and mechanism-orientated methods.
We can compare our results with small-scale SGD studies. Taniguchi et al. 19  Based on their results, the average ratio of SFGD to SGD is less than 1% in both sampling periods. Our results, if using the same unit (treating SGD as a line source along the coastline and dividing the total SGD by the shoreline length of the Bay), ranges from 1662.2 m 3 /m/d to 1937.2 m 3 /m/d, which are slightly smaller than Taniguchi's estimates. The ratio of SFGD to SGD in our study, ranging from 7.5% to 9% as shown in Fig. 4a when r YL # 0.14, however, is much higher than that of their study. Peterson et al. 18 estimated SGD in the same area as Taniguchi et al. 19 using radon and radium isotopes. They estimated a SGD flux of 4.5-13.9 cm/d in September 2006, most of which was recirculated seawater. Our results, if using the same unit (treating SGD as an area source and dividing the total SGD by the area of the Bay), ranges from 8.9 cm/d to 10.3 cm/d, which are within the range of Peterson's estimates for the same season. Since our tracers-orientated estimations of SGD and SFGD are based on the mass balance of radium isotopes and fresh groundwater after they entered the seawater, they are conservative in the sense that if there had no groundwater exploitations in the coastal area of Laizhou Bay 12,38 , their values would be larger.
The Bohai Sea has three major bays: Laizhou Bay, Bohai Bay and Liaodong Bay (Fig. 1). The seawater salinity in Laizhou Bay was significantly lower than that in other two bays 39 . The real reason for this was not investigated in detail because it seems consistent with the freshwater discharge from the Yellow River abutting Laizhou Bay. Our investigation gives a plausible explanation for the abnormally low salinity in Laizhou Bay. We concluded that the SFGD, which accounts for 57% , 88% of the annual flux of the Yellow River in 2012, is a key contributor to the low salinity given that the direct input from the Yellow River to Laizhou Bay is very limited (r YL # 14% in August 2012).
SGD has been widely recognized to be a pathway for enriching coastal waters in nutrients, carbon and metals 1,8,22 . In some areas, nutrient fluxes via SGD were shown to rival those from surface waters 7,40,41 . The nutrient input via SGD in the Yellow River Estuary is at least five times of that via the Yellow River 20 . With new understanding of our assessments of the SFGD, SGD and the Yellow River input into Laizhou Bay, the management of the Bay related to fresh groundwater resources, ecology and environment in coastal and offshore areas should be reviewed.
Although SGD has been estimated in many coastal areas all over the world based on radium isotope tracer methods 29,42 , large scale tracers-orientated SFGD studies have not been correspondingly conducted even if SFGD is as important as, or even more important than SGD. In this study, we proposed a tracers-orientated method (using radium isotopes and salinity in seawater and coastal groundwater as tracers) to estimate SFGD. To the knowledge of the authors, this is the first time to quantify large-scale SFGD using a tracers-orientated method based on field radium and salinity measurements. Our method has potential for application considering that the seawater salinity is a common physical quantity that can be measured easily and a straightforward tracer to indicate freshwater. The proposed method can be readily applied to estimate SFGD in other coastal areas all over the world.
The main limitation of current study is most probably the steady state premise, a common approach used in all the previous tracersorientated SGD studies by radium isotope methods 7,8,22,24,43 . Since the variation of radium storage in the whole bay approaches zero near the time when the radium mass in the whole bay reaches its maximum, and the observation period of our field work is just near such a maximum-time, the steady state assumption is approximately valid. Strict quantification of the error induced by steady state assumption, however, needs not only much more radium isotope measurement data in terms of time series, but also quantification of the seawater flow in the whole bay for a long period. This will be an interesting and challenging work in the future.

Methods
Sampling. Radium samples were collected from and adjacent to Laizhou Bay. Large volume water samples (,60 L for seawater, ,15 L for river water, respectively) were pumped and filtered through a 0.45-mm filter for Ra extraction. The coastal groundwater samples (,15 L) were taken from the nearshore zone within 100 m of the high tide mark with PushPoint samplers inserted into sediments at a depth of ,1.5 m. These samples can capture the chemical components (particularly radium isotopes) of groundwater immediately before it discharges into the sea from the aquifer. The water samples were passed slowly through Mn-fibers (,25 g) produced according to the method proposed by Moore 44 for extracting radium isotopes. The flow rate was controlled not to exceed 1 L/min to ensure complete Ra adsorption on the Mn-fiber. These fibers were then washed thoroughly to remove all particles and salts and taken to the hydrogeological laboratory at The University of Hong Kong for measurements. The salinity, temperature, pH of the water samples were measured in situ using a HI9828 Model probe (HANNA).
Measurements. The long-lived radium isotope, 226 Ra, was determined by a radon-inair monitor (RAD7, Durridge Co.) as proposed by Kim et al. 45 . After the 223 Ra measurements were completed, the Mn-fiber samples were aged for 2 -6 weeks to allow 222 Rn and its daughters to equilibrate with 226 Ra. This is an indirect measurement of 226 Ra based on secular equilibrium between 226 Ra and 222 Rn. In order to improve the test results, the determination of 226 Ra in our study slightly modified the method of Kim et al. 45 . Before measurement we made long fiber incubation time (. 20 d) to ensure secular equilibrium between 226 Ra and 222 Rn and then allowed ,24 h for purging the RAD7 system prior to analysis thereby reducing the background noise. In addition, we reduced errors by using long measurement times 46 . The expected error of 226 Ra measurements is 67%. The short-lived radium isotope 223 Ra was analyzed using a two-channel radium delayed coincidence counting system (RaDeCC) 47 . The expected error of 223 Ra measurements is 612%.
Models. Three models (flushing time, water and salt mass balance models) as well as the 226 Ra mass balance model were used in the study to estimate the flushing time of www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 8814 | DOI: 10.1038/srep08814 coastal water, the freshwater discharge, the Yellow River input, and SGD, respectively. Details were shown in Supplementary Information.