This is a non-peer reviewer preprint submitted to EarthArxiv. The manuscript was submitted for review to Nature Communications on 28 march 2019. The flow of fresh groundwater and solutes to the world’s oceans and coastal ecosystems

The flow of fresh groundwater towards the world oceans may provide substantial inputs of critical nutrients and solutes to the oceans. Here we present the first spatially resolved global model of coastal groundwater discharge to show that, in contrast to most previous estimates, the contribution of fresh groundwater accounts for only 0.2% of the freshwater input and 1% of the solute input to the oceans. However, the coastal discharge of fresh groundwater displays a high spatial variability and for an estimated 20% of the world’s estuaries, salt marshes and coral reefs the flux of terrestrial groundwater is high enough to pose a risk for pollution and eutrophication. Coastal groundwater discharge constitutes a diffuse and largely unmonitored source of nutrients which given the slow rates of groundwater flow may continue to pose a pollution risk for vulnerable coastal ecosystems long after sources of pollution have been addressed. Submarine groundwater discharge, the flow of fresh or saline groundwater to oceans, may be a significant contributor to the water and chemical budgets of the world's oceans. The fresh component of submarine groundwater discharge is critical due to its high solute and nutrient load, has been estimated to be up to 10% of the river discharge to the world’s oceans and to equal the inputs by rivers to the for solutes such as carbon, iron, silica and strontium, and potentially buffers ocean acidification with groundwater alkalinity. Although a large number of studies have estimated fresh submarine groundwater discharge (fresh SGD) locally, it is difficult to derive a global estimate from these point estimates, because they are highly variable, often uncertain and strongly biased towards high discharges. Global-scale estimates of fresh SGD vary by four orders of magnitude, with the most recent estimates ranging up to 10% of the river discharge towards the oceans. Large-scale spatially distributed estimates have only been published for the entire or parts of the US coast. The total flux of groundwater to the ocean (herein called coastal groundwater discharge) can be divided into three distinct fluxes (Fig. 1a): 1) fresh SGD that we consider to be meteoric groundwater discharging below the mean sea level; 2) near-shore terrestrial groundwater discharge (NGD) that is fresh groundwater discharging above the mean sea level in the first hundreds of meters near the coastline; and 3) recirculated sea water. Fresh SGD and NGD are driven by recharge from onshore precipitation and are critical due to their high solute and nutrient load whereas recirculated sea water is driven by waves, tides, storm surges and density-dependent flow. We use density-dependent groundwater models to quantify the partitioning of onshore and offshore groundwater discharge and the sensitivity of coastal discharge to controlling variables such as topography, permeability, recharge and size of coastal watersheds. We subsequently quantify coastal groundwater discharge at a global scale by combining a series of numerical models of coastal groundwater flow with a global geospatial analysis of controlling variables. In addition to submarine groundwater discharge we quantify near-shore terrestrial discharge, a flux that has been overlooked in most hydrological studies but affects coastal water budgets, evapotranspiration and ecosystems. Controls on coastal groundwater discharge Sensitivity analysis using a density-dependent model of coastal groundwater flow shows that fresh SGD is an insignificant part of the total groundwater flux in most coastal groundwater systems, and most groundwater discharges on land (Fig. 1d-g) where it contributes to baseflow to rivers and evapotranspiration. For a coastal groundwater system consisting of relatively permeable rocks (Fig. 1b) 29 % of the total groundwater recharged on land discharges in the oceans. However, for a system with the global median values of permeability and topographic gradient which were derived from a global geospatial analysis of coastal watersheds (see supplementary information section S1) the modeled fresh SGD flux is only 0.5% of the total recharge flux. The discharge flux reaches a maximum of 0.23 m a at the coastline (Supplementary Fig. S13), which is lower than the lowest local value of fresh SGD reported in the literature known to us. While fresh SGD is insignificant in most cases, when permeability exceeds a threshold value of 10 m and topographic gradients exceed 1%, fresh SGD increases rapidly and can become the dominant groundwater discharge component (Fig. 1e and 1g). Somewhat counterintuitively, groundwater recharge and size of the contributing area do not influence SGD in most cases. Once a threshold value is reached where groundwater input equals the flow capacity of the subsurface any further increases in recharge or contributing area only generate additional terrestrial discharge, whereas fresh SGD stays constant (Fig. 1d and 1f). The model results show that offshore freshwater discharge is mirrored by a zone of nearshore terrestrial groundwater discharge (NGD) (Fig. 1b), which is higher than fresh SGD in most settings except at high permeability or topographic gradients (Fig. 1d-g), but has to our knowledge been overlooked in global hydrological assessments. The existence of onshore groundwater seepage near surface water features has long been recognized for lakes, wetland or streams. However, seepage near coastlines has not been quantified at large scales. In the two-dimensional crosssectional models presented here near-shore discharge generates a seepage flux across the land surface that peaks at the coastline and decays exponentially with distance to the coast (Fig. 1b and SI Fig. S13). This modeled seepage flux represents a mixture of evapotranspiration, ponding, surface runoff and lateral groundwater flow perpendicular to the cross-section towards streams.

Submarine groundwater discharge, the flow of fresh or saline groundwater to oceans, may be a significant contributor to the water and chemical budgets of the world's oceans 1,2 . The fresh component of submarine groundwater discharge is critical due to its high solute and nutrient load 3 , has been estimated to be up to 10% of the river discharge to the world's oceans 1 and to equal the inputs by rivers to the for solutes such as carbon 4 , iron 5 , silica 6,7 and strontium 8 , and potentially buffers ocean acidification with groundwater alkalinity 9 . Although a large number of studies have estimated fresh submarine groundwater discharge (fresh SGD) locally, it is difficult to derive a global estimate from these point estimates, because they are highly variable, often uncertain and strongly biased towards high discharges 1 . Global-scale estimates of fresh SGD vary by four orders of magnitude, with the most recent estimates ranging up to 10% of the river discharge towards the oceans 1,2,10 . Large-scale spatially distributed estimates have only been published for the entire 11 or parts of the US coast 12 . The total flux of groundwater to the ocean (herein called coastal groundwater discharge) can be divided into three distinct fluxes (Fig. 1a): 1) fresh SGD that we consider to be meteoric groundwater discharging below the mean sea level; 2) near-shore terrestrial groundwater discharge (NGD) that is fresh groundwater discharging above the mean sea level in the first hundreds of meters near the coastline; and 3) recirculated sea water. Fresh SGD and NGD are driven by recharge from onshore precipitation and are critical due to their high solute and nutrient load whereas recirculated sea water is driven by waves, tides, storm surges and density-dependent flow.
We use density-dependent groundwater models to quantify the partitioning of onshore and offshore groundwater discharge and the sensitivity of coastal discharge to controlling variables such as topography, permeability, recharge and size of coastal watersheds. We subsequently quantify coastal groundwater discharge at a global scale by combining a series of numerical models of coastal groundwater flow with a global geospatial analysis of controlling variables. In addition to submarine groundwater discharge we quantify near-shore terrestrial discharge, a flux that has been overlooked in most hydrological studies but affects coastal water budgets, evapotranspiration and ecosystems.

Controls on coastal groundwater discharge
Sensitivity analysis using a density-dependent model of coastal groundwater flow shows that fresh SGD is an insignificant part of the total groundwater flux in most coastal groundwater systems, and most groundwater discharges on land ( Fig. 1d-g) where it contributes to baseflow to rivers and evapotranspiration. For a coastal groundwater system consisting of relatively permeable rocks (Fig.   1b) 29 % of the total groundwater recharged on land discharges in the oceans. However, for a system with the global median values of permeability and topographic gradient which were derived from a global geospatial analysis of coastal watersheds (see supplementary information section S1) the modeled fresh SGD flux is only 0.5% of the total recharge flux. The discharge flux reaches a maximum of 0.23 m a -1 at the coastline ( Supplementary Fig. S13), which is lower than the lowest local value of fresh SGD reported in the literature known to us 13 . While fresh SGD is insignificant in most cases, when permeability exceeds a threshold value of 10 -12 m 2 and topographic gradients exceed 1%, fresh SGD increases rapidly and can become the dominant groundwater discharge component ( Fig. 1e and   1g). Somewhat counterintuitively, groundwater recharge and size of the contributing area do not influence SGD in most cases. Once a threshold value is reached where groundwater input equals the flow capacity of the subsurface any further increases in recharge or contributing area only generate additional terrestrial discharge, whereas fresh SGD stays constant ( Fig. 1d and 1f).
The model results show that offshore freshwater discharge is mirrored by a zone of nearshore terrestrial groundwater discharge (NGD) (Fig. 1b), which is higher than fresh SGD in most settings except at high permeability or topographic gradients ( Fig. 1d-g), but has to our knowledge been overlooked in global hydrological assessments. The existence of onshore groundwater seepage near surface water features has long been recognized for lakes, wetland or streams 14,15 . However, seepage near coastlines has not been quantified at large scales. In the two-dimensional crosssectional models presented here near-shore discharge generates a seepage flux across the land surface that peaks at the coastline and decays exponentially with distance to the coast ( Fig. 1b and SI   Fig. S13). This modeled seepage flux represents a mixture of evapotranspiration, ponding, surface runoff and lateral groundwater flow perpendicular to the cross-section towards streams.  1. Sensitivity analysis (d-g) of numerical models of coastal groundwater discharge (a-b) demonstrates that the flux of groundwater to the ocean is controlled primarily by topographic gradient (e) and aquifer permeability (g) and is relatively insensitive to watershed length (d) and groundwater recharge (f). The histograms in panels d-f show the distribution of the controlling parameters in all 40,082 global coastal watersheds. Panel b and c show modeled groundwater fluxes over the land surface and seabed (b), flowlines and salinity (c) for a base-case model coastal watershed with global median values for groundwater recharge (0.143 m a -1 ), but with relatively high permeability (10 -12 m 2 ) and topographic gradient (2.5%).
Global coastal groundwater discharge, fresh SGD and NGD were calculated by pairing the results of 351 model runs to geospatial data of 40,082 coastal watersheds (Fig. 2, S8 and S9). The results show that fresh SGD is an insignificant contributor to the water and solute budgets of the world's oceans ( Supplementary Fig. S8). Only 4 (0.04-11) % of the groundwater recharged in coastal watersheds is discharged offshore. The calculated global fresh submarine groundwater discharge is 78 (0.4-210) km 3 a -1 , which amounts to 0.2 (0.001 -0.6)% of the river discharge to the oceans 16 and which is at the lower end of previous global estimates 1,10 . The uncertainty of the global fresh SGD flux is mostly caused by the high uncertainty of the values of permeability that were used. Additional sources of uncertainty are the representative topographic gradient of coastal watersheds, groundwater recharge and the size the area that contributes to fresh SGD. Fresh SGD would be 44 instead of 78 km 3 a -1 when a lower estimate of topographic gradients would be adopted that follows the average gradient of coastal streams instead of the average gradient. The difference between two alternative global models of groundwater recharge results in an uncertainty range of 53 to 78 km 3 a -1 .
When the contributing area is assumed to be twice the size of coastal watersheds the fresh SGD flux is 118 km 3 a -1 . While the uncertainty ranges reported here are large, the best estimates of permeability, topographic gradient and recharge provide a good fit of the model to observed watertable gradients in coastal watersheds ( Supplementary Fig. S11), which suggests that the reported best estimates of coastal discharge fluxes are relatively robust.
The calculated fresh SGD is only a minor fraction, 0.06% (0.0003% -0.2%), of the global total SGD flux, which includes recirculated seawater, and has been estimated as three to four times of the river flux globally based on measured concentrations of radiogenic radon in seawater 17 . This is in line with analytical models that estimated seawater circulation due to tidal and wave forcing to be roughly equal to the estimated total SGD flux 18 . The very low terrestrial contribution to the overall SGD flux means that global SGD consists almost exclusively of recirculated seawater and the netinput of solutes to the oceans by SGD is much lower than previously assumed. A first order estimate of the input by fresh SGD of carbon, nitrogen, silica and strontium based on published compilations of the average solute concentrations in coastal groundwater 7,8,19,20 suggest the contribution of submarine groundwater discharge is ~1% of the input by rivers (Table 1), which is much smaller than the up to 100% contribution suggested by some recent studies that extrapolated global inputs from local and regional-scale estimates [4][5][6] . There are insufficient data available for the concentration of iron in coastal groundwater, but even with relatively high concentrations of 40 mg L -1 that exceed local estimates in the literature 5 the contribution would equal 1% of the river flux 21 . Our estimate of solute fluxes in fresh SGD is a first order, spatially aggregated average that assumes conservative transport. The spatially distributed map of coastal groundwater discharge (Fig. 3) provides the possibility of better future estimates for locations where concentrations of specific solutes are known.
In addition to submarine groundwater discharge, the model results suggest that a substantial part of coastal groundwater discharge takes place onshore ( Supplementary Fig. S9). Global NGD is estimated as 147 (1 -289) km 3 a -1 , which is 0.4 (0.001-0.8) % of the global river discharge. NGD takes place in a zone that on average extends 400 m from the shore (Supplementary Fig. S15). NGD predominantly contributes to evapotranspiration, but exceeds potential evapotranspiration 22 in 28 (0.07-58) % of the global coastline where it contributes to surface runoff and baseflow close to the shoreline. The combined submarine and near-shore groundwater discharge flux equals 224 (1.4 -500) km 3 a -1 and is shown in Fig. 3a. Model experiments demonstrate that the partitioning of coastal groundwater discharge in submarine and terrestrial discharge is highly sensitive to the local topographic gradient (Fig. 1d), and therefore the total coastal groundwater discharge is a more robust estimate than the onshore and offshore components.   3. The flow of groundwater to the oceans is highly variable, dominated by localized hotspots, and can locally be source of water that rivals surface water discharge. Global maps of a) coastal groundwater discharge (CGD) in m 2 a -1 and b) coastal groundwater discharge (CGD) divided by the surface water flux to the oceans.

Hotspots of coastal groundwater discharge
Although the overall contribution of water and solutes by coastal groundwater discharge to the oceans is low, coastal groundwater discharge is highly variable, with 10% of the global coastline contributing 90% of the total discharge. Comparison with data on the discharge of rivers to the oceans 23 shows that coastal groundwater discharge is highly variable and can locally be close to the river input ( Fig. 3a and 3b), which given its relatively high solute load 3 means that it may dominate the solute input in some coastal ecosystems. We define coastal groundwater discharge hotspots as watersheds where the coastal groundwater discharge exceeds 100 m 2 a -1 and 25% of the river discharge. The threshold value of 100 m 2 a -1 reflects a relatively conservative lower bound for reported values at locations with high coastal groundwater discharge (Supplementary Tables S5, S6) and associated ecosystem impacts (Supplementary Table S7). Coastal groundwater discharge hotspots ( Fig. 3a and 3b) cover 9 (0.02-30) % of the global coastline and are predominantly located in areas with a steep coastal topography due to glacio-isostatic rebound, active tectonics or volcanic activity and areas consisting permeable unconsolidated sediments, carbonates or volcanic rocks ( Supplementary Fig. S16). The distribution of hotspots is consistent with documented sites of high fresh groundwater discharge globally that are predominantly located in North America, Europe and East Asia (Supplementary Table S6 and S7). However, at many hotspots, such as Iceland and parts of South America, Africa and south Asia, and many tropical islands coastal groundwater discharge has been unexplored to our knowledge.

Local impacts of coastal groundwater discharge
Due to its high spatial variability, fresh groundwater discharge can locally have strong impacts on coastal hydrology and ecosystems. Coastal groundwater discharge can control the salinity, nutrient budget and productivity of coastal lagoons 24 , salt marshes 25 [34][35][36] shows that 26% (0.4%-39%) of the world's estuaries and 17% (0.3%-31%) of the salt marshes at risk of eutrophication. In addition, 14% (0.1-26%) of the coastline that is located within 500 m of a coral reef is at risk of eutrophication. A review of sites with documented ecosystem impacts of coastal groundwater discharge (Supplementary Table S7) suggests that the threshold values of coastal groundwater discharge and nitrogen input in the adjacent coastal watersheds that we used to define high risks are relatively conservative and that adverse ecosystem impacts may also occur at lower threshold values for nitrogen application or coastal groundwater discharge. In addition to its importance for coastal ecosystems, coastal groundwater discharge can locally also be a freshwater resource that is used as drinking water or for other purposes in a limited number of locations, but has been generally overlooked 37 . However, this resource may also be sensitive to pollution, and exploitation of this resource would need to be carefully managed to avoid salt-water intrusion and adverse impacts on coastal ecosystems.
The values of global coastal groundwater discharge represent a natural undisturbed system and do not include groundwater pumping. However, many coastal groundwater systems may be affected by groundwater pumping. A comparison of the calculated coastal groundwater flux and published values of groundwater depletion in coastal watersheds 38 shows that in most coastal groundwater systems are not associated with depletion (Fig. 4b). Groundwater depletion is concentrated in semi-arid regions and in 13% (8% -19%) of the global coastal watersheds groundwater depletion exceeds coastal groundwater discharge. In these watersheds coastal groundwater discharge has already or will reduce to zero at some stage in the future and instead seawater will start intruding into terrestrial groundwater systems.

Comparison with published large-scale models of coastal groundwater discharge
Our model results estimate a coastal groundwater discharge flux for the contiguous US of 8.5 (0.1-17) km 3 a -1 , which is in the same range as the 15±4 km 3 a -1 estimated by Sawyer et al. 11 . However, the two estimates cannot be compared directly. The estimate by Sawyer et al. is based on the assumptions that groundwater discharge is controlled solely by surface morphology and drainage density, which is a debatable assumption. In contrast our model experiments highlight the role of the flow capacity of the subsurface and the key role of permeability in governing groundwater flow and discharge in coastal groundwater systems ( Fig. 1 and 2). For the Gulf and Atlantic coasts of the USA our model estimates a coastal discharge of 6.6 (0.1-11.3) km 3 a -1 , which is at the lower end of recent estimates by Zhou et al. 39 that range from 9.7 (7.2-12.0) to 27.1 (22.8-30.5) km 3 a -1 . These estimates are based on a series of regional groundwater models 12 , include more detailed hydrogeology and permeability structure and may therefore be more accurate than our model estimates, especially at local scales. On the other hand, these models did not include solute transport and density-dependent flow and used a discretization of 250 m. As a result, the partitioning of groundwater discharge around the coastline may not have been well resolved, and not including the fresh-salt water interface may have led to overestimation of coastal groundwater discharge. Both methods are difficult to apply to global scales. Methods based on surface morphology and hydrography depend critically on the availability of high resolution datasets that are currently only available for part of the global coastline. While there have been recent advances towards global 3D groundwater models 40 , the prohibitive computational expense of including density-riven flow in these models limits their utility for resolving coastal groundwater fluxes.

Conclusions
The assessment of near-shore submarine and terrestrial groundwater discharge reported here provides the first estimate at high spatial resolution of these water fluxes at a global scale. Our analysis shows that the global flux is dominated by a small number of coastal watersheds that are distributed around the globe, including numerous locations that have so far not been explored. In contrast to river discharge, coastal groundwater discharge is frequently unmonitored. However, our global analysis shows that locally coastal groundwater discharge can in many cases pose an equally serious risk for coastal water quality and coastal ecosystems. In addition, groundwater discharge is relatively diffuse compared to surface water discharge and may therefore affect larger areas. Coastal groundwater discharge links onshore groundwater systems with coastal ecosystems, which means that changes in groundwater pumping or land use affect the flow of nutrients to coastal ecosystems 41,42 , and should be taken into account in coastal environmental management. The estimates provided here can help guide future research and monitoring of this water flux and its effect on coastal ecosystems. This is especially important as population pressures and increase in agricultural activity are likely to increase nutrient and contaminant inputs to coastal groundwater in many areas in the future. Furthermore, because groundwater flow rates are typically very slow, measures to improve groundwater quality onshore may take decades before they affect offshore water quality 43 .
Therefore, quantitative estimates of coastal groundwater discharge are of key importance for identifying present and future risks to coastal water quality.

Modeling coastal groundwater discharge
We simulated submarine and terrestrial groundwater discharge in coastal groundwater systems using a numerical model of coupled density-driven groundwater flow and solute transport. The model code, GroMPy-couple, is a Python shell around the finite element code escript 44,45 which has been used to simulate subsurface fluid flow 44 . We implemented an iterative scheme 46 to solve the fluid flow and solute transport equations and the equations of state for fluid density and viscosity.
GroMPy-couple uses an implementation of a seepage boundary condition based on an existing implementation in the model code MODFLOW 15 , which ensures a realistic and numerically stable partitioning between onshore and offshore groundwater discharge in models of coastal groundwater systems 15 . The model code simulates the flow of fresh (meteoric) groundwater, the mixing and recirculation of seawater at the fresh-salt water interface at the coast due to dispersion and the onshore and offshore discharge of groundwater. We did not model transient processes like tidal forcing of groundwater flow and wave setup, that are responsible for the bulk of the recirculation of seawater in coastal aquifers 47,48 . Note that the submarine discharge of fresh groundwater is relatively insensitive to transient seawater circulation offshore 48 . The model code has been validated by comparison with a salt water intrusion experiment 49 , analytical solutions of groundwater discharge 14,15 , and model experiments using the widely-used model code SUTRA 50 . See Supplementary Information section S2 for more details of the model approach.

Geospatial data analysis
The model input was based on a geospatial data analysis of the controlling parameters of coastal groundwater discharge. We analyzed watershed geometry 51 , topographic gradients 52 , permeability 53,54 and groundwater recharge 23,55 for 40082 coastal watersheds globally (see Supplementary Figure 1). We use coastal watersheds as a first order estimate of the spatial scale of coastal groundwater flow systems. While watersheds may be a good first order approximation of the contributing area for coastal groundwater systems, model experiments 56 and water budget analysis of river basins 57 indicate that there may be a significant component of regional groundwater flow that crosses watershed boundaries in some systems. To account for this we used a contributing area that is two times the size of coastal watersheds as an upper estimate for the calculation of coastal groundwater discharge, which covers the upper end of reported percentages of regional flow in previous work 56,57 . We use a modified version of the global permeability map 53,54 where higher values of permeability were assigned to the predominantly karstic coastal carbonates and unconsolidated sediments, which are frequently dominated by coarse deposits in coastal settings (see Supplementary Information section S2). Topographic gradients were calculated by dividing the elevation by the distance to the coastline. Two alternatives for the representative gradient were used that included the average gradient of all points in a coastal watershed and the gradient of points that coincided with a stream only.

Global coastal groundwater discharge fluxes
Global fluxes of fresh SGD and terrestrial discharge were obtained by linear interpolation of 357 model results to the 40082 coastal watersheds (see Supplementary Fig. S8 and S9). The interpolation is based on a comparison of permeability, topographic gradient and recharge input, which is the product of groundwater recharge and the size of the contributing area, for the model results and the geospatial data of coastal watersheds. We used reported ±1 standard deviation uncertainty of permeability 54 , the differences between two alternative recharge datasets 23,55 and the differences in the topographic gradients between elevation grid nodes covered by streams and the entire coastal watershed to calculate minimum and maximum estimates of the discharge fluxes. The interpolated discharge values were compared to published values of surface runoff 23 , potential evapotranspiration 22 and tide and storm surge elevation 58 , which were assigned to each watershed using GIS tools.

Model-data comparison
Comparison between modeled and measured hydraulic gradients in 336 coastal watersheds with water level observations from a global dataset 59 Table S7). However, for the remaining locations the reported values strongly exceed the modeled values, but the reported fresh SGD also exceeds the total groundwater recharge in adjacent coastal watersheds ( Supplementary Fig. S12b), which suggests that these studies strongly overestimate fresh SGD. The frequent inconsistency of fresh SGD estimates with onshore groundwater budgets has been noted by several earlier authors 11,60 and may be due to uncertainties in methods to quantify fresh water discharge 61,62 or biased selection and reporting of study sites for SGD 1 .

Quantifying solute transport, eutrophication risk and groundwater depletion
First-order estimates of the global solute flux that is transported by fresh SGD were calculated by contributed to the eutrophication of terrestrial and nearshore ecosystems 33 . We compared the modeled coastal discharge fluxes to groundwater depletion by using published groundwater depletion rates 38 and multiplying these rates with the area of coastal watersheds.

Acknowledgments:
EL and TG were supported by NSERC and CIFAR. NM was supported by the BMBF project SGD-NUT (Grant #01LN1307A). We would like to thank Lou A. Derry for a very helpful review of an early version of this manuscript.

Author contributions
TG conceived this study. EL and TG designed the methods. EL constructed the model code and conducted the geospatial analysis and numerical modeling. All authors contributed to interpretation and discussion of the results and writing the manuscript.

Data availability
The results of the model sensitivity analysis and the parameter space exploration are available as supplementary data files S1 and S2, respectively. The results of the geospatial analysis and the interpolated coastal groundwater discharge fluxes for the global watersheds are available in shapefile format as supplementary data file S3. The data files are also available on Pangea (url to be supplied).

Code Availability
The model code used to simulate coastal groundwater discharge (GroMPy-couple) is available at

Supplementary information
The supplementary information is divided in six sections: S1) Geospatial data analysis of coastal watersheds, S2) Numerical model of coastal groundwater discharge, S3) Calculating global distributed coastal groundwater discharge, S4) Assessment of solute transport and eutrophication risk, S5) Sensitivity analysis and S6) Model-data comparison. The first four sections provide more details on the methods to calculate coastal groundwater discharge and the last two provide more detailed information on the model results.

S1: Geospatial data analysis of coastal watersheds
We analyzed watershed length scale, topographic gradient, groundwater recharge and near-surface permeability for n=40082 watersheds globally. The results are shown in Fig. S1 and Table S1.

Watershed geometry
We used the geometry of coastal watersheds as a first order estimate of the size of the area that contributes to coastal groundwater discharge. We acknowledge that groundwater flow does not necessarily follow surface water divides. In many cases the majority of groundwater flow is local, which means it discharges in the nearest surface water feature and surface watersheds can be a good approximation of the size of contributing areas 1 . However, in areas with high permeability, high relief or low recharge the watertable can be decoupled from topography and regional flow that bypasses the nearest discharge points can be significant 2,3 . Comparison of recharge and discharge estimates in river basins indicates that for the majority of basins the regional flow component is less than 50% 4 . To cover the uncertainty of the contributing area in our calculations of coastal groundwater discharge, we used a value of two times the size of surface watersheds as a maximum estimate.
First, 40,082 coastal watersheds were selected using global watershed 5

Permeability
Permeability of each coastal watershed was extracted from a global dataset of near-surface permeability (up to ~100 m depth) 7 . The permeability map is based on a high-resolution global map of surface lithology 8 and a compilation of large-scale permeability estimates of near-surface geological units 9 . We made two changes compared to the permeability map with the aim of ensuring that the modeled coastal groundwater discharge is a high, but still conservative and realistic estimate.
In the global permeability map, areas in which the lithology consisted of mixed unconsolidated sediment or unconsolidated sediments with an unknown grain size were assigned a relatively low permeability (10 -13 m 2 ) 9 . This value is below the threshold for generating significant coastal groundwater discharge, except for settings with a very high topographic gradient. However the multimodal distribution of this unit 9 suggests that in many cases this unit contains coarse grained sediments. In layered unconsolidated sediments, the effective permeability is likely to be close to the value of the most permeable sub-unit, whereas the permeability value assigned in the global permeability map is the mean permeability on a log scale. We instead assume that the permeability for coarse-grained unconsolidated sediments (10 -10.9 m 2 ) is more appropriate to simulate regional groundwater flow in coastal aquifers consisting of mixed or unknown unconsolidated sediments.
The permeability of carbonates in the global permeability map is likely underestimated in coastal areas with strong karstification. Coastal carbonates are predominantly karstified 10 , in part because the fresh-salt water mixing zone at the coastline promotes dissolution and the formation of permeable karst conduits 11 . We therefore adopted a higher permeability estimate equal to the value reported in the global permeability map plus one standard deviation (10 -10.3 m 2 ). While we acknowledge that this value is highly uncertain as a global average for coastal karst aquifers, and that in reality their permeability is likely to be highly variable, this value is in line with reported values of regional-scale permeability in coastal karst aquifers [12][13][14] .

Groundwater recharge
Groundwater recharge was derived from a calibrated global groundwater model 15 . We also used an older global recharge model 16 , which generally predicts lower values of recharge and was used as a lower bound for the uncertainty analysis of global submarine groundwater discharge.

Topography
Elevation data 17 was extracted for each coastal watershed. We calculated distance to the coastline and elevation for each point in the elevation raster. We then calculated the average and standard deviation of the topographic gradient by dividing elevation by the distance to the coast. In addition, we calculated the topographic gradient for raster cells that contained a stream. The locations of streams were obtained by grouping raster cells for which the distance to the coast is the same within a range that is equal to the size of one raster cell, and then selected the raster cell with the lowest elevation for each distance value. In addition we extracted the average topographic gradient for the first 1000 m from the coastline.

Additional geospatial datasets
We extracted a number of global datasets for a comparison with the calculated submarine and terrestrial discharge, including river discharge 15 , potential evapotranspiration 18 , the elevation of tide and storm surges 19 and watertable gradient, which was quantified using a global model of watertable depth 20 and global elevation datasets 17 .

S2: Numerical model of groundwater flow in coastal aquifers
We modeled coupled density-dependent groundwater flow and solute transport in coastal aquifers using a new model code, GroMPy-couple, which is a Python shell around the generic finite element code escript 21,22 . Escript is used to solve the governing equations for coupled solute transport and density-driven groundwater flow. The choice for escript was motivated by the fact that the code runs on parallel computing facilities which strongly increases computational efficiency and by the Python interface of escript, which provides a tool for the automated execution and analysis of large numbers of model runs. An alternative like the more well-known model code SEAWAT 23  has been applied to model coastal groundwater flow.

Governing equations
Density dependent groundwater flow was modeled by solving the following equation 27 : Where is storativity (s 2 m -2 ), P is pressure (Pa), t is time (sec), C is solute concentration (kg kg -1 ), is fluid density (kg m -3 ), is porosity (dimensionless), is solute expansion coefficient (dimensionless), k is permeability (m 2 ), is dynamic viscosity (Pa s), g is the gravitational constant  Table 2 and were found using linear regression to published equations of state 29 over the range of salinities of 0.0 to 0.035 kg kg -1 .

Solving the flow and solute transport equations
The solute transport and groundwater flow equations were solved using the finite element code escript 21,22 . We follow an iterative coupling algorithm 30 in which the solute transport equation (equation S2), equations of state (equations S5 and S6) and groundwater flow equations (equations S1 and S4) are solved iteratively until the change in pressure and concentration between two iterations is less than 1 × 10 -4 Pa and less than 1 × 10 -7 kg kg -1 , respectively.

Model geometry
Groundwater flow was simulated in a two-dimensional cross section of the subsurface. While we acknowledge that coastal groundwater flow is a three dimensional process, the computational demands of running large numbers of three dimensional models would be prohibitive, given the high spatial resolution required to accurately model density-driven flow and the constraints on timestep size imposed by numerical stability of modeling advective solute transport 31 .
We assigned a constant linear slope to the terrestrial and marine parts of the model domain.
The linear slope is a simplification. Testing more complex topographies, with for instance a lower slope of the near-shore parts that is often found in sedimentary settings, or conversely high relief and cliffs in erosional settings would increase the number of model runs that are needed to cover parameter space significantly, and would make the computational costs prohibitive. We therefore aimed to cover the first order effects of a linear topography for the study presented here.
The first 1000 m of the model domain are covered by seawater. The length of the landward size of the model was constrained by the watershed length described in section S1. The thickness of the model domain was kept at 100 m for the exploration of the parameter space, and was varied between 100 and 500 m for sensitivity analysis.
We applied a spatial discretization that varied from 3 m in a zone that extends 500 m offshore and 250 m onshore to 10 m on the landward boundary of the model domain. The zone with fine discretization is centered around the fresh-salt water interface that was calculated using an analytical solution 32 : Where y is the depth of the interface below sea level (m), K is hydraulic conductivity (m s -1 ), which is calculated as = / , ρf and ρs are the density of freshwater and seawater, respectively (kg m -3 ), Q is the discharge rate of fresh groundwater (m 2 s -1 ) and x is distance to the coastline (m). The discharge term (Q) in eq. S7 was calculated using a depth-integrated version of Darcy's law: Where b is aquifer thickness (m), and h is hydraulic head (m). In case the calculated discharge exceeded the total recharge input (i.e., the product of the recharge and length of the model domain) into the system, discharge was capped to equal the recharge input.

Initial and boundary conditions
A specified recharge flux boundary and a seepage boundary condition were applied to the upper model boundary at the terrestrial side of the model domain. For the seaward side of the model domain we applied a specified pressure that equals the load of the overlying seawater. Initial salinity was equal to seawater values of 0.035 kg kg -1 under the seabed and in a saltwater toe that extends inland following equation S7. Initial pressures were calculated by solving the steady-state version of the groundwater flow equation (eq. S1).
The exchange of groundwater and surface water or evapotranspiration was simulated using a seepage boundary algorithm 33 . The seepage algorithm was chosen because it represents a more realistic upper boundary than often used fixed specified pressure or flux boundaries, while avoiding the computational cost of explicitly modeling evapotranspiration and surface-groundwater exchange.
The seepage boundary condition was implemented using an iterative procedure. First, initial pressures were calculated by solving the steady-steady-state groundwater flow equation (eq. S1, with the derivative of pressure and concentration over time set to zero). For the first step a specified flux was assigned to the entire top boundary, which represents groundwater recharge from precipitation. Following the first iteration step, a specified pressure boundary was adopted at any surface node where the fluid pressure (P) exceeds 0 Pa. Fluid pressures were then recalculated again by solving eq. S1 using this new boundary condition. Following each iteration step, the flux to the To avoid oscillations in the solution, only the nodes that generate 10% or more inflow compared to the seepage node with the highest rate of inflow are removed from the boundary condition after each iteration. This iterative procedure is continued until the number of seepage nodes reaches a steady value.
The iteratively calibrated steady-state seepage boundary is used as an initial seepage boundary during the transient model runs. At each timestep, the active seepage boundary is inherited from the previous time step. The seepage boundary condition is removed for any node that has become a net source of water into the model domain. Any non-seepage node at the surface where the fluid pressure exceeds 0 Pa is added to the seepage boundary. This implementation of the seepage boundary condition ensures that the hydraulic head never exceeds the surface elevation, and that there is not more inflow than the specified recharge rate at any node at the surface. Both possibilities would be unrealistic, but allowed when either a specified flux or a specified pressure boundary condition would be used for the upper boundary. In addition, the seepage boundary method as implemented here avoids the use of an unknown and uncertain drain conductance parameter that is used in drain boundary conditions, which aim to provide a similar realistic upper boundary as the seepage boundary 33 but use a different algorithm to achieve this.

Assumption of constant thickness and saturation
The model domain was assumed to be fully saturated and the saturated thickness was constant in the model domain and independent of the modeled pressures and hydraulic head. This is a simplification that avoids the numerical instability and high computational costs of modeling We used a constant thickness over the model domain. The thickness was varied between 50 and 500 m in a first sensitivity analysis. For the final model runs we adopted a standard thickness of 100 m, which is equal to the median aquifer thickness that were used to compile data for the global permeability map 7,9 , and is roughly equal to the thickness where the majority of young groundwater and active groundwater circulates following global compilations of radiogenic isotope data of groundwater 34,35 . Adopting a standard thickness of 100 m ensures that the modeled values of transmissivity (the product of permeability and thickness) are consistent with the global permeability map.

Model runtime
The transient models were run until a steady state was reached. We assumed that the model has reached steady state when the change in pressure is less than 1 Pa a -1 and the change in solute concentration that is less than 1 ×10 -4 kg kg -1 a -1 . The initial timestep size was 5 days and was increased by a factor of 1.03 after completing each timestep. To avoid numerical instability in solving advective solute transport, the maximum size of the timestep was adjusted automatically to not exceed the Courant-Friedrichs-Lewy (CFL) condition 36 : where q is fluid flux (m s -1 ), Δt is timestep size (s) and Δx is the size of one element (m) as calculated by escript 21 . We used a relatively low limiting value of CFL=1.0 to ensure numerical stability for the large set of different models that we tested.

Model validation
To assess how well the model can simulate fresh and salt water discharge in coastal aquifer systems we compared the model code with experimental data from a salt water intrusion experiment 37  The model simulates onshore discharge using a seepage algorithm 33 . The modeled extent of the area where seepage takes place was compared to an analytical solution of the seepage boundary 33 , which uses Dupuit assumptions, i.e., the vertical component of groundwater flow is assumed to be zero. The results are shown in Fig. S3. Our results agree show good agreement with the analytical solution, but consistently show a slightly smaller size of the discharge zone. This is in agreement with numerical experiments by Bresciani et al. 41 who shows that the analytical solution overestimates the size of the seepage zone due to the limitations imposed by the Dupuit assumption.
The modeled discharge fluxes over the land surface were compared to an analytical solution of groundwater discharge by Bokuniewicz 38 . The analytical solution assumes a constant linear hydraulic gradient over the land domain instead of the more realistic recharge and seepage boundary that we used in GroMPy-couple. We adopted a specified hydraulic head boundary condition for three model runs that were compared to the analytical solution. The modeled discharge flux matches the analytical solution perfectly (Fig. S4).
Finally, modeled fresh and salt water fluxes were also compared to modeled fluxes by Sutra.
A comparison of modeled terrestrial and submarine groundwater discharge simulated by GroMPycouple and Sutra is shown in Fig. S5. For these model setups the upper boundary condition was simplified to a fixed watertable (i.e., fluid pressure is zero) at the land surface. This is because Sutra does not have the option to simulate a seepage boundary. All other conditions were equal to the model runs that were used to explore submarine groundwater discharge in this manuscript. The results show a good match between Sutra and GroMPy-couple. GroMPy-couple predicts slightly higher recirculated submarine groundwater discharge fluxes, with an average of 12% of the total submarine groundwater discharge being saline in the Sutra model experiments and 15% in GroMPycouple. The cause for this is unknown. However, because the difference is very small compared to the other uncertainties (such as permeability) in modeling coastal groundwater discharge and both models show identical results in the salt-water intrusion benchmark discussed previously we consider GroMPy-couple sufficiently accurate to quantify coastal groundwater discharge.

Model sensitivity analysis and exploration of parameter space
First, we performed a sensitivity analysis of submarine groundwater discharge by varying watershed length, topographic gradient, groundwater recharge and permeability in a range covers the values found in the geospatial analysis (compare Fig. S1 and Table S1). In addition, we tested the effect of a realistic range of values of aquifer thickness and permeability anisotropy and longitudinal dispersivity for which no data was available. The sensitivity analysis consisted of a total of 130 model runs. The parameters ranges are shown in Table 3. The base case model used the median watershed length and groundwater recharge from the geospatial analysis, and a higher permeability and topographic gradient than the median watershed to better show the sensitivity of submarine groundwater discharge to the various parameters.
Second, we conducted a number of model experiments (n=495) to explore the parameter space for recharge input (recharge multiplied by contributing area), topographic gradient and permeability. The ranges of these parameters are shown in Table S4. Permeability anisotropy (the ratio of horizontal over vertical permeability) was kept constant at a value of 10. For fractured crystalline rocks vertical permeability may exceed horizontal permeability, whereas for layered sediment sequences anisotropy can reach a factor of 100 or more 44 . We used a constant anisotropy value of 10 to strike a balance between these two end-members.
Note that a more accurate implementation of permeability anisotropy in our models would require information on the orientations of fractures and bedding in coastal aquifers, which are currently not available at a global scale.

Representative topographic gradient
The cross-sectional models use a simple linear topographic gradient. The topographic gradient is important because it sets a maximum for the hydraulic gradient in each watershed and because it governs the partitioning between terrestrial groundwater recharge, terrestrial discharge and submarine groundwater discharge. We evaluated two metrics as representative topographic gradient: the average topographic gradient of coastal watersheds and the topographic gradient of drainage features (ie., streams) only. To explore which metric provides the best model of coastal groundwater discharge we compared a series of numerical model runs which used a 2D map-view model that included the full topography of coastal watersheds with a series of simplified crosssectional models. The map-view models were conducted using a standard finite difference model code of steady-state depth-integrated groundwater flow implemented in the programming language Python 45 and employed an identical seepage algorithm 33 as described previously for the crosssectional models to simulate groundwater discharge. The map-view models include the full topography of each raster cell in the elevation dataset, and are therefore expected to provide a better estimate of groundwater recharge and discharge processes than the cross-sectional models.
The map-view models simulate depth-integrated groundwater flow in a single layer, which is a common assumption in regional-scale groundwater models. The models used elevation data from a global digital elevation dataset 17 , and recharge and permeability following the results of the geospatial analysis described in section S1. We simulated groundwater flow and discharge in 59 randomly chosen watersheds and compared the modeled terrestrial, near-shore and submarine groundwater flux in the map-view model with the fluxes simulated in the cross-sectional models. For this series of model runs both the map-view and cross-sectional models simulated single-density groundwater flow and did not include solute transport.
The results of a single model run for an example watershed are shown in Figure S6. A comparison of the recharge and discharge fluxes as a function of the distance to the coastline for the map-view and the cross-sectional models are shown in Fig. S6f. In total 36 % of the groundwater that is recharged discharges within 500 m of the coastline for the map view model. Groundwater discharge is distributed along topographic lows that coincide with surface water drainage features in the watershed. When projected to the distance from the coastline groundwater recharge and discharge are distributed unevenly (Fig S6f). In the cross-sectional models recharge is focused in a small area upstream and there is a large area where recharge equals discharge and the net flux over the top boundary is zero (Fig S6f). The cross-sectional model was rerun two times, each time with a different topographic gradient corresponding to the average gradient of the entire watershed and the average gradient of the main drainage channel in the watershed. The modeled discharge for these two cases is 39 and 11% of the applied recharge, respectively. This shows that in this case the cross-sectional model with the average topographic gradient compares best to the modeled discharge in the map-view model.
Comparison of the modeled coastal discharge fluxes for cross-sectional and map-view models for 59 randomly chosen watersheds is shown in Fig. S7. For 44 out of the 59 watersheds the modeled coastal groundwater discharge as a percentage of the total recharge in the map-view model is covered by the three cross-sectional models or the difference is 10% or less. Overall, the average topographic gradient of the entire watershed and the gradient of the streams are both equally good predictor of coastal groundwater discharge, with a coefficient of determination (R 2 ) value of 0.51 and 0.50, respectively. While cross-sectional models clearly do not capture the full variation of groundwater recharge and discharge, overall they area a relatively good first order approximation of coastal groundwater discharge, provided that models using the average topographic gradient and the topographic gradient of streams only are both taken into account.  Table 7). The threshold for nitrogen application is an approximate lower limit for nitrogen application in areas in Europe and the USA where nitrogen concentration in the vadose zone and groundwater have increased strongly over the last century 52 . Given the very high sensitivity of coastal ecosystems and especially coral reefs to nitrogen 57 , this threshold value is likely to result in a relatively conservative estimate of eutrophication risk. The calculated locations of sites with high eutrophication risk was compared with the location of sensitive ecosystems such as estuaries 58 , salt marshes 59 and coral reefs 60 . We acknowledge that eutrophication of marine ecosystems is a complex function of nutrient input, transport, denitrification and mixing with seawater. The eutrophication risk reported here should be considered as a first order estimate that can guide follow up studies.

S5: Model sensitivity
The model runs that were used to calculate the global coastal groundwater fluxes use fixed values of aquifer thickness, permeability anisotropy and dispersivity. The degree to which this may affect model result is shown by results of sensitivity analysis in Fig. S10. Compared to the much more sensitive parameters topographic gradient and permeability, coastal discharge is relatively insensitive to dispersivity and permeability anisotropy. In contrast, aquifer thickness does affect the discharge fluxes significantly. The adoption of a constant thickness for the global model results may therefore introduce significant uncertainties to the global model. We estimate this uncertainty as roughly a factor two, based on global estimates of the depth of young groundwater and active groundwater flow 34,35 . Note that this source of uncertainty is still much lower than the order of magnitude uncertainty in the global permeability dataset that was used in our analysis.

Watertable gradients
We compared the modeled groundwater observed watertable in these watersheds shows a reasonable agreement (Fig. S11). The median ratio is 1.06, which means that our model may slightly overestimate watertable gradients and underestimate permeability and coastal groundwater discharge. The model explains slightly less than half the observed variance in watertable gradients. The coefficient of determination (R 2 ) is 0.57, which suggest that while overall the model matches the gradients well it does not fully capture local variability and should be considered as a first order estimate that is representative for relatively large spatial scales (i.e., watershed size, ~11 km). While a calibration of the coastal groundwater discharge model to match local watertable data would be beyond the scope of this study, the agreement between the mean modeled and observed watertables indicates that the model results presented here are relatively robust.

Comparison with local estimates of fresh SGD
The modeled values of coastal groundwater discharge and fresh SGD were compared to local fresh SGD estimates. While there are abundant estimates of total SGD in the literature studies that provide a robust quantitative estimate of the contributions of the fresh component of SGD are relatively scarce. We compiled fresh SGD estimates that were predominantly obtained using seepage meters, where the fresh component of SGD was estimated using a salinity balance of the water discharging in seepage meters [61][62][63][64][65] or direct sampling of discharge 66 . In addition, a number of estimates were based on a combination of seepage meters and radon or radium isotopes in seawater 67,68 , seawater salinity anomalies to estimate the fresh water contribution 69 , and pore water chemistry and models of solute flux in the seabed 70 . We did not include the relatively numerous estimates of fresh SGD that were calculated using Darcy's law, because these implicitly assume that all water flowing in coastal aquifers discharges directly in the ocean, do not take into account the interaction of fresh and salt water at the coastline and diverge strongly from analytical models 38,71 or numerical models of coastal groundwater discharge like the model presented here.
Comparison of modeled and reported coastal groundwater discharge values in Fig. S12a and Table S5 shows that the modeled discharge is in several cases much lower than the values reported in many locations. Note that the reported values of fresh SGD were compared to modeled coastal groundwater discharge and not modeled fresh SGD, because in contrast to our approach the reported fresh SGD values in the literature often include water discharging above the mean sea level and below the high tide line or lack information to discern the onshore and offshore components of coastal discharge. Five studies 61,65,66,68,70 report values that are in the uncertainty range presented by our model results. The remaining five studies [62][63][64]67,69 show values that are one or more orders of magnitude higher than the model predictions. However, comparison of modeled values of permeability and reported values by five of these studies (see Table S5) show that with the exception of one study these values are in the same order of magnitude, and the misfit with reported fresh SGD values is not likely to be due to a systematic underestimation of permeability in our models.
Reported values of fresh submarine groundwater discharge are often equal or higher than the total recharge input in the system (Fig. S12b). All the fresh SGD estimates are from locations that have mapped rivers in the hinterland, which normally channel substantial part of the overall recharge in coastal aquifer as river baseflow and which means that only a small part of recharged groundwater can contribute to submarine groundwater discharge. Note that the reported fresh SGD values are from humid settings where models of groundwater recharge 15,16 are relatively reliable 72,73 . The contradiction between reported submarine groundwater discharge rates and the amount of groundwater available in coastal watersheds to supply fresh SGD has been noted by several previous authors 74,75 . The reasons for the mismatch may be the difficulty and high uncertainty of separating the fresh and recirculated components of SGD in local studies 76,77 , and potential bias of sampling locations towards sites of relatively high and focused discharge that is visible at the shoreline 78 .
Furthermore, to our knowledge sites with low or no SGD tend to not get reported in the literature.

Comparison with a karstic carbonate aquifer with no surface runoff
Initial estimates of coastal discharge assigned a permeability of 10 -11.8 m 2 to carbonate units following the global permeability map. The resulting estimates of coastal discharge were lower than expected for coastal karstic aquifers. For instance in the Yucatan peninsula surface water features are absent, and the majority of groundwater recharged is expected to discharge near the coast, along karst conduits and coastal and submarine springs 79 . Initial model runs predicted a coastal discharge of only 20% of the groundwater recharge for the eastern part of Yucatan. The modified permeability values where a higher permeability was assigned to carbonate units resulted in near-shore discharge that was equal to 87% of the groundwater recharged in the coastal watershed. This value is in accordance with discharge estimates for this region 79 . However, the total flux may still be underestimated, comparison with stream locations on Google Earth suggests that while overall the watershed and stream database 5 that support our analysis is robust, in areas with extensive karst such as Yucatan the stream density may be overestimated, which means that the representative length scale and the total groundwater input may be underestimated in our models. However, areas with continuous karstic carbonate cover over scales that exceed the scales of coastal watersheds (~11 km) are rare globally 8 .

Qualitative comparison with locations of significant SGD use and reported ecosystems impacts by SGD
A comparison of modeled fresh SGD and coastal groundwater discharge with nine locations where significant use of fresh submarine groundwater has been reported 80 is shown in Table S6. In addition, a selection of seven locations that have reported strong impacts of fresh SGD on coastal ecosystems is shown in Table S7. The latter is admittedly a small sample of the large body of literature on the impacts of fresh SGD, a more extensive compilation is planned as a follow up study. In most locations that report use of fresh SGD the modeled fresh SGD exceeds 50 m 2 a -1 (Table S6), which is much higher than the global median discharge of 0.4 m 2 a -1 . In two locations (Bahrain and Quissico, Mozambique) the modeled fresh SGD is very low. However, in both cases the modeled total coastal groundwater discharge is much higher (50 and 233 m 2 a -1 , respectively) than the fresh SGD component. Since the partitioning between fresh SGD and near-shore terrestrial discharge is highly dependent on topographic gradient the underestimation is likely due to the models not representing coastal topography well enough by using a single linear topographic gradient. The coastal (onshore and offshore) groundwater discharge in locations where effects on ecosystems have been reported (Table S7) is in most cases much higher than the median value of 30 m 2 a -1 . The nitrogen application in the adjacent coastal watersheds are highly variable and in some cases relatively low. This suggests that either other sources of nutrients or pollutants affect the coastal ecosystems or that agricultural nitrogen application has been underestimated in the relatively low-resolution global dataset that was used 56 . In addition, in these local studies it may have been in some cases difficult to separate the effects of surface water and groundwater input. Nonetheless overall the model results successfully match the relatively high discharges required to explain locations with significant use of fresh SGD or where coastal groundwater discharge impact the solute and nutrient budgets of coastal ecosystems.       For 75% of the watersheds the cross-sectional discharge falls between these values or the difference between the modeled discharge is less than 10%.     groundwater recharge (b). Five out of ten reported fresh SGD estimates are much higher than our model results (a), but these estimates are also equal to or exceed the total freshwater input in adjacent coastal aquifers (b). Numbers denote the reference for each estimate.       Grid cell size (m) 3 1 to 10, interval = 1 Table S4. Parameters ranges used for the exploration of parameter space and quantifying global coastal groundwater discharge.