Impacts of Future Climate Change and Baltic Sea Level Rise on Groundwater Recharge, Groundwater Levels, and Surface Leakage in the Hanko Aquifer in Southern Finland

The impact of climate change and Baltic Sea level rise on groundwater resources in a shallow, unconfined, low-lying coastal aquifer in Hanko, southern Finland, was assessed using the UZF1 model package coupled with the three-dimensional groundwater flow model MODFLOW to simulate flow from the unsaturated zone through the aquifer. The snow and PET models were used to calculate the surface water availability for infiltration from the precipitation data used in UZF1. Infiltration rate, flow in the unsaturated zone and groundwater recharge were then simulated using UZF1. The simulation data from climate and sea level rise scenarios were compared with present data. The results indicated changes in recharge pattern during 2071–2100, with recharge occurring earlier in winter and early spring. The seasonal impacts of climate change on groundwater recharge were more significant, with surface overflow resulting in flooding during winter and early spring and drought during summer. Rising sea level would cause some parts of the aquifer to be under sea level, compromising groundwater quality due to intrusion of sea water. This, together with increased groundwater recharge, would raise groundwater levels and consequently contribute more surface leakage and potential flooding in the low-lying aquifer.


Introduction
Groundwater recharge is an important process in maintaining groundwater levels and the sustainability of groundwater resources.Future climate change (e.g., in terms of mean annual temperature and precipitation) would influence hydrological systems, groundwater recharge and groundwater resources.The increasing trend for global warming is predicted to continue in the future and is expected to have greater impacts on hydrological systems, resulting in more vulnerable water resources [1,2].
In permeable, unconfined aquifers located in low-lying coastal areas, changes in sea level and potential future sea level rise may affect groundwater quantity and quality due to salt water intrusion [3][4][5].Rasmussen et al. [6] reported a combination of significant changes in groundwater recharge, sea level rise, groundwater abstraction and the state of drainage canals, due to future climate change, as being crucial for sea water intrusion and groundwater quality in low-lying coastal aquifers west of the Baltic Sea in Denmark.Expected changes in surface leakage due to changes in the relationship between recharge, sea level rise and groundwater level rise are thus of interest.Surface leakage represents groundwater that leaves the aquifer in form of a discharged water to land surface whenever the altitude of the groundwater table exceeds land surface [7].If surface leakage were to increase in the future, this would have direct implications e.g., land use planning.Increased surface leakage could also cause flooding and contamination of wells.For example in Finland, surface water intrusion into wells owing to flooding or to changes in the relationship between groundwater level and surface water level has been observed to increase e.g., the bacterial content in groundwater [8].In addition, sea water intrusion into aquifers can increase salt concentrations or salinities of fresh water aquifers.Therefore there is an urgent need to assess the potential impact of future climate change on low-lying aquifers that are vulnerable to changes in groundwater conditions.
A number of recent studies have assessed the impact of climate change on groundwater using groundwater flow modelling as an assessment tool, and have found that recharge estimation is one of the most challenging parts of groundwater modelling [9][10][11][12][13][14][15][16].Many previous studies have linked the results of recharge models to groundwater flow models in a one-way approach in which time-dependent recharge is calculated first, basically for different soil columns where hydraulic permeability of soil and average groundwater level are known.The time-dependent groundwater recharge value is then implemented into a groundwater model (such as MODFLOW) as a time-dependent flux boundary.This is a valid approach when the groundwater level fluctuation is small compared with aquifer thickness and when the capillary rise is insignificant, which is usually the case in highly permeable glaciofluvial deposits (e.g., eskers in Finland).However, despite the insignificant capillary fringe, the recharge can still be affected by the groundwater level in a shallow groundwater system.If the groundwater level rises close to or at the ground surface, less recharge and more surface leakage and even runoff may occur.This is an important issue which needs to be taken into account when simulating the impact of climate variations and change in very shallow parts of aquifers.Hence, the state of the groundwater level in different time periods needs to be better represented than by the one-way approaches recently developed.One alternative approach could be solving the Richards' equation in three-dimensions (3D).There are a few numerical models available that are able to solve the Richards' equation in 3D, for example HydroGeoSphere (HGS) [17], OpenGeoSys [18], ParFlow [19,20] and FLUSH [21].However, solving the Richards' equation in 3D, especially in a regional-scale model, can be computationally demanding [22] and hence more robust approaches should be used in practical applications.These could involve e.g., coupling UZF1 [7,23] or HYDRUS [24] to MODFLOW [25,26].HYDRUS was developed by coupling the 1D-Richards' equation with MODFLOW in order to balance computational speed and model efficiency, and has been shown to provide an accurate characterization of the vadose zone flow process based on the Richards' equation [27,28].However, coupling UZF1 and MODFLOW was found to provide 25%-50% faster computation time than coupling HYDRUS and MODFLOW [29].Moreover, the HYDRUS1D-MODFLOW model can suffer from numerical oscillations, which requires an extended period to spin-up initial conditions to avoid numerical artifacts.This also makes it difficult to calibrate the HYDRUS-MODFLOW model [30].The UZF1 package solves the 1D kinematic-wave approach for unsaturated flow, by considering only the gravity-driven downward fluxes in the unsaturated zone, and does not account for the effect of capillary force, thus providing more time-efficient simulations than the fully 3D solution of the Richards' equation [22].UZF1 is implicitly coupled to MODFLOW by iteratively balancing the dependence of groundwater recharge and discharge on groundwater head.It considers head-dependent groundwater discharge out of the groundwater system, which allows groundwater discharge to the land surface (such as a stream or a lake) whenever the groundwater level exceeds the land surface elevation.Shallow groundwater can result in groundwater discharge at the surface in UZF1 [28,29].In addition, it has been pointed out that on a regional scale at large space and time discretization, the solution of the Richards' equation such as in HGS may not be more inaccurate than the simplified process of coupling UZF1 with MODFLOW [22].
In recent years, fully coupled approaches have also been developed and used to simulate sea-aquifer interactions in order to determine the density variation in groundwater and sea-aquifer interface.For example, coupled MODFLOW and MT3DMS have been used to simulate the variable-density groundwater flow and solute transport in the SEAWAT computer program [31].Yang et al. [32] used HGS to simulate the effects of tides and storm surges on a coastal aquifer.Another study summarized the processes, measurement, prediction and management of seawater intrusion in a coastal aquifer including a number of computer codes that are capable of simulating seawater intrusion into the coastal aquifer [33].The fully coupled sea-aquifer model in conjunction with a 3D Richard's equation would probably be the most accurate tool to estimate the interaction and the interface between seawater and groundwater.However, the density of the Baltic Sea is relatively low, on average 1.005 kg m −3 , compared with that of normal oceanic water (average 1.025 kg m −3 ) [34].In the Gulf of Finland, the density of seawater varies from 1.001 to 1.006 kg m −3 (average 1.003 kg m −3 ) and salinity from 3.0% to 10.23% (average 6.6%) [35].The salinity value varies with depth, from 3.0% to 5.5% in shallow water [36] up to 10.23% at a maximum water sampling depth of 80 m [36].These values seem to remain stable, as they are similar to values determined for water samples taken more than 30 years ago in this area [37].In future climate change scenarios, salinity in the central Baltic Sea is predicted to decrease by about 2.0%-2.5% [38], with a decrease of 8%-50% from the present for the whole Baltic Sea region [39], due to the predicted increase in freshwater inflow into the Baltic Sea.Moreover, tidal variations are small along the Finnish coast, where the annual seawater level varies between 2.0 and −1.3 m [40], while there is reported to be no tide in the Baltic Sea [41].The climate change scenarios do not predict the direction of change in Baltic Sea density, which makes the prediction of density-dependent models arguably quite trivial.The Baltic Sea level rise and fluctuation estimated for various climate change scenarios can be used in the first instance to estimate the impacts of Baltic Sea level fluctuation and its impacts, together with those of changes in climate variables, on groundwater and interactions between the Baltic Sea and aquifers connected to it.
The main objective of this study was to assess the impact of climate variation and change, and of sea level variation on groundwater recharge, groundwater level and surface leakage in shallow, unconfined, low-lying coastal aquifer in southern Finland.Because of the relatively low density of the Baltic Sea and the lack of predicted data on density changes in this water body under different climate change scenarios, density-dependent models were not used.Instead, the following modelling approach was used in this study: First, the snow and potential evapotranspiration (PET) models were used to assess the infiltration rate to the UZF1 model.Second, UZF1 was coupled to MODFLOW-2005.Changes in sea water level were imported as time-dependent constant head boundary conditions in MODFLOW.The impact of soil frost on water infiltration was considered to be minimal, because according to [42] infiltration into the Hanko aquifer can take place even at a surface temperature of −2 °C.In cold, snow-dominated regions, the impact of soil frost on water infiltration is generally an issue, but in permeable, sandy aquifers recharge can occur even in partially frozen soil [43,44].The groundwater simulation results based on future climate (precipitation, temperature, and sea level rise) scenarios (A1B and B1) were then compared with present day climate conditions, in order to evaluate the potential change in groundwater resources in the future.

Study Area
The study area is located on the Hanko peninsula, on the southern tip of Finland, at approximately 59°20' N 23°05' E, and covers approximately 117 km 2 (Figure 1).The local economy consists of 60% services and 38% industry and the population of Hanko in 2013 was 9282 (www.hanko.fi).Hanko is also a popular summer resort, so the population increases greatly during the summer due to the arrival of holiday home owners and tourists.During the period 1971-2000, mean daily minimum, mean and maximum temperature at the study site during the year was −4.2, 5.7 and 16.6 °C, respectively, and mean annual precipitation was 620 mm [45].Scots pine (Pinus sylvestris) forest is the main land use in the aquifer area.The Hanko aquifer forms part of the First Salpausselkä ice-marginal formation, the most extensive glaciofluvial formation in Finland.The elevation varies between 10 and 14 m above sea level (m a.s.l.) along the proximal (northern) side of the First Salpausselkä Quaternary ice-marginal formation, decreases to <2.0 m a.s.l.along the northern coastline, and on the distal (southern) side gradually decreases in the south and south-east to an average elevation of 5-7 m a.s.l.The distal side is partly covered by swamps and peatland (Figure 1).In the east, the area is covered by sand dune terrain from Aeolian deposits.The main groundwater area consists of Quaternary sand and gravel deposits above the Precambrian crystalline igneous and metamorphic bedrock.The thickness of the Quaternary sediment varies from <1.0 to 75 m (average ~13 m).Groundwater level varies between 2 and 10 m below the ground surface in the inland area to <2 m below the ground surface in the coastal area, where the groundwater discharges into the Baltic Sea (Figure 2).Groundwater recharge occurs mainly twice a year, during spring (late March to early April) and late autumn (November to early December), following infiltration of snowmelt and rainfall.Groundwater flows mainly northward to the coast, but some flows south-southeast to the peatland and swamp area and some east to the Baltic Sea.

Groundwater Level Data
Daily groundwater level data during March 2009 and March 2010 were obtained from 12 monitoring wells along the west coast of the main groundwater area.Historical groundwater level data since 1972 were taken from the OIVA database at the Finnish Environment Institute (SYKE).However, monitoring data were not available for the whole period for all wells, as only 18 of approximately 240 wells have been measured regularly once or twice a month since 1972 and 1974, while the remaining wells were measured in particular periods.

Climate Data and Climate Change Scenarios
Daily precipitation and temperature data during the period 1963-2010, measured at the Tvärminne weather station, were obtained from the Finnish Meteorological Institute (FMI).Snow water equivalent data for the same period, measured at a station in the Santala area, were obtained from SYKE.Trend analysis was carried out for the temperature and precipitation data during the period 1968-2010 using the Mann-Kendall and Sen non-parametric tests [46].Climate model data on two greenhouse gas emissions scenarios, A1B and B1, for the period 2001-2100 and a simulation of present climate 1971-2000 from CLM (Climate Limited-area Modeling community) were obtained from the World Data Center for Climate, Hamburg [47].The A1B emissions scenario predicts a possible future world of very rapid economic growth, global population peaking in mid-century and rapid introduction of new and more efficient technologies, with a balance across all energy sources.The B1 scenario envisages a convergent world with rapid change in economic structures towards a service and information economy and the introduction of clean and resource-efficient technologies.
The data has a grid resolution of 0.2 × 0.2 Decimal Degrees where four data points fall over the Hanko area.Daily precipitation and temperature data were derived from the average values of those four points and represent the climate data over the model domain area.The climate simulation data for the 30 year period 1971-2000 were used for comparison with the measured data for Hanko obtained from FMI.The Delta approach was used for the transfer process between present climate and scenario data [12,44,48,49].For comparison, two climate parameters (temperature and precipitation) from the A1B and B1 scenario data were grouped into two 30 year time spans: 2021-2050 and 2071-2100.

Sea Level Data and Sea Level Rise Scenarios
Daily sea level data for the period 1971-2010 were obtained from Marine Research, FMI.During 1971-2000, the mean sea level varied between −0.49 and +0.58 m, with an average level of −0.02 m.The sea level rise scenarios used were A1B and B1 (high and medium regionalized) from the CLIMBER model, with compensation for the vertical post-glacial land movement component for the period 2000-2100 with the baseline 1995 [50,51].By the end of the 21st century, sea level in the Hanko area was predicted to rise by 22 cm, to 81 cm above the baseline 1995 value depending on the scenario.After calibration with sea level in the Hanko data, the mean sea level in Hanko was +0.09 m in the B1 scenario (medium regionalized) and +0.51 m in the A1B scenario (high regionalized) (Figure 3).By the end of 2050, the mean sea level in Hanko was predicted to be −0.03m in the B1 scenario and +0.13 m in the A1B scenario, while during the period 2021-2050 it was predicted to be −0.10 m and −0.01 m in the B1 and A1B scenarios, respectively.The average predicted sea level for both scenarios was −0.06 m, which is lower than the average for 1971-2000 (−0.02 m).Trend analysis was also carried out for the seawater level data during the period 1971-2010 using the Mann-Kendall and Sen non-parametric tests.

Methods
Recharge estimation and groundwater flow modelling were performed in three consecutive simulation processes.First, the snow and PET models were used to calculate the surface water available for infiltration from the precipitation input data used in the UZF1 package.Then infiltration rate, flow in the unsaturated zone and groundwater recharge were simulated using the UZF1 package run in MODFLOW-2005, whereupon the 3D groundwater flow model MODFLOW-2005 was used to study the impact of climate change on groundwater resources in the Hanko area.The simulated groundwater models from scenarios A1B and B1 were compared with the simulated model from present data , in order to assess the impact of climate change based on those climate scenarios on groundwater resources.

The Snow and Potential Evapotranspiration (PET) Models
The water available for infiltration (UZF1 FINF infiltration rate) was estimated by the following framework: The input parameters were daily precipitation and temperature.Snow water equivalent was used as a calibration target.The precipitation was apportioned between rainfall and snowfall as [52]: where Ptot is the observed precipitation; Cs is the correction coefficient for snowfall; Cr is the correction coefficient for rainfall; and Fr and Fs are the fraction of rainfall and snowfall, respectively.By setting Tmin and Tmax, the threshold temperatures at which all precipitation is assumed to fall in the form of snow and rain, respectively, Fr and Fs were estimated as [52]: (3) The degree-day model, which takes into account freezing of the snowpack, was used to calculate daily snowmelt Md [52,53]: where Md is the daily snowmelt depth; Km is the degree-day factor; Ta is the daily air temperature; Tm is the temperature at which the snow begins to melt; Fd is the rate of freezing in the snowpack, Kf is the degree-day freezing factor; Tf is the threshold temperature for freezing; and e is a coefficient indicating the non-linear relationship between refreezing and temperature [54].
The liquid water retention capacity of snowpack (WSP) is related to liquid water of ice in the snowpack (WIP) as follows [52]: where R is a retention parameter.Mass balance of the snowpack was calculated as: where SWE is the snow water equivalent.
The PET (mm/d) was predicted using a temperature-based method described previously [55].In this method, mean daily temperature (Ta, °C) and day length (D, hour) are needed as input parameters: where ea is saturation vapour pressure (kPa) at mean daily temperature.

Model Calibrations for the Snow and PET Models
The calibration was performed independently of UZF1 and MODFLOW.The calibration target was observed snow water equivalent values during the period 1963-2001.The model was validated during the period 2002-2012.In order to reduce the calibration parameters, the threshold temperatures Tmin, Tm, Tf and the exponent e were set to constant.The calibration parameters were Tmax, Cs, Cr, Km, Kf, and R and they were constrained between values found in the literature [52,53,[56][57][58].The correlation coefficient (R 2 ) between observed and modelled snow water equivalent values was set to 0.64 and 0.72 for the calibration and validation period, respectively.The plot of observed and modelled snow water equivalent values for Hanko 1963-2012 is shown in Figure 4. Mean monthly value of snow water equivalent was calculated from the calibrated data during 1971-2000, which represented present data for comparison with the scenario data.Then, using the same parameter sets, temperature and precipitation from the A1B and B1 scenario data during 2021-2100 were used to simulate the snow water equivalent to be used as the initial infiltration rate (UZF1 FINF) value for UZF1 and MODFLOW-2005.

Unsaturated Flow (UZF1 Package)
Water infiltration into soil was simulated by the UZF1 package [7,23], which is particularly suited to MODFLOW-2005 [26].Its purpose is to simulate vertical water movement in the unsaturated zone, groundwater recharge rate, and storage in the unsaturated zone using a 1D kinematic wave approximation to the Richards' equation that ignores capillary forces.Direct evaporation of groundwater and evapotranspiration are also included.UZF1 applies the Method of Characteristics [59], which simulates unsaturated flow through homogeneous sediment between the land surface and the groundwater level.In this study, it was assumed that soil frost did not affect water infiltration, i.e., the impact of soil frost on hydraulic conductivity was not included.The parameters needed for the UZF1 and MODFLOW-2005 simulations are listed in Table 1.

Groundwater Flow Model and Model Calibration
Model construction and conditions were as described by [60,61].A 3D geological model of the Quaternary deposit in the Hanko aquifer was constructed in order to provide the geological framework for the groundwater flow model using various geological and geophysical data, including gravimetric survey and ground-penetrating radar (GPR) data, sediments from drilled wells, groundwater level data, Quaternary maps, and topographical Digital Elevation Model (DEM) data.The software used was an integration of ArcMap/ArcInfo and the 3D geological module was from Groundwater Modelling System (GMS).The bedrock surface was first identified and then a solid model of the Quaternary deposit was constructed between topographical and bedrock surfaces.The Hanko groundwater flow model was simulated using the finite-difference MODFLOW-2005 [25], under the GMS graphic environment.The model was a single layer model with size 71 km 2 and consists of a rectangular grid of 186 rows and 269 columns.Grid cell size was 50 m × 50 m, decreasing to 5 m × 5 m for cells near the main groundwater intake wells.The top of the model was the DEM surface and the bottom the bedrock surface and the thickness of the model varied from 1 to 75 m (average 25 m).The bedrock has very low hydraulic conductivity and was regarded as an impermeable layer forming the bottom boundary of the model.
The model boundaries were assigned based on the 3D geological model, including specific head boundary, general head boundary (GHB), drain and seepage boundary [62], and no flow boundary (Figure 5).The specified head boundaries are located along the interface of the aquifer and the Baltic Sea coastline.The GHB is located in the north-east of the model domain, to describe the conditions between the aquifer and the external source/sink along the fractured bedrock area.Drain and seepage boundaries are located in the east, south-east, and west coast boundaries, where drain and seepage flows occur.In these areas the bedrock surface is at shallow depth and the thickness of Quaternary sediment is generally less than 1 m.The no flow boundary is located in the impervious bedrock area in the west of the model, next to Hanko town.Groundwater extracted in the study area is used mainly for domestic purposes, drinking water, and water supply to industry.The water intake well in the Hanko area consists of five pumping stations belonging to the Hanko Water and Wastewater Works and other pumping wells owned by local industries and private individuals.The groundwater pumping rate is on average 3900-4630 m 3 d −1 , with the maximum permissible pumping rate being 9620 m 3 d −1 for the public water intake well and about 400-1700 m 3 d −1 for local industries.Similar groundwater pumping rates have been recorded in the database since 1998.The future population in Hanko is predicted to remain the same as at present, although there may be an increased demand for water supply in the future from the local industries.Thus the same groundwater pumping rates as at present were applied in the groundwater flow models of the scenario data.
All available groundwater level data from 240 observation wells including 12 daily monitoring wells during 2009-2010 and 18 monthly monitoring wells during 1972-2010 were used in the calibration process using PEST [63] and trial and error by manually adjusting the K values until the best fit was obtained between the observed and simulated data.The steady-state flow models were first simulated for the groundwater head during April and November 2009, representing wet and dry conditions, respectively.The correlation coefficient (R 2 ) between observed and simulated head from both steady-state flow models was 0.99.The water balance change between outflow and inflow was less than 0.01%.Then the monthly time steps transient flow model was calibrated to groundwater levels during 2001-2010.The RMSE (root mean squared error), MAE (mean absolute error), and model bias were estimated based on [64] and [44].The mean RMSE for residual between simulated and observed heads during 2001-2010 was 0.56 m, MAE was 0.49 m, and bias was −0.11 m. Figure 6 shows the spatial distributions of mean hydraulic head residuals-differences between simulated and observed heads from observation wells during 2001-2010.

Climate Variables and Sea Level
Trend analysis of the temperature and precipitation data during the period 1968-2010 using the Mann-Kendall and Sen non-parametric tests showed increasing trends, with a rate of 0.0419 °C yr −1 for mean temperature and 0.0054 mm yr −1 for mean precipitation.The Baltic Sea level data for Hanko during the 40 year period 1971-2010 showed a decreasing trend of 1.33 mm yr −1 .All scenarios showed an increase in mean temperature of at least 1.0 °C for all periods compared with the present period 1971-2000 (Figure 8a), and the mean temperature increase was higher in winter than in summer.Scenario A1B (2071-2100) gave the greatest increase in mean temperature (range 2.9-4.6 °C), with a mean annual increase of 3.4 °C compared with 1971-2000.The maximum increase in mean temperature (3.7-4.6 °C) took place in winter and the minimum (2.6-2.9 °C) in summer.Scenario B1 (2071-2100) gave a similar mean temperature pattern to A1B (2071-2100), but with approximately 1.0 °C lower mean temperature.In the period 2021-2050, the mean temperature in the A1B and B1 scenarios was similar, with a mean temperature increase of 1.0-1.4°C in summer and 1.3-2.0°C in winter.There was a mean annual increase of 1.3-1.4°C from 1971 to 2000.

Snow-PET Model
The amount of water available for potential infiltration rate (FINF) is presented in Figure 9a.
The simulations indicated increased annual potential infiltration rate for all scenarios, by 5%-34% compared with 1971-2000 (Figure 9b).Seasonally, changes in infiltration rate varied depending on scenario.The B1 (2021-2050) scenario gave the lowest infiltration rate of all scenarios and the lowest infiltration rate occurred in winter and spring, with an approximately 7%-45% decrease from 1971 to 2000.The other scenarios showed increased infiltration rate in all months except April.The highest increase (21%-69%) in infiltration rate occurred in autumn and winter, while in summer (July) the infiltration rate increased by approximately 6%-12% compared with 1971-2000.In April, the A1B (2071-2100) and B1 (2071-2100) scenarios gave a reduction in infiltration rate of 14%-16% compared with 1971-2000.Table 2 summarizes mean monthly initial infiltration rate over 30 years for the present period  and the A1B and B1 scenarios during the periods 2021-2050 and 2071-2100.

UZF1-MODFLOW
The water balance variables obtained from the simulation with the UZF1-MODFLOW-2005 model including UZF1 infiltration, UZF1 ET, surface leakage, GW ET, groundwater recharge, change in sea level, and change in water storage, are summarized in the following sections.The simulation results for UZF1 infiltration, UZF1 ET, GW ET and groundwater recharge showed slight or no differences between the high and medium sea level scenarios and therefore only the simulations from the high scenarios are presented in diagrams (Figure 10a-d).

UZF1 Infiltration and UZF1 ET in the Unsaturated Zone
Simulated average UZF1 infiltration (Figure 10a) over the period 1791-2000 was 478 mm per year, or approximately 70% of FINF (676 mm), and showed the same pattern of seasonal variations as the FINF values from the snow model (Figure 10a).The simulations showed increased mean annual UZF1 infiltration in all climate scenarios, with the increase ranging from 5% to 32% of present (Figure 10a), and with the greatest increase in A1B (2071-2100) and the smallest in B1 (2021-2050).Seasonally, the simulations for most scenarios showed increased UZF1 infiltration, although in B1 (2021-2050) UZF1 infiltration decreased by approximately 16%-40% from present in winter and spring due to the decrease in precipitation in the corresponding period compared with present (Figure 10a).Simulated average UZF1 ET for the present period  was 194 mm per year, with the highest ET of 42 mm in summer (August) and no ET in winter.Simulation results for all scenarios in all periods showed an increase in annual UZF1 ET of between 19% and 42% relative to present, with A1B (2071-2100) giving the highest value and A1B (2021-2050) the lowest (Figure 10b).Seasonally, the simulations showed the highest UZF1 ET in summer (August) and the lowest in winter.However, the greatest change in UZF1 ET (greater than 100%) relative to present occurred in winter.

Evapotranspiration in the Saturated Zone (GW ET) and Groundwater Recharge
Simulated GW ET was 12 mm per year, which was quite small compared with UZF1 ET (194 mm per year).The simulations indicated that in all scenarios except B1 (2021-2050), mean annual GW ET increased by 78%-98% from present, whereas in B1 (2021-2050) it decreased by 25% from present.The pattern of seasonal change in GW ET pattern was similar to that in UZF1 ET, with the simulations generally showing the highest values in summer and the lowest in winter.However, in B1 (2021-2050), the simulations showed a decrease in GW ET of approximately 26%-52% from present in late spring to early autumn (May-October) (Figure 10c).
Simulated groundwater recharge during the period 1971-2000 was 283 mm per year, or approximately 42% of FINF (676 mm).Groundwater recharge at present  takes place twice a year, during late autumn-early winter and during late spring.The greatest recharge occurs in December (44 mm) and the lowest recharge in July-August (4-5 mm).Simulated mean annual groundwater recharge in scenarios A1B and B1 changed by between −14% and +33% compared with present.Simulated recharge in all scenarios except B1 (2021-2050) was higher than at present, with recharge values of 356, 376 and 351 mm, respectively, for A1B (2071-2100), A1B (2021-2050) and B1 (2071-2100), or an approximately 26%, 33% and 24% increase over present.In contrast, B1 (2021-2050) gave a decrease in recharge value of 243 mm, or 14% from present.
Seasonally, the simulations of groundwater recharge showed a change in recharge patterns from present depending on climate scenario and period.During 2021-2050, the groundwater recharge patterns from both the A1B and B1 scenarios were still similar to present, with recharge taking place twice a year, during late spring and late autumn-early winter (Figure 10d).However, based on scenario data, recharge started about a month earlier during autumn compared with present, and A1B (2021-2050) gave increased recharge for all seasons compared with present, with the highest increase in late autumn.In contrast, B1 (2021-2050) gave a decrease in recharge during winter and spring of 13%-55% compared with present.During 2071-2100, the simulated groundwater recharge showed a similar pattern and value in both the A1B and B1 scenarios.Overall, groundwater recharge pattern showed a shift in both scenarios, with recharge starting about one month earlier in autumn and continuing to winter and with higher total recharge compared with present.No recharge took place during late spring (April).Instead, after its highest increase in winter, the recharge continued to decline in late spring and there was low recharge in summer, with summer starting about one month earlier than at present (Figure 10d).

Surface Leakage
The simulation data showed a 18%-65% increase in mean annual surface leakage compared with present in all scenarios except B1 (2021-2050), which gave an approximately 10% decrease compared with present.Seasonally, the simulated surface leakage showed similar seasonal variation patterns to recharge, with high surface leakage during the recharge periods in autumn and winter and low surface leakage in summer (Figure 10e,f).In addition, the simulated surface leakage showed the influence of sea level rise.During 2021-2050, the difference in sea level between the high and medium scenarios was small, with predicted sea level varying between −0.03 and +0.13 m, and the simulated surface leakage based on those sea level values showed differences of less than 6% (Figure 10e).During 2071-2100, when the predicted sea level was +0.53 and +0.17 m for the A1B high and medium sea level scenarios, respectively, and +0.37 and +0.09 m for the B1 high and medium scenarios, respectively, the simulated surface leakage differed by 23% for the A1B scenarios and 16% for the B1 scenarios (Figure 10f).

Relative Change in Groundwater Level and Baltic Sea Level
A positive inflow indicates a higher sea water level than groundwater level at the boundary and sea water flow into the model domain or the aquifer.On the contrary, a negative outflow indicates groundwater flow out of the model domain because of a higher groundwater level compare to sea water level.During 1971-2000, the mean sea level in the Hanko coastal area varied between −0.49 and +0.58 m (average −0.02 m).Based on the A1B (high) and B1 (high) scenarios, by the end of 2100 mean sea level was predicted to be at +0.51 m and +0.37 m, respectively.Relative changes between the simulated groundwater level and Baltic Sea level gave a 27.3% and 19.7% increase in inflow of sea water to the aquifer in the A1B (high) and B1 (high) scenarios, respectively, relative to present (Figure 11a).The lowest relative change between simulated groundwater level and the Baltic Sea level, −16.7%, was found in the A1Bmed (2021-2050) scenario (Figure 11a), where the mean sea level was predicted to be −0.03 m.This indicates more groundwater flow out from the model domain than at present.The simulated surface leakage showed a positive correlation with both groundwater recharge and inflow, which in turn had a direct response on sea level.Table 3 summarizes the correlation coefficients between surface leakage and recharge, and surface leakage and inflow.Figure 12 presents the spatial distribution of changes in surface leakage rate under climate scenarios A1B and B1 and their high and medium (med) sea level scenarios during the periods 2021-2050 and 2071-2100 compared with present .

Groundwater Level
Simulated groundwater levels in the A1B and B1 scenarios were measured from grid cells at the observation wells.Mean groundwater level over 30 years for each scenario for the periods 2021-2050 and 2071-2100 was compared with present .The simulated groundwater levels based on scenarios A1B and B1 and their high and medium sea level scenarios showed different groundwater patterns for wells located in the middle of the aquifer and wells located in the coastal area.In the middle of the aquifer (Figure 13a), the simulations showed a decrease in groundwater level only in the B1 (2021-2050 high and medium) scenarios, of approximately 0.81 and 0.85 m, respectively, from present.The other scenarios showed increases in groundwater level of approximately 0.77-0.98m from present, with the A1B (2021-2050 high) scenario showing the highest increase in groundwater level (0.98 m).In the coastal area (Figure 13b), the simulations also showed a decrease in groundwater level only in the B1 (2021-2050 high and medium) scenarios, by approximately 0.16 and 0.23 m, respectively, from present.The other scenarios showed increases in groundwater level of approximately 0.13-0.39m from present.The changes in groundwater level in the coastal area were a direct response to changes in sea level, with the highest change in groundwater level for A1B (2071-2100 high), the highest sea level, and the lowest for B1 (2021-2050 medium), the lowest sea level.Figure 14 shows the spatial distribution of changes in groundwater level for climate scenarios A1B and B1 and their high and medium sea level scenarios during the periods 2021-2050 and 2071-2100 compared with present.

Discussion
The transient groundwater flow model applied here, utilizing UZF1 coupled with MODFLOW-2005, provided information not only on groundwater recharge into the aquifer system, but also on surface water and groundwater interactions, such as groundwater discharge to sea water and to low-lying areas.It thus provided a more realistic picture of the groundwater recharge process, taking into account the position of the groundwater level in every time step, than previous studies [9][10][11][12][13][14][15][16]44].Understanding and information on this process are important for sustainable water resource management and vulnerability assessments in low-lying aquifers.
Assessment of climate impacts on groundwater recharge and levels using numerical groundwater modelling as a prediction tool contains many kinds of uncertainties, some of which are associated with the modelling techniques and methods used to predict and calibrate the numerical model to fit with geological data in flow modelling.There are also uncertainties deriving from the different kinds of future climate scenarios used [16].The simulation results showed increased groundwater recharge in all future climate scenarios except B1 (2021-2050) and showed a strong correlation between groundwater recharge and the amount of precipitation, compared with the other variables, e.g., temperature and evapotranspiration.This may indicate that precipitation has a greater impact on groundwater recharge than temperature.However, an increase in winter rainfall rather than snowfall was caused by the increased temperature during winter for both future climate scenarios (A1B and B1).The results of the snow model also indicated a decrease in snow water equivalent, while the amount of precipitation increased.Similar results have been reported previously [44,48,65,66].The large predicted increase in temperature, particularly during winter and spring, had a great impact on recharge pattern during 2071-2100.The increase in temperature in winter and spring caused more snowmelt earlier and a major shift in infiltration period, starting earlier in early spring.With a reduction in precipitation and a continued increase in temperature and evapotranspiration during late spring (May), there would be a reduction in groundwater recharge and groundwater resources during summer.
Changes in recharge in terms of quantity and recharge pattern affected other water balance components, including surface leakage, whereby groundwater was discharged to surface drainage in low-lying areas, altering the interaction of sea water and groundwater along the coastline and water storage.The simulations showed positive correlations between groundwater recharge and surface leakage.The shallow aquifer in Hanko is relatively flat and in many parts of aquifer the shallow groundwater level is currently close to the ground surface.With a predicted increase in groundwater recharge due to increased snowmelt and winter rainfall under a warmer climate, the groundwater level will rise even more to breach the ground surface threshold and discharge into the low-lying area or surface water bodies as surface leakage.This suggests an increased risk of potential flooding during high recharge from autumn to winter and early spring.
Trend analysis of the sea level data for the Hanko area was carried out for the period 1971-2010 using the Mann-Kendall and Sen non-parametric tests [46,61].Assuming a linear trend, the sea level in the Hanko area decreased by 1.33 mm yr −1 over the last 40 years .This value is in contrast to the regional sea level rise of 1.7-2.07mm yr −1 reported elsewhere [67,68].It has to be considered, however, that the time series taken into account in this study started in 1971 and was rather short.Based on sea level observations starting in 1888, the decreasing trend in Hanko and in the whole of Finland was considerably slower after 1970 than before.The long-term sea level values in the Baltic Sea can be well explained by the factors of post-glacial land uplift and changes in global mean sea level, while the recent deviation from a steady linear trend can be explained by the water balance of the Baltic Sea and a strong correlation with the North Atlantic Oscillation index [69].
Sea level had an impact on the relationship between inflow and outflow along the Hanko coastline.A high influence of sea level rise on groundwater level was observed in wells along the coastline.In contrast, the observation wells in the middle of aquifer area showed small or insignificant influences of sea level on groundwater level.Previous studies on the impacts of future sea level rise on groundwater quality and groundwater resources have reported that these will mainly take place along the coastline, main rivers and drainage canals in low-lying aquifers [5,6,70,71].Those studies also reported the influence of stream types and tides on the groundwater level of aquifers along the coastline, with higher groundwater levels in response to sea level rise causing more discharge into the streams.In the present study, sea level rise had an impact on the groundwater level and surface leakage in the Hanko aquifer.For instance, in response to a rise in sea level, the groundwater level was 0.6-0.7 m higher in well 1 (approximately 60 m from the coastline) than in observation well 2 (approximately 120 m from the coastline).The highest increase in sea level (0.51 m a.s.l.), predicted by the end of 2100 in the A1B high (2071-2100) scenario, would push the groundwater level along the coastline to increase by approximately 0.39 m above the present level.Increases in groundwater recharge and sea level will increase the groundwater level, which will increase surface leakage and cause to more discharge low-lying inland areas and along the coastline.In addition, a rise in sea level to 0.51 m a.s.l. by the end of 2100 (based on the A1B-high scenario) would cause some parts of the aquifer area to be under sea water, especially at observation well 1, where the well head elevation is approximately 0.58 m a.s.l.This result agrees well with findings in a previous study that in the Hanko area sea level changes have a direct and very fast effect on groundwater level [72].Another previous study also reported positive correlations between seawater level and groundwater level in wells located in the coastal area in Hanko [61].
Previous studies suggest that sea water intrusion may result in contamination of the groundwater and a reduction in aquifer size [3,4,71].In the Hanko area, although there is evidence that seawater pushes the groundwater level in wells located at the coastline, the physical and geochemical data for those wells show insignificant contributions of seawater from seawater intrusion into the aquifer [73].This may be because of the low salinity of the Baltic Sea [36], which may not pose a high risk of seawater intrusion compared with that in other parts of the world.However, the salinity concentration exceeds the drinking water standard and in the long term or during storm surges over the low-lying area [33], flooding or inundation of seawater can result in saltwater (with higher chloride concentration) infiltrating into the soil and/or damaging infrastructure in the area.The Hanko area experienced the highest sea water level since 1887 (+1.24 m a.s.l.) during a storm surge on 9 January 2005, while pushed the sea level high above the normal sea level.With increased climate change and sea level rise in the future, the frequency of storm surges is also expected to increase, which can pose a risk to groundwater quality and the water supply system, e.g., pipeline corrosion.Increased groundwater recharge in the future would cause the groundwater level in Hanko to increase closer to the ground surface or discharge to low-lying areas.A shorter infiltration depth or an existing wetland area would pose a contamination risk for groundwater quality.A high groundwater level would pose a risk to infrastructure and water supply systems.
Although most of the eight climate scenarios (except B1 (2021-2050)) showed increased groundwater recharge relative to present, six of the scenarios showed changes in water storage ranging from −4.3% to +3.8% from present and the other two B1 scenarios (2021-2050) showed decreases in water storage of −17.8% from present.Changes in water storage less than 5% are probably within the range of uncertainty, but a −17.8% decrease from present indicates a critical impact of climate change on groundwater resources.
The approach used in this study could be applied in the other low-lying coastline aquifers such as in Finland, where approximately 300 of the total classified 6000 Finnish shallow groundwater areas locate less than 100 m from the Baltic Sea shoreline [73], or in the Baltic Sea Regions e.g., Denmark or Lithuania [74].

Conclusions
This study assessed the impact of future climate variations and sea level rise on groundwater recharge and groundwater level in the shallow, unconfined, low-lying coastal Hanko aquifer in southern Finland, using groundwater flow modelling as an assessment tool.The unsaturated zone flow UZF1 package was coupled with MODFLOW-2005 to simulate the flow from the unsaturated zone through the aquifer, using the infiltration water values produced by the snow and potential evapotranspiration (PET) models.The A1B and B1 climate and sea level rise (high and medium) scenarios were applied for two periods, 2021-2050 and 2071-2100.The water balance in the model domain was computed and compared against the reference period of present data .
The simulations showed the impacts of changes in groundwater recharge and sea level rise on groundwater storage and surface leakage in the Hanko aquifer.Changes in the groundwater recharge pattern from present were predicted in the A1B (2071-2100) and B1 (2071-2100) scenarios, with peak recharge not occurring during late spring but earlier, during winter and spring.A future increase in winter temperature could increase winter rainfall, cause more snowmelt, and reduce the snowpack, which would allow more water to infiltrate into the soil and further enhance groundwater recharge during winter and early spring.However, a continued increase in air temperature and evapotranspiration during late spring and a reduction in precipitation would cause a reduction in groundwater recharge and pose a risk of water shortages during summer.The B1 (2021-2050) scenario gave a −17.8% decrease in water storage from present, while the other scenarios gave changes of less than 5%.Any true decrease in water storage could cause critical water shortages in Hanko.
The simulated surface leakage showed strong correlations to groundwater recharge and sea level.Increased groundwater recharge and sea level in a future climate would increase groundwater levels and cause more surface leakage to low-lying areas, increasing the risk of flooding from surface overland flow during winter and early spring.Moreover, rising sea level in the future would cause some parts of the aquifer to be under sea water, with the associated risk of compromising groundwater quality.

Figure 1 .
Figure 1.Location and Quaternary deposit map of the study area.A cross-section along line A-A' is presented in Figure 2. Contour line labels indicate meters above mean sea level (m a.s.l.).

Figure 3 .
Figure 3. Average predicted sea level in the Hanko area according to the A1B and B1 (high and medium regionalized) scenarios, 2020-2100.Simulation results for two future 30 year periods (2021-2050 and 2071-2100) were used for comparison with current 30 year data (1971-2000).

Figure 5 .
Figure 5. Boundary conditions of the model domain and locations of the observation wells.

Figure 6 .
Figure 6.Spatial distribution of mean residuals-differences between simulated and observed heads from observation wells during 2001-2010.

Figure 7 .
Figure 7. Comparisons between simulated (solid line) and observed (dashed line) groundwater levels for the selected observation wells, 1971-2000.Locations of observation wells (Obs) are shown in Figure 5.

Figure 8 .
Figure 8. Change in (a) mean monthly temperature and (b) mean monthly precipitation in scenarios A1B and B1 for the periods 2021-2050 and 2071-2100 compared with 1971-2000.

Figure 12 .
Figure 12.Spatial distribution of changes in surface leakage rate (%) in climate scenarios A1B and B1 and their high and medium (med) sea level scenarios for the periods 2021-2050 and 2071-2100 compared with present (1971-2000).

Figure 13 .
Figure 13.Changes in mean annual groundwater level in simulation results for the A1B and B1 scenarios and their high and medium sea level scenarios during the periods 2021-2050 and 2071-2100 compared with present (1971-2000) at wells (a) in the middle of the Hanko aquifer; and (b) along the coastline.

Figure 14 .
Figure 14.Spatial distribution of changes in groundwater level for climate scenarios A1B and B1 and their high and medium (med) sea level scenarios during the periods 2021-2050 and 2071-2100 compared with present (1971-2000).

Table 1 .
UZF1 and MODFLOW-2005 parameter used in the study.

Table 3 .
Correlation coefficient (R 2 ) for the cross-plot between surface leakage and recharge, and surface leakage and inflow.