Geological and groundwater flow model of a submarine groundwater discharge site at Hanko (Finland), northern Baltic Sea

Three-dimensional geological and groundwater flow models of a submarine groundwater discharge (SGD) site at Hanko (Finland), in the northern Baltic Sea, have been developed to provide a geological framework and a tool for the estimation of SGD rates into the coastal sea. The dataset used consists of gravimetric, ground-penetrating radar and shallow seismic surveys, drill logs, groundwater level monitoring data, field observations, and a LiDAR digital elevation model. The geological model is constrained by the local geometry of late Pleistocene and Holocene deposits, including till, glacial coarse-grained and fine-grained sediments, post-glacial mud, and coarse-grained littoral and aeolian deposits. The coarse-grained aquifer sediments form a shallow shore platform that extends approximately 100–250 m offshore, where the unit slopes steeply seawards and becomes covered by glacial and post-glacial muds. Groundwater flow preferentially takes place in channel-fill outwash coarse-grained sediments and sand and gravel interbeds that provide conduits of higher hydraulic conductivity, and have led to the formation of pockmarks on the seafloor in areas of thin or absent mud cover. The groundwater flow model estimated the average SGD rate per square meter of the seafloor at 0.22 cm day−1 in autumn 2017. The average SGD rate increased to 0.28 cm day−1 as a response to an approximately 30% increase in recharge in spring 2020. Sensitivity analysis shows that recharge has a larger influence on SGD rate compared with aquifer hydraulic conductivity and the seafloor conductance. An increase in recharge in this region will cause more SGD into the Baltic Sea.


Introduction
Submarine groundwater discharge (SGD) is the flow of groundwater from the seafloor to the coastal sea (e.g. Burnett et al. 2003;Moore 2010). It has been observed as a source of nutrients, trace metals and other potentially harmful substances from land to coastal ecosystems (e.g. Moore 2010; Szymczycha et al. 2012Szymczycha et al. , 2016Luijendijk et al. 2020). Growing recognition of the importance of SGD on the environmental state of coastal sea areas has resulted in an increasing number of studies, aimed at estimation of the rates of SGD and associated fluxes by different kinds of methods (e.g. Burnett et al. 2003Burnett et al. , 2008Schlüter et al. 2004;Peterson et al. 2008;Gleeson et al. 2013;Tait et al. 2013;Schubert et al. 2014;Sadat-Noori et al. 2015;Szymczycha et al. 2012Szymczycha et al. , 2016Krall et al. 2017;Idczak et al. 2020). Among those methods, numerical groundwater flow models have been developed as a means to estimate and predict the SGD and associated fluxes under current and future conditions (e.g. Andersen et al. 2007;Langevin 2003;Kaleris et al. 2002). Three-dimensional (3D) geological models represent the geometry and stratigraphy of depositional sequences, and provide the geological framework for groundwater flow models. Successful estimates of the SGD and associated fluxes to the coastal sea by numerical flow modeling require good understanding of the aquifer 3D geometry, the distribution of porosity in the host sediments, and the groundwater flow pathways.
Shallow groundwater in northern Europe typically resides in the thin overburden of late Pleistocene and Holocene glacial and post-glacial sediments that immediately overlie the Precambrian crystalline bedrock. These groundwater areas are hydrologically complex with heterogeneous sedimentary cover including till, glacioaquatic gravel and sand, glaciolacustrine and post-glacial silt and clay, and in some areas, reworked littoral gravel, sand and clay (Saarnisto and Saarinen 2001;Donner 2010;Räsänen et al. 2009). Major glacigenic landforms with coarse-grained sediments are important groundwater reservoirs in the region, including the First Salpausselkä ice-marginal formation (SSI) that continues tens of kilometers along the Hanko Peninsula in the south coast of Finland, and to the Baltic Sea (Häkkinen 1982;Fyfe 1990;Kielosto et al. 1996). For many locations along the Hanko Peninsula, groundwater flow models imply fresh groundwater discharge from the SSI to the coastal Baltic Sea (Luoma and Okkonen 2014). Indeed, Virtasalo et al. (2019) recently discovered the first SGD site in Finland, which is connected to the SSI. The groundwater discharge is associated with pockmarks on the seafloor at seawater depths between 4 and 17 m, with estimated SGD rates from 222 Rn measurements of approximately 0.4-1.2 cm day −1 . Comparable and larger SDG rates have been documented for other sites in the Baltic Sea (e.g. Schlüter et al. 2004;Szymczycha et al. 2012Szymczycha et al. , 2016. However, the significance of SGD and associated fluxes on the overall eutrophic environmental state of the Baltic Sea is poorly understood. The objectives of this study were to construct a 3D geological model of the Hanko SGD site by integrating previously collected and new geophysical and geological data, and, based on this knowledge, to develop a numerical groundwater flow model in order to provide estimates of the SGD rate to the coastal sea. The 3D geological model is built under a geographic information system (GIS) platform and the 3D stratigraphic module in the Groundwater Modeling Software (GMS). The groundwater flow model is constructed by using MODFLOW 2005 code (Harbaugh et al. 2017). Due to the low density and low salinity of the northern Baltic Sea, the density-dependent model is not applied.

General setting
The study area is located in the Isolähde-Lappohja groundwater area in the municipality of Hanko, south Finland, at approximately 59°89″ N 23°21″ E, and covers 12 km 2 of the mainland and the Baltic Sea (Fig. 1). The shallow aquifer in the study area consists of porous gravels and sands of the SSI ice-marginal formation and is bounded by the Baltic Sea. The topography in the study area varies considerably between onshore and offshore. A steep slope forms at the boundaries between aeolian hills from the main land-14-28 m above sea level (asl)-and the shoreline, and between the edge of the shallow shore platform (1-4 m bsl) at 100-250 m distance from the shoreline and the sea floor at 14-17 m bsl. The study area belongs to the temperate coniferous-mixed forest climate zone with cold, wet winters. The mean annual temperature measured at the Tvärminne weather station, approximately 5 km south of the study area, during 1981-2010 was approximately +6°C (Pirinen et al. 2010). The lowest daily temperatures are generally recorded in January and February, and the highest during July and August. The mean annual precipitation was approximately of 634 mm .
The sea is essentially nontidal, but irregular water level fluctuations of as much as 2.1 m take place as a result of variations in wind and atmospheric pressure (Wolski et al. 2014). The annual mean sea surface salinity ranges between 4.5 and 6.5 PSU, and the annual mean sea surface temperature between 4 and 9°C during the period 1927-2011 (Merkouriadi and Leppäranta 2014). The low salinity results from the high riverine runoff from the large Baltic Sea catchment area, and from the long distance to the narrow connection to the North Sea through the Danish straits.

Geology and hydrogeology
The aquifer in the Hanko area is underlain by the basement of Precambrian crystalline igneous and metamorphic rocks, and covered by late Pleistocene and Holocene glacial and post-glacial deposits. The Precambrian bedrock, which mainly consists of granite, quartz diorite and granodiorites, forms a sharp contact with the sedimentary cover, with some rock exposed in the area (Kielosto et al. 1996). The studied aquifer is situated in the First Salpausselkä icemarginal formation (SSI; Fig. 1), which was formed in the course of the last Weichselian and Holocence deglaciation of the Fennoscandian continental ice sheet (Fyfe 1990(Fyfe , 1991Kielosto et al. 1996). The SSI was deposited over a period of 217 years, by ice-margin retreat, ca. 12,100 varve years ago (with reference to year 2000), during the Younger Dryas climatic event (Sauramo 1923;Saarnisto and Saarinen 2001). The ice-marginal formation consists of glacial till, gravel, sand and clay, together with postglacial littoral gravel, sand and clay. The primary ice-marginal formation was deposited as a narrow ridge of contiguous meltwater fans and local feeding eskers that were formed along the ice-margin grounding line in an ice-contact lake that was more than 100 m deep at Hanko (Fyfe 1990). When the ice sheet withdrew from the area, the SSI was successively covered by glaciolacustrine rhythmically alternating (varved) silt and clay, postglacial lacustrine poorly bedded clay, and brackish-water organic-rich mud (Virtasalo et al. 2007(Virtasalo et al. , 2014. The sea level has been regressive since deglaciation due to glacioisostatic land uplift (today 4 mm/year; Kakkuri 2012), except for a shortlived transgression at the beginning of the brackish-water phase ca. 7,600 years ago (Virtasalo et al. 2007). The icemarginal formation was exposed to waves and eventually to wind as it gradually rose from the sea. The original cover of glaciolacustrine and postglacial lacustrine silts and clays was removed, and ridge morphology was truncated and flattened from the top (Fyfe 1990;Kielosto et al. 1996). The topographic landform varies between 14 and 28 m asl along the northern ridge of the SSI and its elevation decreases to less than 2 m asl along the northern and southern coastlines.
The shallow aquifer in the study area is unconfined with the thickness of the Pleistocene and Holocene deposits varying from less than 1 to 68 m, with the average thickness of approximately 25 m. The water table in the Isolähde-Lappohja inland area is at 9.6-13.2 m asl, and falls close to the sea level at the shoreline, where groundwater discharges to the Baltic Sea. The results of well testing and sediment sample analysis show that the hydraulic conductivity of the aquifer varies from 0.3 to 4.8 m day −1 in the silty sand and fine sand, to up to 27.7 m day −1 in the sand and gravel (Luoma and Okkonen 2014;Luoma et al. 2017). Groundwater recharge occurs mainly twice a year during spring (late March-early April) and late autumn (November-early December) from the infiltration of snowmelt and rainfall, respectively (Luoma and Okkonen 2014). Groundwater flows mainly towards the Baltic Sea in the northnorthwest and in the south-southeast.

Materials
The geological and groundwater flow models of the shallow groundwater in the study area were constructed by using all available geological, geophysical and hydrogeological data that were produced during this study and from previous studies. Locations of monitoring points that yielded data used for the interpretation are presented in Fig. 2. Data produced during this study consist of the following-12 km of groundpenetrating radar (GPR) profiles acquired in October 2017; 12 km of shallow reflection seismic profiles from the offshore area acquired in May 2017; 0.9 km of refraction seismic survey along the shoreline acquired in April 2018; and 3 km of a walking survey of electrical conductivity (EC) and temperature measurements of near-bottom water in the shallow area along the shoreline carried out in June 2017. The GPR profiles were collected using a GSSI SIR 4000 control unit with a GSSI antenna operating at 200 MHz central frequency. Data were recorded using a 220 ns time window. Location was recorded using a Trimble GeoXH 6000 handheld GPS receiver with VRS network correction. Location, topography, and GPR profiles were combined and processed using Geodoctor 3.2 software. Signal processing methods applied were background removal, low-pass and high-pass filtering, and linear gain. A relative permittivity value of 6 was used, which corresponds to sandy dry soil (Annan 2009). The GPR profiles were interpreted following Neal (2004) and were calibrated with the groundwater levels from four observation boreholes measured at the same period, and the drill logs of 12 groundwater observation wells. Details of acquisition, processing and interpretation of GPR and offshore seismic survey data are presented in Virtasalo et al. (2019). Data from previous studies consist of the gravimetric survey during 2004 obtained from the Geological Survey of Finland (GTK 2020) and courtesy of the city of Hanko (Sito 2018); information on groundwater observation boreholes (e.g. drill logs, groundwater levels) were obtained from the SYKE-POVET database (SYKE 2020); a LiDAR-based digital elevation model (LiDAR DEM) was obtained from the National Land Survey of Finland (NSL 2020). Groundwater level and temperature were monitored by autonomous measurements every 30 min at the observation borehole HP101 by pressure transducer from November 2018 to April 2020. Weather data (rainfall, snow thickness and surface temperature) during the same period were obtained from the FMI Tvärminne weather station (FMI 2020), about 5 km south of the study area ( Fig. 1).

Geological and conceptual model
The 3D geological model of the late Pleistocene and Holocene deposits of the shallow groundwater area in Isolähde-Lappohja was constructed under a GIS platform by using ArcGIS/ArcMap for two-dimensional (2D) spatial analysis and the stratigraphic module in the Groundwater Modeling Software (GMS) for 3D interpretation and visualization. In order to provide the geological framework for the groundwater flow model, geological and conceptual models of the aquifer area are described in this section. The late Pleistocene and Holocene deposits sharply overlie the crystalline Precambrian bedrock, which is clearly visible from the gravimetric survey. A grid map of the bedrock topography was produced from the interpolation of the bedrock surface elevation data obtained from the various data sources, e.g. bedrock outcrops, drilled bedrock, gravimetric, seismic and GPR data using the kriging and the inverse distance weighting (IDW) interpolation methods in ArcMap. In this study, the bedrock surface topography was identified throughout the study area as shown in Fig. 3. The bedrock surface topography varies between −68.3 and + 27.4 m asl with the average elevation of −14.5 m asl. High elevations of the bedrock topography are located in the west, the east and the northern coastline, where the bedrock surface is partly exposed to the surface. Low elevations of the bedrock topography are seen in two main depression areas: in the NW-SE direction across the SSI, and in the NE-SW direction parallel to the SSI ridge, where the deepest part is found in the Baltic Sea adjacent to the SSI formation. The thickness of the Pleistocene and Holocene deposits as obtained from the subtraction between LiDAR-DEM and bedrock surface topography is presented in Fig. 4, with the variable thickness between less than 1 and 68 m. The thickest sedimentary cover is found in the bedrock depression areas. A 3D geological model of the sediment units was constructed based on the drill logs of 12 groundwater observation wells that reached the bedrock, GPR profiles and the background geological setting of the area. GPR profiles provide useful information of the internal structure of the sediment units; however, the penetration depth of the GPR was limited to the uppermost approximately 20 m, which was insufficient to cover the full sedimentary cover. In addition, the LiDAR-DEM landform was also used to support the interpretation of the sediments at the ground surface.
Typical GPR profiles of the depositional units in this area are shown in Fig. 5, which presents the GPR profile f6 at the location between the profile f1 and the shoreline (see Fig. 2 for location). The GPR profile f6 illustrates steeply inclined planar cross-bedded sands and gravels, dipping to the SE direction. This unit represents the primary deposit of the SSI which progressed into the distal area in the SE. The upper parts of the inclined beds were truncated and overlain by a subhorizontal layer, which consists mainly of gravel. This erosional surface can be observed in the GPR profiles throughout the study area at almost the same level at approximately 17 m asl. The truncating layer was observed in nearby gravel pits to be clastsupported, rounded and well-sorted gravel, indicating the winnowing of sediments by high energy wave action during emergence. It was interpreted as a lag deposit of the primary SSI deposit from the littoral process in the beach ridge areas. The gravel layer was overlain by a aeolian fine sand deposit, which is clearly discernible in the LiDAR-DEM (Fig. 2). The LiDAR DEM in the Isolähde area reveals a large area of the ancient shorelines (beach ridges) which feature as a parallel ridge landform of coarse-grained sediments (gravels and sands) dipping toward the north-west, the proximal direction ( Fig. 2). Some parts of the beach ridges were deformed by aeolian processes and redeposited as dunes. Figure 6 presents a GPR profile f9 illustrating a channel-like feature with the thickness of 6-7 m next to the wp46. It was interpreted as an outwash channel-fill of coarse-grained sediments that are incised in finer-grained sediments and fully saturated with fresh groundwater. The strong reflectors at wp49 and wp52 were interpreted as the saturated coarse-grained deposits. The same features are found in the GPR profiles along the shoreline between the pockmarks E and B.
The late Pleistocene and Holocene depositional succession from bottom (the oldest) to top (the youngest) is composed of the following units: (1) till; (2) primary SSI deposit of coarsegrained sediments (gravels and sands); (3) glaciolacustrine fine-grained sediments (rhythmically alternating of silt and clay couplet layers); (4) post-glacial reworked coarsegrained littoral and aeolian deposits; and (5) post-glacial lacustrine silty clay and brackish-water organic-rich mud from the late depositional phases of the study area. The littoral deposit consists mainly of gravels and coarse-sands, while the aeolian deposit consists of fine sands. However, due to the sparse drill log data and uncertainties of top depths information, the littoral and the aeolian deposits were combined into a single post-glacial coarse-grained unit. Figure 7 presents the 3D visualizations of those depositional units. The hydraulic properties and the bulk volumes of each unit in the model domain are presented in Table 1.
The aquifer is composed of primary sands and gravels of the SSI and littoral reworked deposits that extend approximately 600 m in to the Baltic Sea. However, the greatest    Fig. 2 for location). Numbers 40-42 represent wp locations thickness of the sediments was found in the areas 100-250 m seaward from the shoreline. Further offshore, the aquifer was covered by thick glaciolacustrine rhythmite (silt and clay unit) and post-glacial lacustrine and brackish-water muds (maximum thickness of 43 m), which pinch-out on the shore platform slope. Based on the 3D geological model, the hydrogeological units of the groundwater area can be delineated as shown in Fig. 8. Boundaries of the groundwater areas are clearly identified based on the bedrock topography, the existence of fine-grained sediments and groundwater divide. Elevated bedrock topography separates groundwater areas in the west and partly in the east between Isolähde and Lappohja. The thick fine-grained layer in the southwest separates groundwater areas between Isolähde and Koverhar due to the poor hydraulic connectivity. A groundwater divide is located in the middle of the study area around HP100 and 2 . 5 9 × 1 0 7 7.89×10 7 6.09×10 7 8.74×10 7 2.65×10 7 Effective porosity (%) 10 20 1 25 1 Pore volume (m 3 ) 2 . 5 9 × 1 0 6 1.58×10 7 6.09×10 5 2.19×10 7 2.65×10 5 Explanation of symbols is shown in Fig. 2 HP106, where the groundwater flows into two directions; to the north-northwest and to the south-southeast. Thick finegrained layers in the Baltic Sea area overlie and constrain the extent of the aquifer, especially in the deeper part where the aquifer is partly confined or semiconfined. The very poor hydraulic conductivity of the thick fine-grained layers prevent SGD as well as the intrusion of seawater into the aquifer. Due to the complex depositional origin, the hydraulic conductivity (K) values of the primary deposits are expected to be quite low. High K values could be found locally in the channel-fill outwash coarse-grained sediments. Better estimates of K could be obtained from calibration of the groundwater flow model.

Model discretization
The groundwater flow model was constructed by using MODFLOW 2005 code (Harbaugh et al. 2017) under the GMS graphic environment. It is a two-layer model, covering an area of 2.9 km 2 . The model domain is discretized into 5 m × 5 m grid size and consists of a rectangular grid of 420 rows and 500 columns. Of the total model domain area, 1.1 km 2 (37.9%) is in the offshore area, where 0.76 km 2 (26.2%) of the total is overlain by post-glacial fine-grained (silt and clay) layer and 0.34 km 2 (11.7%) of the total is exposed on the seafloor and in direct contact with seawater.
Based on the conceptual model and hydrogeological information, individual layers are described as follows: Layer 1 consists of gravels and sands of the upper part of the Pleistocene and Holocene deposits, and covers sediments in the entire unsaturated and part of the saturated zones. The top of layer 1 is defined by the LiDAR DEM for the onshore area, and by seismic profiles for the offshore area. MODFLOW is a finite-difference groundwater flow code for the saturated flow zone. To avoid dry cells that might occur in the model domain during simulation, an additional thickness of 5 m was added to the groundwater level and was assigned as the bottom of layer 1. In addition, a fine-grained (silt and clay) layer was included in layer 1, and simulated as a seafloor layer in the MODFLOW River (RIV) package. The thickness of the model varied from 1 to 47.7 m (average of 16.6 m).
Layer 2 is placed between the bottom of Layer 1 and the bedrock surface topography. There is insufficient data to identify the distribution of till throughout the aquifer area. Thus, till was combined in layer 2, where the bedrock topography was well identified. The crystalline Precambrian bedrock has a very low hydraulic conductivity and was regarded as an impermeable layer forming the bottom boundary of the model. The thickness of the model varied from 1 to 50 m (average of 17.3 m).

Boundary conditions
The MODFLOW processing packages applied in this study include the General Head Boundary (GHB), no flow, River (RIV) and Recharge packages (Fig. 8). No water intake well was available during this study. The model boundary conditions were assigned based on the geological and hydrogeological conditions. The GHB was assigned for most parts of the boundary areas including the groundwater divide, the bedrock elevation areas in the west, north and north-east, and the fine-grained unit in the offshore boundary area. GHB allows groundwater to flow either into or out of the model domain depending on groundwater elevation changes along the boundary. Those areas can basically be assigned as noflow boundaries; however, to minimize effect of dry cells and nonconvergence with MODFLOW, the model boundaries in those areas were assigned as a GHB with very low conductance to ensure that the impact of in-and outflow along those boundaries are insignificant, and the simulation still performed smoothly without convergence problems. The areas inside the model domain that were truncated by the bedrock surface, were automatically assigned as a no-flow boundary. In the offshore area, the fine-grained layer was assigned as a seafloor layer for the RIV package. The conductance of the seafloor was assigned based on the hydraulic properties of the fine-grained materials. In the pockmark areas, the estimated discharge rates from Virtasalo et al. (2019) were used as the control points of the simulation. In the area that has no finegrained layer cover, the top of the submerged aquifer was assigned as the GHB with the variations of the sediment conductance depending on the sediments at the seafloor. The head stages of GHB and RIV were assigned at the average sea level (at zero meters).

Model input parameters
Hydraulic parameters, including horizontal saturated hydraulic conductivity (Kh), anisotropy-a ratio of horizontal (Kh) and vertical hydraulic conductivity (Kv)-specific yield and porosity, were initially assigned to the model layer based on hydrogeological data of the study area. The initial Kh values for the model input were obtained from the soil analysis data and the slug test analysis performed in the Hanko area (Luoma and Okkonen 2014). The estimated K values corresponded to the stratigraphy, varying from 0.3 to 4.8 m day −1 in silty sand and fine sand, and 5 to 30 m day −1 in sand and gravel. The spatial distributions of Kh values were assigned corresponding to the majority of the sediments in the model domain. The Kh values were finally adjusted during the model calibration process and were within the range of K values in the aquifer test data. No Kv values are available-instead, the anisotropy values were assigned as 10 and 5 for the fine-grained and coarse-grained areas, respectively. Likewise, the specific yield values of 0.20 and 0.22 and the effective porosity values of 0.1 and 0.25 were assigned for the fine-grained and coarsegrained areas, respectively, throughout the model domain area.
Groundwater recharge estimation depends on many factors including precipitation, soil type, vegetation, evapotranspiration, topography, depth to the water table, and land use pattern. The study area is part of the Natura 2000 nature conservation area, with Scots pine (Pinus sylvestris) forest. The aquifer is unconfined and groundwater recharges directly from the infiltration of snowmelt and rainfall. In this study, groundwater recharge was assigned as a net recharge, which is the total amount of infiltration that arrived at the water table. It was estimated from the water balance method based on the amount of precipitation and evapotranspiration, and the infiltration coefficients of different soil types: where P is precipitation (rainfalls and snowmelts), and PET is potential evapotranspiration. The weather parameters, including daily surface temperature and precipitation, were measured at the Tvärminne weather station, approximately 5 km south of the study area (FMI 2020). The daily potential evapotranspiration (PET, mm day -1 ) was estimated using a temperature-based method (Hamon 1963) as follows: where D is day length (hour), T is mean daily temperature (°C), and 'ea' is saturation vapour pressure (kPa) at mean daily temperature; ea ¼ 0:6108 exp 17:27 T = T þ 237:3 ð Þ ½ ð 3Þ The initial infiltration coefficients were assigned based on soil types, e.g. 0.4 was applied for the sand and gravel areas, and 0.1 for the fine-grained area. In the steady-state simulation, under the assigned net precipitation (P -PET), the infiltration coefficients were adjusted by association with the K values of the aquifer materials and the conductance of the seafloor to obtain the final recharge values. In the transient simulation, the K values of the aquifer materials and the conductance of the seafloor were fixed as the same condition as in the steady-state. The net precipitation (P -PET) was calculated based on the weather data, and the infiltration coefficients were adjusted during the calibration to obtain the final recharge values.

Model calibration
Due to the sparse groundwater observation wells in the model domain, groundwater levels from GPR profiles were also used for the calibration. Groundwater levels from GPR profiles were calibrated with water levels from the observation wells prior to calibration with the simulated groundwater level. Altogether four groundwater levels from the observation wells and 93 data points from the GPR profiles measured in autumn 2017 were used for the calibration of the steady-state flow simulation. In addition, the estimated discharge rates, derived from the 222 Rn measurements of the SGD pockmarks from Virtasalo et al. (2019), were used for the calibration. The calibration process was carried out by using the automatic parameter estimation in PEST (Doherty 2010) and trial-anderror by manually adjusting the K values, groundwater recharge, and conductance of the GHB and river bed (or seafloor in this study) boundaries until the best fit was obtained between the observed and simulated data. The transient flow model was constructed by using a daily stress period and time step of 406 days from February 2019 to April 2020. The calibrated groundwater level from the steady-state flow was used as an initial head for the transient flow. In addition, the calibrated K values of the aquifer materials and the conductance of the seafloor from the steady-state flow were used for the transient flow model with variations of the recharge estimated from the precipitation data during February 2019 to April 2020. The simulated transient flow was calibrated with the groundwater level monitoring data at the same period. As mentioned earlier, the net precipitation (P -PET) values were calculated based on the weather data and the infiltration coefficients were adjusted during the calibration to obtain the final recharge values, so that the simulated groundwater level agreed with the observation data. However, the monitoring data were available only from HP101 observation well. The changes (in percentages) of infiltration coefficients from the steady-state condition in areas around HP101 were also applied to the infiltration coefficients in the rest of the model domain. The calibrated flow field obtained from each simulation was used to calculate the particle tracking groundwater path lines by using the MODPATH program (Pollock 1989).

Sensitivity analysis
Sensitivity analysis was performed to examine the contribution of single parameters to the SGD rate by arbitrarily changing the values of a parameter by certain percentages (10-50%) from the initial condition, while the other parameters remained unchanged, and vice versa. In this study, groundwater recharge, Kh value and the seabed conductance of the GHB and RIV models were selected for the sensitivity analysis for the whole of the model domain. In addition, a local sensitivity analysis was conducted for the offshore area, in order to investigate the effect of changes in the GHB and RIV conductances on the SGD rate between high and low permeability regions of the seafloor. The SGD rates were calculated based on the results of the simulations and compared with the SGD rate for the autumn 2017 condition.

Field investigation and groundwater level monitoring
Based on data from the observation boreholes and the GPR profiles, groundwater discharged to the sea at elevations of approximately 0.0-0.5 m bsl. The EC of groundwater from the observation wells in the inland area varies between 6 and 10 ms m −1 (freshwater values), and along the shoreline varies between 17 and 30 ms m −1 (freshwater values), while the EC of seawater is 1,030-1,080 ms m −1 (brackish water values). Temperature of groundwater is quite constant all the time at approximately 7.0-7.3°C during the 1.5 years monitoring at HP101, while temperature of seawater varies during the year and was approximately 16-19°C during the measurements in summer 2017. No obvious groundwater influence was detected during the walking survey of the EC and temperature measurements in the shallow water area along the shoreline, especially in the potential areas where the seafloor consists of cobbles, gravels and coarse sands. The SGD probably was washed away and mixed with seawater by waves and currents immediately above the seafloor.
Groundwater level data from the observation well HP101 and weather data (surface temperature and precipitation) measured during November 2018-April 2020 are presented in Fig. 9. The total annual precipitation measured at the Tvärminne weather station was 404, 669 and 410 mm during the years 2018, 2019 and 2020 (as of July 2020), respectively. The mean annual precipitation in this area during 1981-2010 was approximately of 634 mm (Pirinen et al. 2010).
Groundwater levels varied between 11.35 and 13.12 m asl and corresponded to the change of weather data. In spring 2020, groundwater level increased by 1.70 m from 11.42 m asl on 20 October 2019 to 13.12 m asl on 2 April 2020, due to increased precipitation and reduced evapotranspiration and snow cover during that period.  Table 2 presents water balance changes between outflow and inflow, which were less than 0.01% for both simulations. The total amounts of the SGD fluxes were the sum of outflow from the seabed leakage (RIV package) and the general head boundary (GHB). Groundwater level and recharge in autumn 2019 are similar to autumn 2017, which share the common low groundwater levels during dry season in this area. Figure 11a,b present the spatial distributions of the calibrated horizontal hydraulic conductivity values (Kh values) of model layer 1 and layer 2, respectively. Figure 11c,d present the estimated recharge from the calibrations with the observation data in autumn 2017 and spring 2020, respectively. In autumn 2017, the total groundwater recharge to the aquifer Fig. 9 Groundwater level data from the automatic monitoring in observation borehole HP101 between autumn 2018 and spring 2020. Precipitation and surface temperature data from the same period are from the Tvärminne weather station, 5 km south of the study area (FMI 2020) is approximately 758.1 m 3 day −1 ( Table 2). The total amount of SDG to sea is approximately 773.2 m 3 day −1 , with an average of 0.22 cm day −1 per square meter of seafloor during the same period. The calibrated K values of the aquifer materials vary between 0.01 m day −1 in the finer-grained sediment area and 35 m day −1 in the outwash coarse-grained gravel and sand areas. In spring 2020, the estimated recharge increased by approximately 30% compared to autumn 2017 and autumn 2019. The amount of SGD also increased to the average of 0.28 cm day −1 .

Groundwater flow modeling
Figure 11e,f present the calibrated groundwater levels in autumn 2017 and spring 2020, respectively. Figure 11g,h present the areas of SGD and the SGD flux rates estimated from MODFLOW simulations during autumn 2017 and spring 2020, respectively. Figure 12 presents a 3D visualization of the model grids illustrating groundwater flow paths of the autumn 2017 simulation. The flow is directed from the groundwater divide in the north-west toward the discharge areas in the Baltic Sea in the south-east. Groundwater discharges directly to the sea in the shallow water area near the shoreline, whereas in the deeper part of the shore platform and on the slope of the platform, groundwater is confined by the fine-grained unit, resulting in steep upward flow paths with focused discharge to the sea at locations of high hydraulic conductivity such as pockmarks. Table 3

Model limitations
The model was developed with sparse data; thus, some limitations and uncertainties remain as described in the following: 1. Owing to the relatively low salinity (ranges between 4.5 and 6.5 PSU) and low density (average 1.005 kg m −3 ) of the Baltic Sea, compared with that of oceanic water (salinity of 35 PSU and density of 1.025 kg m −3 ), the contribution of seawater to the aquifer in the Hanko area was reported to be less than 0.3% and has no significant impact on the nonpumping coastal aquifer (Luoma and Okkonen 2014). Therefore, density-dependent groundwater flow and transport were not included in the MODFLOW model. 2. Groundwater recharge was estimated using the simplified water balance method, which is sufficient for the modeling and was finalized by calibrating with the observation data. The vadose zone is quite thick in some areas, and estimation of the recharge can be improved with more detailed site investigation. GWL groundwater level, GHB general head boundary Fig. 10 Comparisons between simulated and observed groundwater levels for the a steady-state model; and b transient flow model 3. Due to the sparse groundwater observation boreholes in the model domain, groundwater levels from GPR profiles were also used for the calibration. Special care must be taken when deriving groundwater levels from GPR profiles compared to observation boreholes, as the GPR reflection is impacted by the capillary transition zone (Igel et al. 2013). Moreover, the transient simulation was calibrated with the groundwater level data from one observation well. The data applied in the calibration process based only on this monitoring data may not be representative for the whole of the model domain. 4. The SGD mainly took place in the shallow water areas where the measured discharge rates were not available during the field investigation. The SGD rate was estimated and was entirely depended on the simulation results. A validation of the results could be done by detailed site investigation at a later stage.

Discussion
Coarse-grained outwash channel fills are important seaward groundwater conduits in the glacigenic SSI ridge. A NW-SE oriented bedrock depression, which is the interpreted direction of the glacial meltwater paleocurrent, could represent the main route of sediment transport to the ice contact lake, as high proportions of coarse-grained sediments are found in this area. In many places along this depression, the coarse-grained sediments were deposited as high-porosity, channel-fill outwash sediments perpendicular to the SSI ridge morphology (Fyfe 1990). On the other hand, depressions of the bedrock topography in the NE-SW direction, parallel the SSI ridge, contain thick glacial and post-glacial fine-grained sediments that  confine the aquifer. Fyfe (1991) described a channel-fill outwash of highly porous gravels and sands with thickness of approximately 7 m that cut across the SSI in a gravel pit in Tammisaari, approximately 20 km north-east of the study area. These channel-fills could be the main pathways for the groundwater flow, especially in the shallow shore platform where the fine-grained layers are thin or absent. Virtasalo et al. (2019) documented a number of pockmarks on the shore platform slope, where sandy sediments are exposed on the seafloor. The simulation of the groundwater flow model was performed with sparse data and many constraints. The recharge was associated with the K values during the calibration. Based on the assigned recharge, the calibrated K values of the aquifer materials are quite low, generally less than 5 m day −1 . However, the low K values could be caused by the origins of the sediments, especially in the deeper part of the aquifer (layer 2), where the sediments are a mixture of coarse and fine-grained sediments of the primary deposits with increasing proportions of fine-grained sediments in the distal part. Virtasalo et al. (2019) concluded that the aquifer sediments, at least in the offshore area, are mainly composed of fine sands, and represent the distal part of a subaqueous icecontact fan. Interbeds of coarser sand and gravel provide conduits of higher hydraulic conductivity for the groundwater flow, and lead to the focusing of the flow to the pockmark sites.
The simple groundwater discharge rate was estimated based on the Darcy equation between two known groundwater level locations as follows: where V is the seepage velocity of the groundwater flow [L T −1 ], K is hydraulic conductivity [L T −1 ], i is the hydraulic gradient [L L −1 ], and n is the porosity [L L −1 ]. Simple calculations of the SGD rates at different parameters and locations are presented in Table 4. Head 1 is groundwater level at the groundwater divide during autumn 2017, and head 2 is groundwater level at the shoreline. The results show that the SGD rates vary depending on the hydraulic gradient, K values and porosity values. The groundwater flow simulations estimated the SGD to the Baltic Sea at the average rates of 0.22 cm day −1 (range 0.0-1.21 cm day −1 ) and 0.28 cm day −1 (range 0.0-1.60 cm day −1 ) per square meter of the seafloor for autumn 2017 and spring 2020, respectively. Virtasalo et al. (2019) estimated the SGD rates based on the 222 Rn measurements of groundwater samples from borehole HP101, and pockmarks E, D and B during autumn 2017. The SGD rates varied between 0.4 and 1.2 cm day −1 between pockmarks E and B, respectively. These rates are consistent with the results of this study at the K values of 0.1 m day −1 (Table 4), which is expected for fine sands in the pockmarks. However, higher SGD rates may take place locally on the shore platform and slope areas with higher K values. Virtasalo et al. (2019) observed more active pockmarks in the eastern part of the shore platform. Sediment samples collected from pockmarks D and B in the east are composed of fine sand, whereas in pockmark E the fine sand is covered with soft organic-rich mud. The organic-rich mud layer could continue to the west, where the thick layer of fine-grained sediments is found in the area around SW3-Ova8 (Fig. 2). This thick fine-grained layer could prevent SGD in the western part of the shore platform. A similar setting was described by Andersen et al. (2007) where the SGD was largely controlled by the geological structure of the aquifer and the fresh groundwater discharge predominantly occurred within a narrow zone of the upper 10-15 m of the intertidal zone along the shoreline.
The total submarine groundwater discharge was estimated at approximately 773.2 m 3 day −1 in autumn 2017 and at 986.7 m 3 day −1 in spring 2020 ( Table 2). The values are modest compared to the estimated annual average direct groundwater discharge from the whole Finnish coastline to the Baltic Sea, which is approximately 1,000,000 m 3 day −1 (Peltonen 2002). Direct groundwater discharge from comparable glacial meltwaterdeposited landforms in Finland to the sea has been estimated at 173,000 m 3 day −1 (2 m 3 s −1 ; Mälkki 2003). For comparison, the mean annual discharge of the nearby Karjaanjoki River is 1,650,000 m 3 day −1 (Saura et al. 2010).
The simulated flux rates of SGD for autumn 2017 and spring 2020 showed high SGD flux rates in the shallow shore platform along the shoreline, where the coarsegrained sediments are exposed directly to the sea. These agree with the GPR profiles and drill log data where the reflections of the coarse-grained materials were observed in the shallow part of the aquifer. Also, the wave reworking of the shore platform surface may have removed finegrained sediments and increased the hydraulic conductivity of the aquifer surface materials, thereby enhancing groundwater discharge to the sea. However, groundwater influence was not detected during the field investigation of the shallow water area. The actual SGD rate could possibly be very low and the walking survey measuring the EC and temperature in the shallow sea area may not have been suitable for the measurement under high wave conditions. On the other hand, this may imply seasonal variation of the SGD. The simulation indicated strong relationship between the SGD and recharge. Based on the monitoring data, the highest recharge occurs during spring (April to early May). Virtasalo et al. (2019) observed the melted spots on the sea ice along the shoreline during spring, which could be caused by the discharge of warmer groundwater (average temperature of 7°C). A similar situation was described by Szymkiewicz et al. (2020), where the rate of SGD varies seasonally and in relation to recharge, with peaks in late winter to early spring. In this study, the measurement of EC and temperature at the sea bottom along the shoreline was carried out in summer when the groundwater level and recharge are low. As a consequence, the SGD rate to the seafloor could have been low during that time. Based on the EC data, the submarine fresh groundwater portion of the total SGD flux can be estimated by using the following equation: where Qf is simulated fresh groundwater discharge per unit area of seafloor per unit time [L 3 L −2 T −1 ], and Qt is simulated total groundwater discharge [L 3 L −2 T −1 ]. ECm, ECf and ECs are the EC values of the measured shore-platform water, fresh groundwater and the Baltic Sea water, respectively. The measured EC values on the shore platform varied between 1,030 and 1,080 ms m −1 . Given the EC of the Baltic Sea water of 1,080 ms m −1 , the EC of fresh groundwater from observation well HP101 (inland) of 10 ms m −1 , and the lowest measured value of 1,030 ms m −1 on the shore platform, the Qf was estimated at 0.0467 Qt, which is a very small fraction of freshwater discharge to the seafloor during summer.
Groundwater level monitoring data are available from only one observation well (HP101) and may not represent changes in the groundwater level throughout the study area. This caused high uncertainty in the transient simulation. By applying the same configuration of the aquifer materials, e.g. the same K values and seafloor conductance, the transient simulation results could provide the estimated SGD under different seasonal recharge rates. The long-term groundwater level monitoring data contributed to the flow modeling approach, making it possible to estimate the SGD under the sparse data condition. However, it should be noted that in the absence of any data on discharges, simultaneous calibration of permeability and recharge leads to inherently nonunique results.
The sensitivity analyses using the steady-state simulation of the autumn 2017 condition reveal that at the model domain scale, changes in groundwater recharge cause larger changes in the SGD rate compared with the aquifer K value and the seafloor conductance. This indicates that groundwater recharge has a larger contribution to the variation in the SGD rate in this area. An increase in groundwater recharge will cause more SGD into the coastal sea. In addition, the weather and groundwater monitoring data indicate a positive correlation between groundwater recharge and precipitation. Based on the future climate change scenarios, precipitation is predicted to increase in winter, whereas evapotranspiration is predicted to increase during summer in southern Finland (Olsson et al. 2015). This would increase the winter time SGD to the Baltic Sea. The sensitivity analysis of the offshore area reveals the low impact of parameter changes on the SGD rates in the low-permeability seafloor areas. However, in the specific high-permeability seafloor regions, such as the areas around the pockmarks B, D and E, the SGD rates decrease as the seafloor conductance at the GHB increases. The increase in conductance of the seafloor could create more uniform flow in a large area, which could reduce the SDG rate in the pockmark areas. This indicates that any changes that enhance the seafloor conductance, e.g. wave erosion, could alter the SGD rates and patterns. In addition, analysis shows the high sensitivity of the SGD rate in the permeable seafloor regions, indicating the need for site investigations and sufficient groundwater discharge data for model calibration and validation.

Conclusions
A three-dimensional geological model of the shallow coastal aquifer belonging to the First Salpausselkä ice marginal formation in the northern Baltic Sea was carried out in Isolähde-Lappohja area, southern Finland. The late Pleistocene and Holocene depositional succession in the study area consists of till, glacial coarse-grained and fine-grained sediments, postglacial fine-grained deposits, and reworked coarse-grained littoral and aeolian deposits. The aquifer is composed of the glacial and postglacial coarse-grained sediments, which were deposited in a NW-SE oriented bedrock depression, in the direction of glacial meltwater discharge. The aquifer is exposed on the shallow shore platform that extends approximately 100-250 m offshore from the shoreline, where the unit slopes steeply seawards and becomes covered by a thick layer of glacial and post-glacial muds. Groundwater flow preferentially takes place in the channel-fill outwash coarse-grained sediments and in sand and gravel interbeds that provide conduits of higher hydraulic conductivity, and have led to the formation of pockmarks on the shore platform edge and slope at water depths between 4 and 17 m. A two-layer groundwater flow MODFLOW model was constructed based on the results of geological and hydrogeological data for both steady-and transient states. The steady-state flow model was run with the calibration of 97 groundwater level data points from the observation wells and GPR profiles measured in autumn 2017 (correlation coefficient R 2 of 0.99 and RMSE of 0.37 m). The transient flow model was calibrated with the daily monitoring groundwater level data taken during February 2019 to April 2020 (R 2 of 0.96 and RMSE of 0.10 m). The simulation results estimated the average SGD rate to the Baltic Sea at 0.22 cm day −1 in autumn 2017. The average SGD rate increased to 0.28 cm day −1 as a response to an approximately 30% increase of recharge in spring 2020. These values are consistent with the previous SGD rate estimates based on 222 Rn measurements in this area. Results of the simple calculation using hydraulic gradient show that the SGD rates vary with the changes of the hydraulic gradients, hydraulic conductivity and porosity of the aquifer media. Furthermore, the sensitivity analyses reveal that at the model domain scale, recharge has a larger contribution to the variation in the SGD rate compared with aquifer K value and the seafloor conductance. While at the local scale in the offshore area, the seafloor conductance in the highpermeability area in the shallow shore platform has a larger impact on the SGD rates than in the low-permeability regions in the deeper offshore area.
The groundwater flow models were successfully developed with sparse data but still contain uncertainties. More site investigations with emphasis on recharge estimation, and groundwater discharge data, are required for the calibration and validation of models in the future.