Ocean Rise Effects on Shallow Groundwater in Coastal Areas: A Case Study in Juelsminde, Denmark

1. Introduction

The global sea level rise is accelerating due to the melting of ice sheets and the thermal expansion of the oceans as a result of climate change [1,2,3]. This will affect millions of people living in low-lying coastal areas worldwide [1,4,5] and many coastal urban areas are already experiencing the consequences of increased precipitation, more frequent storm events and sea level rise SLR [6].

While SLR poses direct challenges to coastal areas, such as flooding and erosion [7,8,9], groundwater inundation (GWI) is an indirect and increasingly problematic consequence of SLR [8]. Ground water level (GWL), especially in the permeable ground, responds to tidal forces. This narrowing of the unsaturated space between the GWL and infrastructure may lead to GWI. Extremely high tides and SLR will intensify these floods due to a higher groundwater table and reduced unsaturated space to store water from storm events [6,10,11,12,13]. There may be many contributors that can influence GWI. It may consist of a single parameter or multiple parameters coinciding, such as rainfall and high tides during a storm event. Furthermore, low-lying coastal zones are prone to compound flooding, where marine and/or surface flooding can happen simultaneously with GWI [4,9]. Coastal barriers designed to protect against surface water are not effective in preventing flooding from below, making GWI especially problematic. This leaves buildings, basements and both underground and surface infrastructure at risk of flooding [6,8,11,12,14,15]. Consequently, it can lead to contamination of the surface and groundwater with sewage [10].

The role of groundwater is important to fully understand the effect of SLR in coastal areas; however, it is often inadequately addressed in studies related to this topic. The effects of SLR are often assessed using numerical modelling, but there are few studies on modelling and forecasting GWI. This lack of information on GWI in coastal areas means that it is not being adequately considered in planning or adaptation initiatives [6,8,9,12,16]. Modelling GWI requires consideration of all parameters that may influence the GWL during extreme weather events, as well as the risk of compound flooding [9]. So far, many of the models used to investigate these systems have been based on physical principles. However, these models require substantial amounts of information about the system, and even when the data is available, they can be difficult to calibrate. Therefore, machine learning models are gaining popularity among hydrologists because they perform well and require less data input [6].

Existing research papers on the effects of SLR on GWL portray a variety of methods and perspectives, reflecting the complexity of the topic.

Habel et al. 2017 [10] used a groundwater flow model to assess the impact of SLR and high tides on GWI in Waikiki, Honolulu, Hawaii and concluded that an SLR of only approximately 1 m would result in a significant increase in the surface area experiencing GWI and its corresponding ramifications. Bjerkelie et al. 2012 [12] also found that SLR causes groundwater rise (GWR). Using a 3D groundwater flow model, researchers discovered that in New Haven, Connecticut, a 0.91 m SLR scenario resulted in a corresponding 0.91 m increase in the near coast GWL. Furthermore, the model indicated that even a GWL ranging from 5.2–7.3 m above sea level responded to the SLR. Results from Knott et al. 2019 [15] using a groundwater flow model showed that, in coastal New Hampshire, USA, there is an estimated mean rise in GWL of 66% of SLR between 0 and 1 km from the coastline. Additionally, the study found that there is a response in GWL up to 5 km inland (3% of SLR) from the coastline. Similar among these three studies is the conclusion that the models must include more parameters and/or data to enhance their reliability [10,12,15].

In Bowes et al.’s 2019 [6] study, two machine learning models, namely long-short-term memory (LSTM) and recurrent neural networks (RNN), were investigated for their effectiveness in forecasting and modelling GWL in Norfolk, Virginia. A comparison of the models showed that LTSM networks are useful in real-time operational forecasts, and the authors predict that GWL forecasting will become a valuable tool in the future management and modelling of coastal flooding.

Although the described studies have addressed this topic in different localities around the world, no systematic research has been conducted on the impact of SLR on GWL in Denmark.

Therefore, this study aims to investigate how the local groundwater table responds to fluctuations in sea level and to uncover regional variations in the dependence on sea level and precipitation in a small town in Denmark. Finally, the predictability of groundwater fluctuations will be tested using a simple machine-learning model based on the information mentioned above.

2. Materials and Methods

2.1. Study Area and Geological Setting

Juelsminde is located on the eastern coast of Jutland, Denmark, on a peninsula bordered by the Kattegat Sea to the east and Vejle Fjord to the south (Figure 1). Juelsminde is prone to flooding and has been identified as one of the top 10 cities in Denmark that are most at risk of flooding. The nearshore areas are dominated by vacation homes, while the harbour is located in the northern part of the city. Further inland, residential houses dominate. A dike with an elevation of 2.6 m above sea level (masl) protects the buildings located south of the harbour from flooding caused by storm surges. The groundwater level is only 0 to 1 m below the ground surface (mbgs), and large parts of the city are drained through a series of ditches, from which water is pumped over the dike into the Kattegat Sea.

In Figure 1, the near-surface soil types are shown. The Juelsminde area is characterised by glacial sediments in an elevated area to the west, covered by postglacial sediments in a lower lying area to the east. Figure 2 displays a geological cross-section based on well data. The glacial till observed in drilling B2 and towards the west in Figure 1 was deposited during the Pleistocene glaciations when the area was repeatedly covered by ice sheets [17]. The postglacial sediments are dominated by marine sand and, to a lesser extent, organic clay (Gyttja). The marine sand was deposited during the Holocene period, when ice melted and retreated, following the latest glaciation. The sea level rose, causing low-lying areas in Denmark to flood, by the combined effect of sea level rise, and the land still being depressed from the weight of the ice sheets. The land has a much slower response to reaching equilibrium than the ocean. The flooding peaked during the Littorina transgression, which occurred approximately 7,500 years ago. This event is responsible for the lower marine sand layer observed in drilling B2. Marine conditions were established and the gyttja layer found in drilling B2 was deposited. The thickness of the lower marine sand is approximately 1 m, but it may be greater in areas where no organic clay is present and a distinct boundary between the upper and lower marine sand can be distinguished. The organic clay layer has a thickness of 3 m and is composed of silt (Gyttja) containing plant remains and shells.

As the land rose following the ice retreat, the shoreline shifted towards its current position. During this regression, the upper layer of marine sand was deposited. It is up to 10 m thick and described as fine- to medium-grained, well-sorted sand containing shells.

Anthropogenic fillings in the area, typically 1–4 m thick, are expected to be local occurrences related to constructions, and as expected the greatest thickness is found in the harbour area.

2.2. Data

Since 2017, six Rotek loggers have been measuring the GWL every 10 minutes. They are located along to coast-perpendicular profiles, with distances from the coastline between 45 and 210 m (Figure 3).

The northernmost profile is located in an open area devoid of significant buildings, whereas the southern profile extends from the beach just south of the harbour and into the urbanised region behind the dike. In the last quarter of 2021, 14 additional IoT groundwater loggers were installed to investigate the dynamics along the coastline and gain a more comprehensive understanding of spatial variation. The tide gauges were installed along the coastline and in residential areas near the coast, as shown in Figure 3. The Danish Meteorological Institute (DMI) conducts sub-hourly measurements of sea level in Juelsminde harbour. Precipitation amounts in Juelsminde are also recorded by DMI. The data is freely available.

2.3. Data Preparation

The various sources of data require different levels of data processing, as shown in Figure 4. The data obtained from the six Rotek loggers were all corrected for barometric pressure by 100%. Where only a few values were missing, they were calculated using linear interpolation. Major outliers were removed, and staggered sections were aligned with the rest of the data series. The sources of sudden jumps (offsets of sections) in the data are unknown, which means that there is a risk of correcting the wrong section of the data. This uncertainty is one of the reasons why this study focuses on the relative changes in the groundwater table. Therefore, the mean of the data series has been subtracted. The data retrieved from the 14 IoT loggers were corrected based on the terrain, and significant outliers were removed by the company that provided the groundwater loggers. The sea level and precipitation data from the DMI required very little correction. Only a few outliers were removed. All data series that were compared were first standardised to the same resolution before analysis, using either linear interpolation or downscaling.

In this study, two different analyses were carried out. (1) A manual analysis was conducted to examine the groundwater responses to high tide and storm events separately. Based on the results of the manual analysis (2), a machine learning model is constructed to assess the predictability of the groundwater response.

Manual Analysis of Sea Level and Groundwater Interaction

Rises in shallow groundwater levels resulting from large, sudden changes in sea level are typically more distinct and easier to identify and separate from each other than minor variations caused by tidal fluctuations. Additionally, they can be traced further inland, making them ideal for evaluating the systemic response of shallow groundwater to sea level rises. The most significant sea level rises are identified and extracted for analyses.

The analysis was conducted on the six data loggers (Rotek), with the longest available time series, placed in two profiles perpendicular to the coast. The sea level rise is determined by subtracting the mean sea level from the level recorded over the past two days. Through automated peak identification, time periods with a peak prominence over 1.2m and SL above 0.5m are selected and listed in Table 1. Information regarding the scale of the SLR and the time interval covering the period of five days before the peak sea level and the subsequent seven days was gathered. This selection resulted in 30 high tide events with relative amplitudes ranging from 1,22m to 1.96m. An example of automatic selection of high-sea-level events is shown in Figure 5.

The minimum threshold of 1.2 m for relative water level change ensures a high enough signal-to-noise ratio inland while also including a sufficient number of high tide events to evaluate possible correlations. The process of peak selection and subsequent analysis methods were inspired by Bowes et al., 2019 [6]. Modified versions of the available code have been used.

2.4. Response Time from Sea Level Rise

To determine the lag time between the rise of sea level and groundwater, an automated algorithm was used to calculate the time from the peak of sea level rise to the peak of groundwater level rise. Each storm event was analysed by extracting the time interval from both sea level and groundwater level data, with the sea level data being interpolated to match the higher resolution of the groundwater data. The raw data was compared to two filtered datasets using two different methods: (1) a Butterworth low-pass filter and (2) a running mean. These methods were employed to eliminate high frequencies, such as tidal changes, from the data. However, filtering affects the peak amplitude and, therefore, was only used to determine the signal delay. Applying a running mean proved to perform best on this dataset. The delay was found using cross correlation of single loggers with the sea level data. All the results were verified manually.

2.5. Attenuation of Groundwater Peak

The inland amplitude damping of the sea level rise indicates the extent to which a storm surge of a certain size affects the groundwater table. This attenuation was calculated for each event by determining the relative increase in groundwater levels at each logger location. The relative increase in GWL was calculated by subtracting the average GWL measured from 48 to 24 hours prior to the GWL peak from the GWL at peak time.

2.6. Predictive Model

A delayed signal from the SLR and corresponding GWR enables the prediction of the GWL based on the sea level. The approach was to use a quick and simple model to gain insight into its predictive capabilities.

The primary objective of developing a predictive model is to estimate future increases in groundwater levels by considering current and forecasted changes in sea levels. The aim was to construct a model that could provide precise forecasts, and thus, the focus was primarily on calculating changes in groundwater levels one day in advance.

Before developing the model, a specific model type was selected. Based on initial tests using the same input in several simple machine learning models, support vector regression (SVR) provided the most accurate results. An SVR-model is a regression model which determines the best fit based on a defined, acceptable error [18]. The model presented in this article is based on a linear kernel.

The model was partly based on historical and real-time data, with sea level and groundwater observations used as input. Precipitation and sea level forecast data provided by the DMI were used. DMI creates weather and sea level forecasts up to five days in advance. These forecasts can be utilised to generate predictions for groundwater levels, and based on these predictions, a warning system can be developed for areas at risk.

The model was trained on data from each of the six loggers individually, and the results were calculated separately for each logger. Hence, in this study, the model functions as six separate models. The model’s predictive capabilities include 24-hour, 3-day and 5-day forecasts. The input parameters were selected in part based on the initial manual analysis of the data and the response of the groundwater.

In the final model, the input parameters are as follows:

• Groundwater level from the past 4 days
• Sea level from the past 4 days and 1 day into the future (using the DMI 1-day forecast)
• Precipitation from the past 24 hours and 1 day into the future (using the DMI 1-day forecast)
• Sine and cosine curve representing seasonality

The model is trained on 60% of the available data, while the validation and cross-validation data each consist of 20%.

3. Results and Discussion

3.1. Groundwater Wave Characteristics

In a homogeneous medium, we expect to find a linear relationship between distance and time for the propagation of a wave peak. Figure 6 displays a box plot representing the observed delay times at the loggers.

The loggers along the northern area, which are not affected by dikes or buildings, appear to exhibit a linear relationship. The two most coastal positions on the southern profile also appear to reasonably follow the trend, while the farthest logger on the southern profile seems to be influenced by other factors. This can be explained by the natural, high variation in the time lag response of the GWL, as the lag time is expected to be dependent on the magnitude of the sea level rise. The peak response in the GWL is clearly visible in the coastal loggers, but the algorithm had difficulty identifying the resulting rise in the loggers located further away from the coast. This is due to a low signal-to-noise ratio in the logger data, which makes it difficult to identify. Moreover, the further inland the loggers are placed, the greater the possibility that the signal from the SLR is drowned out by responses from precipitation events. Furthermore, the further inland, the more likely that local pumping and drainage may also affects the local groundwater level.

Figure 7 and Figure 8 shows that the amplitude of the groundwater level increases as the distance from the coast increases, while the sea level rise decreases. We expect the amplitude damping in a free aquifer to follow a logarithmic function [19]. With a logarithmic time-axis, the relationship can be estimated by fitting a straight line. However, there are outliers, especially when very large sea level rises occur, resulting in extra high amplitudes. Note that event 14 and 25 have outliers caused by possible mispicks. However, there has not been sufficient arguments to exclude them from the dataset. Figure 9, indicates a semi-linear relationship between the rise in sea level and the corresponding increase in the GWL. This relationship is particularly evident in the most prominent SLR events and in Logger 7, which is closest to the coast. Although, the relationship is also detectable in Logger 3 and Logger 6, it is generally significantly less well defined in the two farthest loggers due to signal reduction.

However, these trends are not entirely straightforward. There is significant variation among different events in terms of both the determination of phase shift and the change in amplitude. This is due to the simplification of the groundwater system in Juelsminde and the expectation of complete reliance on the ocean, particularly as the distance from the coast increases. The farther from the coast, the greater the influence of additional factors.

3.2. Sensitivity and Dependency Analysis

The 14 new loggers did not have time series long enough to contain sufficient events for the previous analysis. However, the data were used to conduct an overview analysis of the groundwater’s dependency and sensitivity to sea level rises and precipitation events (Figure 10). The dependence on sea level fluctuations is determined not only by the storm surge responses but also by the trace of tidal changes in the loggers. In these cases, precipitation has a reinforcing effect on the signal, but it is not the primary cause (blue circles). When the GWL cannot be correlated with sea level changes but mainly responds to precipitation, it belongs in the red (triangle) category. In general, loggers located closer to the coast tend to be more vulnerable to the impacts of sea level changes compared to those situated further inland. However, in certain areas, this principle may not apply to individual loggers, and it is presumed that local conditions may influence which parameters have the greatest impact on the groundwater level in a logger. Therefore, it may be difficult to determine exactly how much different parameters affect the various wells, but it can be generally stated that precipitation affects the GWL in all wells to varying degrees, while the impact of sea level is highly dependent on the distance to the coast. The sea level’s dependence is also an indication of the hydraulic connection between the well location and the ocean. Well no. 20 is located relatively far inland, but it clearly responds to tidal changes. It is located adjacent to a wetland that is connected to the ocean, which explains the rapid response.

3.3. Predictability of Groundwater Fluctuations

Results from the 24-hour predictions are displayed in Figure 11 and Figure 12, with the distance to the coast increasing downward. The time period depicted covers a high-water event and the consequent rise in GWL in six wells. The solid grey line represents the observed change in GWL, while the dotted blue line represents the model’s prediction for 1 day ahead. In addition, precipitation and SL fluctuations are visualised below.

The average correlation coefficient for this high-water event is approximately 90%, which indicates the percentage of variations that can be accounted for by the model. The model also performs well, with a high accuracy rate of 78%, in predicting groundwater levels. The average mean deviation is just below 1cm.

The model, however, has its challenges. Three of the major issues are: (1) an overestimation of the contribution of precipitation; and (2) challenges in predicting the extent of extreme GWL rises and (3) relying heavily on the previous groundwater measurements to calibrate the level in the current well. Overestimating the effect of precipitation on the signal can result in periods of a few hours with significantly higher modelled values than those observed. On the contrary, the inability to model extreme events results in modelled values that are lower than the observed values. Both of these factors result in a significant deviation between the predicted and measured GWL. The reliability on previous measurements results in delayed replications of sudden changes and adjust the level retroactively.

Table 3 provides a brief summary of the model’s overall performance in predicting changes in groundwater levels, specifically for 1, 3 and 5 days ahead.

The indicators adjusted R² (R²_adj) and predicted R² (R²_pred) are used to compare the model performance. R²_adj indicates how well terms fit a curve or line, but it adjusts for the number of parameters in a model. This means that if a parameter is included in the model but does not contribute to the model, the R²_adj will decrease. If the parameter contributes to the model, R²_adj will increase. R²_adj will always be less than or equal to R². R² is the residual sum of squared errors divided by the total sum of squared errors.

$$R^2=1-\frac{S{S}_{RES}}{S{S}_{TOT}}=1-\frac{{{\displaystyle \sum }}_i{\left(y_i-{\hat{\mathrm{y}}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_i{\left(y_i-{\stackrel{-}{\mathrm{y}}}_{\mathrm{i}}\right)}^2}$$

(1)

where y_i is the dependent variable for observation i.

$$R_{adj}^2=1-\frac{\left(1-R^2\right)\left(n-1\right)}{n-k-1}$$

(2)

where n is the sample size and k is the number of variables in the model.

R²_pred is also known as PRESS statistics and is used to determine how well a regression model makes predictions. If the data consists of a lot of noise, then R²_pred will be low as it is not possible to predict random noise. R²_pred is useful to avoid overfitting.

$$R_{pred}^2=1-\left(\frac{PRESS}{S{S}_{TOT}}\right)\ast 100=1-\frac{{{\displaystyle \sum }}_i{\left(y_i-{\hat{\mathrm{y}}}_{-\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_i{\left(y_i-{\stackrel{-}{\mathrm{y}}}_{\mathrm{i}}\right)}^2}$$

(3)

where ŷ_(-i) is the predicted value of the response variable for this observation found from the fitted regression equation [20,21].

The numbers show that the modelled values correspond well with the observed values (R²_adj), but there are still difficulties with making independent predictions (R²_pred). However, the low average is due to a few data anomalies that arise from error correction rather than a rising water table. In addition, the average deviation from the observed values is only a few centimetres. When attempting to predict groundwater changes three days into the future, the R²_adj decreases while the ability to make accurate predictions (R²_pred) increases. This may be because the model was originally slightly overfitted, resulting in a greater deviation from the observed data but a less constrained fit of the model when calculating the predicted values. On the other hand, it seems unlikely that the model can generate valuable outcomes when trying to forecast groundwater fluctuations five days ahead. Here, both correlation coefficients are negative.

4. Conclusions

Comparing sea level and groundwater data has demonstrated a need for greater emphasis on addressing flooding caused by groundwater following sea levels rise. Especially in countries with extensive coastlines and low-lying coastal infrastructure. Although the impact of a storm surge is most noticeable in the wells located closest to the coast, the response to the strongest storm surges were recorded in a well located 210 m from the coastline, with four rows of houses between its location and the ocean. The response signal experiences a delay that ranges from 16 hours and 45 m near the coast to more than 30 hours 210 m away from the coast. In several cases, it has been difficult to identify the high-water signal in the loggers located far inland, possibly because the signal weakens logarithmically as it moves away from the coast. This means that the impact of factors such as rainfall becomes significant enough to have a comparable amplitude, which can mask the signal from the sea. Meanwhile, the groundwater level closer to the coast is more heavily influenced by the impact of the sea. The greater the magnitude of a sea level rise, the greater the corresponding, but delayed, rise in the groundwater table.

The delayed response allows for the calculation of the predicted response of the shallow groundwater ahead of its onset. The support vector regression model shows promising results when predicting the groundwater response in a logger 24 hours in advance using sea level and rainfall data. However, it would be interesting to introduce further improvements and move from single-step prediction to predicting entire time sequences. Additionally, expanding the model to a multi-output model calculating the groundwater table in all wells simultaneously would eliminate the need for calibrating for each well.

This study highlights the significance of incorporating the influence of sea level on groundwater variations in the development of predictive models and the examination of future groundwater levels. In certain regions, the interaction between the sea and the shallow groundwater table intensifies, underscoring the relevance of including this factor in analyses. Furthermore, when determining appropriate climate adaptation strategies, it is crucial to account for the natural dynamics associated with these phenomena.