Hydrogeological Characterization of Coastal Aquifer on the Basis of Observed Sea Level and Groundwater Level Fluctuations: Neretva Valley Aquifer, Croatia

Abstract

Hydrogeological data availability is often limited to local areas where usual in situ tests or methods are applied (slug/bail or pumping tests, Electrical Resistivity Tomography (ERT)). Because most problems (e.g., saltwater intrusion mitigation) require problem analysis on larger scales (catchment or sub catchment), hydrogeological identification of global character is preferable. This work leads to the determination of aquifer hydrogeological parameters on the basis of observed sea level, groundwater piezometric head found inland, and barometric pressure. When applied to observed signals, the approach led efficiently to final hydrogeological characterization. After identification of dominant tidal constituents from observed signals, barometric efficiency was successfully determined. Following available information on geological settings, an appropriate conceptual model was applied and updated to count for polychromatic signals. Final determination of hydrogeological parameters relied on root mean square error (RMSE) minimization and led to determination of (i) presence of three stratigraphic units: unconfined sandy aquifer on the top, a confining layer made of clay, and a confined gravel layer; (ii) existence of the clay layer under the sea with a total length of 1400 m; (iii) a clay layer has been identified as confining one by both spectral analysis and determined leakance value; and (iv) estimated confined aquifer specific storage ranging from 2.87 × 10⁻⁶ to 4.98 × 10⁻⁶ (m⁻¹), whereas hydraulic conductivity ranged from 7.0 × 10⁻⁴ to 7.5 × 10⁻³ (m s⁻¹). Both range intervals corresponded to previous in situ findings conducted within the area of interest.

1. Introduction

A wide variety of methods to determine aquifer parameters have been established and conducted on real sites in the past. Starting from the simplest and most cost-effective slug or bail tests, over pumping tests, or borehole drilling combined with laboratory analysis of samples taken from cores, local hydrogeological conditions are usually determined. Compared to slug/bail tests as a low cost-way of determining aquifer parameters and enabling hydraulic conductivity estimation, pumping tests offer more spatially global insight into aquifer parameters such as transmissivity and storativity [1]. Borehole drilling gives reliable insight but with high demanding cost, and is often used as first order data needed for verification of aquifer heterogeneity and geological structure when applied together with geophysical methods [2]. Geophysical methods are most commonly used at an up-to-date level and result in effective and temporally efficient application with easily adjustable spatial domain coverage [2]. Despite their attainability, due to economical demands, they are mostly used for small-to-medium scale areas, giving locally characterized information [3].

The first basis on tidal methods to define coastal aquifer parameters were introduced by Jacob [4,5] and Ferris [6] during the 1940s and 1950s. From that point on, these methods were accepted by the hydrogeological community and were further developed to increase their applicability potential and the reliability of results. Carr and Van den Kamp [7,8] demonstrated a successful application of the tidal methods to conduct the coastal aquifer hydraulic diffusivity in the 1970s. In its origin, the tidal method stems from the partial differential equation (PDE) of transient front propagation through the coastal aquifer, with seaside boundary condition defined as long-term sea level oscillations. General solution for the piezometric head oscillation within the confined aquifer induced by changes in sea level can be derived easily in the analytical form, as shown by many authors [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].

The simplest solutions have been derived with the following assumptions: (i) the aquifer is homogeneous and isotropic; (ii) the beach is assumed as a vertical line, perpendicular to aquifer geological layers; (iii) the effect of the outlet capping and subsea aquitard or semi-permeable layer roof extension are neglected; (iv) loading efficiency is set up to unity; and (v) flow is assumed as 1D, perpendicular to the vertical coastline [13,17]. Improvements in the method have been implied as follows: (i) leakance effect through the semi-permeable overlying layer [10,11,12,20,21,28], (ii) confining layer extending under the sea [10,11,15,19,21,28], and (iii) outlet capping influence [16,20,23,26,29].

Jiao and Tang [17] derived a solution for the coastal aquifer conceptual model by involving the capacity of the solution to capture the leakance effect, mainly for tidal efficiency $TE$ and time lag $\Delta T$ determination. Their solution was upgraded by Dong et al. [13] to enable the use of a real sea level signal instead of its dominant constituents. Geng et al. [14] investigated the effects of outlet capping elastic storage and leakage on piezometric head signal induced by tides. Results emphasized the influence of outlet capping thickness to piezometric head features. Influence of both leakage and storage of the semi-permeable layer on the tide-induced groundwater fluctuations have been analyzed by Li and Jiao [20]. The solution by Wang et al. [25] generalized several existing solutions and emphasized the importance of water table loading effects in enhancing the amplitude and reduction of phase lag of the piezometric head inland. Additional generalization of the piezometric head solution and its transient nature was presented by Chuang and Yeh [12]. The role of outlet capping when combined with multilayered coastal aquifer was investigated by Guo et al. [16]. In 2019, Bakker [9] presented a novel analytical solution to define transient nature of groundwater head within the multilayered aquifer. Li and Jiao [21] demonstrated the importance of leakage effect, separately for offshore and inland aquifer systems. Their results indicated a vice versa effect of leakage on groundwater level fluctuations inland and offshore. Confined aquifer thickness has been shown to have a crucial role in the determination of groundwater signals. Hereby, Cuello et al. [30] demonstrated that the wedging drastically increases the amplitude of the groundwater piezometric head inland. In 2007, Xia et al. [26] presented an analytical solution of piezometric head inland, which incorporated several parameters in tidal methods: (i) subsea roof length extension, (ii) permeability of the outlet capping, (iii) outlet capping thickness, (iv) variable loading efficiency ($LE$), and (v) partial effects of the leakance within the offshore and inland sub-domains.

Besides the derivation of appropriate tidal method-based solutions, their application has been recognized as a crucial when applied for aquifer hydrogeological characterization. Tidal methods were applied to Konan aquifer in Japan [31] to determine phreatic aquifer parameters. On the basis of two monitoring wells, the study elaborated effects of tidally induced fluctuations 1 km inland from the shoreline. Potential of assessment of transmissivity and storage from time lag and tidal efficiency has been demonstrated for Beihai aquifer in China [32]. In [33], the authors presented a framework to assess the hydraulic diffusivities of both confined and unconfined aquifers. Beihai aquifer was of interest once again with the multilayered conceptual model [34] with submarine outlet capping included. Chen et al. [35] applied a solution [36] to determine coastal aquifer parameters (hydraulic diffusivity, beach slope, and aquifer thickness). Beihai coastal aquifer has been further investigated using tidal methods in the work by Xun et al. [37]. Hereby, the authors demonstrated the possibility of work with more than one monitoring well simultaneously. Analysis of Oualidia aquifer in Morroco was elaborated in [38], where hydraulic diffusivities were determined via tidal methods separately for spring and neap occasion.

Although most of the abovementioned papers were focused upon monoharmonic tidal solutions, due to the PDE linearity, final solution can be expressed as a linear superposition of separate contribution of each dominant sea level frequency constituent, thus enabling the work with multiharmonic signals.

Piezometric head signals observed at the site consist of several components, of which some are aperiodic and some periodic. Although periodicity stems from tidal mechanisms, sea or earth tides, and periodic barometric influence, aperiodic phenomena are mainly caused by incident hydrological conditions and tectonic activity, among other factors [15,39,40]. Independently of the observation technique used, monitored time series are a result of a simultaneous effect of both tidally induced changes and barometric changes, thus defining water level found within the well open to atmosphere. Because some of the tidal methods imply the knowledge of loading efficiency [9,14,26], barometric efficiency, $BE$, has to be determined to characterize tidally induced changes inland from the shoreline. The topic of excluding the effects of barometric pressure from observed groundwater level was of interest in several publications [41,42,43,44,45,46,47,48,49,50,51], whereas $BE$ determination can be approached by time domain or frequency domain analysis [50].

Finally, available publications offer a wide range of demonstration of the tidal method’s capability in defining aquifer conceptual models and its parameters such as hydraulic diffusivity, transmissivity, hydraulic conductivity, storage, and confinement conditions, among other factors [47,52,53,54,55,56,57].

In this paper, we present an analysis of observed piezometric head induced by both sea level and barometric pressure at the River Neretva Valley aquifer in Croatia, and finally define the following hydrogeological parameters of the aquifer: hydraulic diffusivity, specific storage, hydraulic conductivity, subsea extending roof length and aquifer confinement degree. Following past and present in situ investigations and features arising from existing monitoring systems, an appropriate conceptual model was selected. Furthermore, on the basis of previous findings [15,26,41,46] and their partial improvements, we present an integrated and systematic procedure of hydrogeological parameter determination when observed signals are constituted of both barometric pressure and sea level tidally induced fluctuations. To our best knowledge, this is the first paper dealing with the problem of aquifer characterization by involving available geological data and observed time series of sea level, barometric pressure, and groundwater piezometric head inland when applied to the River Neretva Valley aquifer.

2. Materials and Methods

2.1. Study Area

River Neretva Valley is located in the very southeastern part of Croatia, at the coast of Adriatic Sea. Until the 1950s, the area was mainly a marsh area in the central valley, and the coastline was located approximatively 1 km inland from today’s shoreline position. In the late 1950s and early 1960s, the first attempts of melioration were performed through two dominant phases. The first phase considered the building of the embankment called Diga at the location of today’s shoreline. After the successful end of the first phase, the second phase meant the construction of a system of pumping stations fed by water by a large number of connected channels to enable the forced delineation of the water level within the area of interest and to ensure the groundwater level to be found beneath the pedological layer. Those measures were assumed to enable the positive agricultural conditions and consequently enlarge crop productivity.

In its present state, as shown in Figure 1, the River Neretva Valley covers 4500 hectares of cultivated area. To the north it is bordered by karstic hills of a maximum 650 m.a.s.l. height, as well as in its southern and southwestern part. The Adriatic Sea represents the valleys western border where the shoreline fully follows Diga. The overall length of the shoreline is equal to 2200 m, representing the area between two rivers, respectively, Neretva to the north and Mala Neretva towards the south. This area, starting from the sea, is known as Opuzen–Ušće. The Vidrice area is located south of the River Mala Neretva and is bounded by karstic hills on its eastern and western side. The only corridor passing out of the karstic area is a 1500 m long corridor close to the town of Opuzen, called Luke. Terrain elevation within the area of interest ranges between −0.50 and +1.90 m.a.s.l. with reference to “Nula Trsta” vertical datum. The mid part of the valley is divided into two areas, Crepina and Jasenska.

The River Neretva is characterized by typical inter-seasonal discharge fluctuations. With an average precipitation of 1200.78 mm year⁻¹ observed at Ploče meteorological station and 1219.91 mm year⁻¹ observed at Opuzen meteorological station for the period 2009–2018, discharge was found between 40 m³ s⁻¹ up to 2700 m³ s⁻¹, as observed after the presence of extended duration rainfalls at the river catchment area, located mostly at the territory of neighboring Bosnia and Herzegovina. The usual range of discharge found at the study area ranged from 40 m³ s⁻¹ up to 1800 m³ s⁻¹. The abovementioned lower limit is caused by controlled evacuation of the biological minimum 47 km upwards from the Opuzen at the hydro power plant Mostar. The mean long-term average annual discharge corresponds to 290–330 m³ s⁻¹. Significant differences in discharge values are present during the hydrological year. The summer season is identified as the dry season, with on average no more than 100 mm month⁻¹ of rainfall, mainly concentrated within a few days, thus representing isolated rainfalls. The discharge dominating the period from mid-May until late September slightly ranged from 40 to 200 m³ s⁻¹. Isolated discharge peaks were found mainly in period of November to April, as a consequence of intensive rainfall along the catchment area. These peaks mainly went above 1600–1800 m³ s⁻¹.

The present climate is typical Mediterranean with significant differences in precipitation and atmospheric temperature in summer and winter conditions. The minimum hourly temperature observed in the area during the 10 years period from 2009 to 2018 was −7.80 °C, whereas maximum observed temperature was equal to 38.80 °C. Inter-annual climate conditions can be characterized as follows: (i) winter period is characterized by the lowest temperature values followed by the highest percentage of precipitation observed; (ii) what characterizes the summer season is extremely dry periods with high temperatures, especially during July and August.

The first attempts to characterize the geological settings of the valley stem from the middle of the last century, which precedes the melioration epoch. In the 1960s, the first in situ geological works were performed [58,59,60]. Due to the limited technical capacities and constraints in the fulfillment of the full geological characterization, the results are shown to not be reliable after comparison with up-to-date results [2,61,62,63,64].

The results shown in Figure 2 [2,61,62,65,66,67] imply the presence of five basic geological layers: (i) a layer made mostly of fine sands, with local presence of less than 1.0 m thick sub-layers of clay and mud. This layer is present all over the study area and mimics the upper unconfined aquifer with an average thickness of 7–10 m; (ii) below the unconfined sandy layer with on average of 18–20 m thickness, a compact clay layer is found. This layer is characterized by a lower thickness of 10–12 m at the east, close to the Opuzen town area, whereas its thickness increases towards the sea where it is detected with 25 m thickness. The roof of this layer does not vary significantly, whereas thickness changes influences its bottom elevation. A compact clay layer continues beneath the seabed towards onshore from the shoreline; (iii) confined layer made of mostly coarse gravel with the partial presence of mud and fine sands. The depth of layers of the roof varies from 20 m to the east to 35 m to the west; (iv) in the central part of the study area at approximatively 40–45 m depth, one finds a conglomerate layer whose thickness ranges from 1 to 3 m. It is important to emphasize that this is of local character without being detected at all by available boreholes, thus emphasizing its discontinuity; (v) below this layer, there is an extension of layer (iii), mostly made of coarse-to-medium gravel, fulfilled by mud, which represents a compact mixed gravely confined aquifer with an overall depth of up to 130 m unless bedrock is found. Bedrock has been identified with geolectrical sounding [2] and by boreholes [61,62] at the absolute depth of zero at the edging valley area, to a maximum of 160 m below ground level found at the very central part of the valley, known as Crepina.

2.2. Monitoring System

The monitoring system of the piezometric states and salt water intrusion parameters was established in 2009. Up-to-date improvements have been made by performing new gauges, installing new piezometers, and enabling real-time data acquisition. For the purpose of this study, the monitoring system was performed to capture the sea level state and piezometric head states found within the piezometers along the study area. Sea level states were observed every 60 min, collecting one value per hour. THALIMEDES OTT was installed at the sea level gauge station to capture long-term sea level oscillations. The measuring range of the gauge was ±19.999 m with a resolution of 0.001 m and accuracy of ±0.002 m. Within the piezometers at the area of interest, TE Connectivity TRUBLUE vented gauges with a 0 to 300 psi range and accuracy of ±0.05% Total Error Band (TEB) were installed. All gauges were temporally synchronized relative to the unique data acquisition system. Piezometric head values were observed at three locations, piezometers D1, D2, and D4, with the same sampling frequency as the sea level. For the D1 piezometer, overall piezometer depth was set up to 37.50 m, with a perforation height of 1.20 m within the confined aquifer. D2 piezometer overall depth was equal to 34.39 m, with a perforation height of 3 m within the gravelly media. For the D4 piezometer, the overall depth was equal to 25.24 m, with a perforated depth of 2.0 m within the confined aquifer. For all piezometers covered by the study, the filter layer was built of material of higher conductivity than natural material at the location in order to avoid signal attenuation. Barometric pressure hourly values were observed by the local meteorological station Ploče, operated by the National Hydrometeorological Institute.

In total, three datasets were used within this paper, as follows: (i) 1 August to 31 August 2015 for the purpose of the hydrogeoleogical aquifer characterization, (ii) 1 July to 31 July 2015 for results verification, and (iii) 30 April to 4 May 2019 a 4 day dataset was performed with 5 min sampling time to conduct the convergence test of the sampling and averaging time analysis and its influence on the major piezometric head signal parameters. All observed values are shown with reference to “Nula Trsta” vertical datum.

2.3. Time Series Analysis

Temporal changes in observed signals can be expressed as [39,44]

h_OBSERVED = h₀ + p_A + h_TIDAL + ε (L),

(1)

where h_OBSERVED mimics the piezometric signal as observed by the gauge; $h_0$ is a trend found within the observed signal including the changes of the piezometric head due to the extreme hydrological conditions, or simply isolated tectonic activity, human-induced activities such as groundwater pumping system operation, and rainfall or evapotranspiration; $p_A$ is the effect of the barometric pressure after detrending procedure, whereas $h_{TIDAL}$ represents sea tides and ε represents noise. Following the monitoring sampling time, we assumed the high frequency noise was negligible compared to other relevant components of the signal; thus, the observed signal was defined as a superposition of three components, as follows: (i) trend caused by aperiodic changes, (ii) effect of barometric changes to observed piezometric head, and (iii) sea-induced tides. To eliminate for trends from the observed signal, a polynomial fit was applied following the least squares method. For all signals covered by the study, sixth order polynomial fit was detected to obtain the lowest root mean square error (RMSE). To eliminate for long term variations, the observed signal was primarily detrended and scaled to zero mean value so that the obtained signal was defined as

$$h_{OBSERVED}{}^{\prime }=h_{OBSERVED}-h_0\approx p_A+h_{TIDAL} \left(L\right).$$

(2)

The signal $h_{OBSERVED}{}^{\prime }$ is constituted of both sea tide and barometric pressure components and was used for further analysis from this point on. To further elaborate, the observed barometric pressure was primarily divided into two components, $p_A^{\prime }$ considering long term barometric pressure changes corresponding to the trend, and $p_A^{\prime \prime }$ representing pressure changes excluding the $p_A^{\prime }$; thus, the following equation can be written: $p_A=p_A^{\prime }+p_A^{\prime \prime }$. Long-term periodicity in barometric pressure is usually explained as the consequence of the spatial displacement of air mass [42,46,56], whereas $p_A^{\prime \prime }$ captures diurnal and semidiurnal barometric changes due to the effects of air heating and cooling following the transition from day to night and opposite, as well as the all changes, except for the trends [45,50]. In this manner, detrended groundwater head signal was defined as temporally dependent:

$$h_{OBSERVED}{}^{\prime }\left(t\right)=p_A{}^{\prime \prime }\left(t\right)+h_{TIDAL}\left(t\right) \left(L\right).$$

(3)

Transformation of the $h_{OBSERVED}{}^{\prime }\left(t\right)$ from the temporal to spectral domain was done via application of Fast Fourier Transform (FFT) to obtain amplitude spectra characteristics and separate sea tide effects from the changes in piezometric head caused by barometric effects.

After the main tidal constituents were determined from $h_{OBSERVED}{}^{\prime }\left(t\right)$, $h_{TIDAL}\left(t\right)$ was calculated as a linear superposition of sine functions whose total number corresponded to the number of tidal spectral constituents $n$ with appropriate parameters being [18,28,37], amplitude a, period $T$, and phase $\theta$, defined for each constituent denoted by $j$:

$$h_{TIDAL}\left(t\right)={{\displaystyle \sum }}_{j=1}^n{a}_j\sin \left[\frac{2\pi }{T_j}t+{\theta }_j\right] \left(L\right).$$

(4)

The difference $h_{OBSERVED}{}^{\prime }\left(t\right)-h_{TIDAL}\left(t\right)$ was actually the residual found within the detrended observed signal relative to tidally induced change in groundwater head, and corresponded to the influence of barometric changes. This was the main premise in terms of defining barometric efficiency, $BE$, by applying different methods ranging from the spectral to the temporal domain [45,50].

2.4. Tidal Efficiency, Time Lag, Barometric Efficiency, and Loading Efficiency

Tidal oscillations found inland from the shoreline remained sinusoidal with the period identical to the driving force (sea level), but with exponential decrease in tidal efficiency ($TE$) and linear increase of time lag ($\Delta T$) [5,6,68]. Hereby, $TE$ was calculated directly from the observed signals as the ratio of the standard deviations of piezometric head and sea level. The advantage comes from using full signals for $TE$ determination, rather than only maximum amplitudes or peak values, as suggested by [31,33,68,69]. To determine $\Delta T$ values, the $h_{OBSERVED}{}^{\prime }\left(t\right)$ signal was amplified by $TE$ and shifted over the time span unless the highest correlation between the sea level and piezometric head was observed [69]. By amplification, we considered

$$h_{\Delta T}{}^{\prime }\left(t\right)=\frac{h_{OBSERVED}{}^{\prime }\left(t\right)}{TE} \left(L\right),$$

(5)

whereas $\Delta T$ was determined on the basis of the minimization of the term

$${\displaystyle \sum }{\left[h_{\Delta T}{}^{\prime }\left(t\right)-h_{TIDE}{}^{\prime }\left(t-\Delta T\right)\right]}^2 (-),$$

(6)

where $h_{TIDE}{}^{\prime }\left(t\right)$ fully corresponds to the definition of $h_{OBSERVED}{}^{\prime }\left(t\right)$ but refers to sea level instead of piezometric head.

Barometric efficiency, $BE$, introduced in 1940s and 1950s [4,5,6] is successfully kept up-to-date to represent the effects of the barometric pressure changes to water level in piezometers tapping the aquifer. If the external stress represented by barometric pressure is transmitted to the aquifer, reaction stress can be divided between the sediments that the aquifer is made of and water contained within the aquifer pore volume. Because $BE$ represents the percentage of stress that is accommodated by the sediment matrix [43,44,50], a zero $BE$ value corresponds to no response case, whereas a unity $BE$ value respresents a fully efficient response to barometric changes. In case of a zero $BE$ value, total barometric pressure-induced stress is accommodated by the pore water, whereas for a unity BE value, total stress is accommodated purely by sediment matrix.

Within this paper, $BE$ was determined by applying the method of Rahi [41], which was shown to be a modification of Clark’s method [42]. The overall procedure to determine $BE$ can be explained as follows:

$$\Delta b_i=b_i-b_{i-1},$$

(7)

$$\Delta h_i=h_i-h_{i-1},$$

(8)

$$IND=\Delta b_i \times \Delta h_i,$$

(9)

$$S_b^i=S_b^{i-1}+\left|\Delta b_i\right| if IND<0 and \left|\Delta b_i\right| >\left|\Delta h_i\right|,$$

(10)

$$S_h^i=S_h^{i-1}+\left|\Delta h_i\right| if IND<0 and \left|\Delta b_i\right| >\left|\Delta h_i\right|,$$

(11)

$$i=1,\dots ,n$$

where $\Delta b_i$ is the incremental difference of the barometric pressure, $\Delta h_i$ represents incremental difference of the piezometric head, $IND$ is a secondary variable, $S_b^i$ is a sum of barometric changes, and $S_h^i$ is a sum of piezometric changes. Summation was done only for physically interpretative barometrically induced changes in piezometric head [41] up to a maximum n of arranged pairs of barometric and piezometric increments. Finally, barometric efficiency was calculated as

$$BE=\frac{S_h^n}{S_b^n}.$$

(12)

Once BE was determined, loading efficiency, LE, could be defined following [7,70] as

$$LE=1-BE.$$

(13)

2.5. Determination of Aquifer’s Hydrogeological Parameters

Available results on geological structure are updated with geoelectrical sounding [2] to verify previous results locally (Figure 2). Results are also verified with available borehole logs to identify final geological settings. This leads to the determination of the conceptual model on the basis of geological layer identification consisting of three dominant layers: unconfined sandy aquifer, impermeable clay layer, and confined gravel layer, plus bedrock. Geological data offer a good insight into the clay layer, which extends under the sea but with an unknown extending length. The seabed is covered by fine sand (sediments generated by the River Neretva), whereas geometrical properties with clear sea level boundary conditions are found at the western boundary of the valley, at Diga. On the basis of available data, an appropriate conceptual model was selected [26]. To define piezometric head inland, a modification of Equation (13) by Xia et al. [26] was used. The total number of constituents used in modified Equation (13) from [26] corresponded to the number of tidal constituents found within the observed sea level signal, thus enabling the use of polychromatic signals. Generated results of the piezometric head data series were compared to h_TIDAL (t) to obtain the minimum RMSE for combinations of parameter values.

Globally, the conceptual model deals with 11 parameters, as follows: $LE$, $BE$, $D$, $K_P$, $m$, $L$, $L_S$, $S$, $K$, $b$, and $x$, where $K_P$ represents the outlet capping conductivity (m s⁻¹) and $m$ presents its thickness (m), $L$ is the subsea extending roof length (m), $L_S$ is the leakance (day⁻¹), $S$ holds for storativity (-), $K$ is hydraulic conductivity (m s⁻¹), $b$ is aquifer thickness (m), and $x$ mimics the distance between the piezometer location and the shoreline.

During May 2019, three samples of sand were collected from the sea bottom and used for a Darcy experiment to define $K_P,$ whose average value was 3.1 × 10⁻⁴ m s⁻¹. Outlet capping thickness is defined as 7 m on average, starting from the coastline and following bathymetric features along the aquatorium on the western side of the study area.

Starting from 11 unknowns and by using the hydraulic diffusivity expression:

$$D=\frac{K \times b}{S}=\frac{K}{S_S} ({\mathrm{m}}^2 {\mathrm{s}}^{-1}),$$

(14)

the number of unknowns was reduced to four: $D$, $L$, $L_S$, and $x$.

The problem was solved by combining different values of the abovementioned parameters until the RMSE reached its minimum value between (i) the signal generated by the modified Equation (13) from [26], and (ii) signal of tidally induced piezometric head filtered from the effects of barometric pressure.

3. Results

Datasets used for the purpose of this study are shown in Figure 3. Following the sea level observations, piezometric head values within the piezometers D1, D2, and D4 were observed with the same sampling frequency as the sea level. Barometric pressure shown in Figure 3 represents a salt water column equivalent of the signal observed during August 2015. In general, all the time series used within this paper were (i) observed for same duration period, (ii) observed by same observation frequency, and (iii) observed continuously without any missing value.

Although the sea level time series is characterized by the largest amplitudes and shows mixed semidiurnal character, the piezometric head states found within the piezometers D1, D2, and D4 showed significant similarities with the sea level signal with exceptions in the lower amplitude and presence of temporal delay in peak values relative to the sea level signal. Besides these two major differences, which are explained as the aquifer’s diffusivity filtering effect on the piezometric head observed inland from the sea level boundary condition, Figure 3 reveals differences in mean values. This is in accordance with previous findings arising from hydrogeological investigation at the study area and determined gradient of the piezometric head [65]. All piezometric head values are expressed with reference to “Nula TRST”. As can be seen, D4 was characterized by the highest piezometric head, whereas D1 corresponded to the lowest one. The mean D4 piezometric head for the dry period was equal to 0.52, D2 was equal to 0.32, whereas D1 was equal to 0.17 m.a.s.l., thus emphasizing dominant groundwater flow direction from the east and southeast towards the west and northwest, thus defining a gradient towards the shoreline. On average, the piezometric head characterizing the dry period at the area of interest for all boreholes used for geoelogical settings determination is shown in Figure 2. As can be seen, the deeper aquifer shows its artesian nature, thus implying the confinement conditions.

Insight into the barometric pressure implies the typical functional relationship between sea level and barometric pressure, identified at the scale of several days and more. This aperiodic feature in barometric pressure stems from the spatial displacement of the air mass over the location of interest during time period of a few days to a week or even more [41,46]. When compared to sea level and piezometric head, the increase on barometric pressure led to a decrease of both sea level and piezometric head and vice versa.

To enable separation of the effects of the barometric pressure and tidal components from available time series, firstly a long-term trend was recognized. Following Figure 3, time series were detrended to eliminate present trends. The usual frequency upper limit to eliminate for trends ranges from 0.5 cycles per day [56] to 0.8 cycles per day [46,50]. Hereby, we did not limit the cutoff threshold, since the sixth order polynomial fit was applied to minimize RMSE when fitting the observed signal, and thus the cutoff was not same for different time series.

Figure 4 gives insight into the detrended and scaled signals of sea level; D1, D2, and D4 piezometric head; and barometric pressure. After the detrending procedure was completed, one could observe the presence of both neap and spring tide within the sea level signal, as well as the transition from diurnal to semidiurnal tidal components, which guarantee tidally induced phenomena were not lost from the signal due to the detrending and scaling process.

Total observed amplitude during neap tide is approximately 50% lower compared to spring tide amplitude. The general feature of periodic piezometric head signals from piezometers D1, D2, and D4 was explained via tidal efficiency, $TE$, and time lag, $\Delta T$. $TE$, as defined in the previous section, fully corresponded to the definition of [7,46,70], and was interpreted in that manner as the ratio of the standard deviations of the piezometric head inland and sea level, whereas $\Delta T$ corresponded to temporal delay of piezometric signal relative to that of the sea level. Following the results in

Table 1. Tidal efficiency, $TE$, and time lag, $\Delta T$, values determined for detrended signals and signals made of four dominant tidal constituents for the period 1 August to 31 August 2015.

Time Series	Detrended		Tidal Constituents
Piezometer	TE (-)	ΔT (h)	TE (-)	ΔT (h)
D1	0.558	1	0.557	1
D2	0.340	2	0.332	2
D4	0.320	3	0.311	3

, D1 was characterized by TE equal to 0.558, whereas $\Delta T$ corresponded to 1 h. By going more inland, $TE$ was decreased, whereas $\Delta T$ increased to 2 h for D2 and 3 h for D4. The lowest $TE$ value was 0.32 and corresponded to D4. The obtained values implied the effect of the aquifer’s hydraulic diffusivity to attenuate the sea level inland from the shoreline and to increase its temporal delay with the increase of distance from the shoreline. Besides the inspection of signal parameters, visual inspection revealed a correspondence in piezometric head and sea level signal in (i) the presence of neap and spring tide amplitudes, (ii) effects of quadrature and syzygy, and (iii) similar periodicity phenomena, which demonstrate the interconnection between sea level and inland piezometric states. Compared to sea level and piezometric head, barometric pressure after detrending still showed the presence of longer than day periodicity feature, larger fluctuations and deflection from its mean value.

After detrending was done for the raw data series observed at the location of the sea level gauge; barometric pressure gauge; and piezometers D1, D2, and D4, in order to eliminate the effects of long-term trends, the observed time series were transformed to frequency domain by applying FFT, as shown in Figure 5. This enabled an insight into the dominant harmonic constituents of the analyzed variables. Barometric pressure was identified by the presence of semidiurnal S2 and diurnal K1 constituents corresponding to changes of temperature (day–night, heating–cooling) and, consequently, the pressure changed. Besides S2 and K1 constituents, barometric pressure amplitude spectra emphasized the presence of low frequency periodicity.

Following the observed spectral components of detrended sea level and piezometric states found within the piezometers, the aquifer indicated the dominant influence by the sea level, as was previously identified from visual inspection of the signals, their $TE$ and $\Delta T$. To support this additionally, it was mandatory to start with the periodicity features of the sea level. Following Figure 5, sea level was characterized by four dominant constituents, two semidiurnal M2 and S2, and two diurnal ones, O1 and K1, respectively: principal lunar semidiurnal constituent (M2), principal solar semidiurnal constituent (S2), lunar diurnal constituent (O1), and lunisolar diurnal constituent (K1). Although amplitudes and periods were defined directly from the spectra, phase lags for each analyzed signal were determined by applying the least squares method of computation package MATLAB (Mathworks Inc., Natick, MN, USA) and are presented in

Table 2. Harmonic constituent variables for the detrended and scaled sea level and piezometric head signals for the period of 1 August to 31 August 2015.

Variable	Sea Level				D1				D2				D4
Harmonic Constituent	M2	S2	O1	K1	M2	S2	O1	K1	M2	S2	O1	K1	M2	S2	O1	K1
Amplitude (m)	0.0942	0.0729	0.0201	0.075	0.0519	0.0387	0.0161	0.0427	0.03	0.0214	0.0086	0.0283	0.0277	0.0198	0.0084	0.0268
Period (h)	12.400	11.999	25.649	24.008	12.400	11.999	25.649	24.008	12.400	11.999	25.649	24.008	12.400	11.999	25.649	24.008
Phase (°)	−29.55	−3.33	159.80	−93.45	−54.77	−27.66	129.30	−117.50	−91.10	−64.24	96.01	−145.20	−108.80	−83.90	84.71	−156.10

. Although a limited duration signal of 31 days was used, the dominant constituent parameter values observed were in correspondance to typical values for the southeastern Adriatic mareographic stations Ploče and Dubrovnik [71] and contributed to 95% of the full tidal signal at the area of interest [46].

From spectral function features of piezometric head observed at piezometers D1, D2, and D4, one could observe the same period values compared to those of the sea level with differences in (i) the reduction in amplitude for each constituent and (ii) phase lag values compared to ones representing sea level signal. Although same period values imply the sea level as a driving force mechanism causing piezometric state dynamic nature of the aquifer, thus emphasizing connection of the aquifer with the sea, general reduction of the amplitudes was noted as decreasing with increase of the distance inland from the shoreline, which corresponded to the definition of $TE$. The effect of the time lag between the observed sea level and piezometric signals is also visible in Table 2 through the phase lag values of constituents, leading to the conclusion that it follows the generally linear increase of the time lag with increased distance from the shore [4,5,6]. Lunar diurnal constituent (O1) resulted in the largest $TE$ values for all three piezometers covered by the study. The following one was the lunar semidiurnal constituent (M2), whereas the $TE$ of solar harmonic constituents strictly underestimated the $TE$ calculated from the signal. This can be explained as the dominant barometric influence to S2 and K1 solar constituents.

General spectrum features found at all piezometers covered by performed analysis were utilized to identify the type of aquifer in the River Neretva Valley. Although significant dominance of the principal lunar semidiurnal constituent M2 is increased with increase of aquifer confinement [52], additional criteria, such as constituent dominancy and constituent presence criteria can be used to identify the aquifer type. Our findings emphasize dominancy of the principal lunar semidiurnal M2 and lunisolar diurnal constituent K1 in the piezometric head amplitude spectrum with the presence of both principal solar semidiurnal constituent S2 and lunar diurnal constituent O1. Because the low permeable clay layer with on average 18–20 m thickness was identified throughout the available in situ experiments and investigations at the site, a confinement of the deeper aquifer layer was assumed. The confinement hypothesis of the aquifer stemmed from (i) available boreholes from the area of interest [61,62,66], (ii) interpretation of ERT profiles [2], (iii) interpretation of the geolectrical sounding [2], and finally, (iv) performance features of the observed piezometric head signals within the study area. For all three piezometers located inland, one could observe identical spectral features with amplitude reduction with an increase in distance from the sea. The presence of the four signal constituents: two diurnal—O1 and K1—and two semidiurnal—M2 and S2—where M2 dominated over the other constituents, led to the conclusion that D1, D2, and D4 piezometers are penetrating the confined aquifer. These findings are in agreement with previously published studies [41,57].

To eliminate tidal effects from the detrended and scaled signals, identified tidal constituents found within sea level and piezometric head signals were used. With four dominant constituents identified, sea level tidal changes were described by linear superposition of four sine functions by using Equation (4), following their parameters—amplitude, period, and phase shift, as shown in Table 2. The difference between the detrended signal and the four constituents signal was identified as residual value, which was assumed to correspond to the influence of barometric pressure changes characterized by the frequencies higher than the frequency cut-off eliminated by the detrending procedure. Piezometric head residuals and barometric pressure detrended values are shown in Figure 6a–d.

The idea of barometric efficiency determination relies on the Rahi method [41], which is a physically based improvement of Clark’s method [42]. Following the features of signal residuals, the Rahi method was applied to determine the $BE$ value for each of the monitoring piezometers at the site. The results are shown in Figure 7 and give good insight into the aquifer’s water level response to changes in barometric pressure. Because the usage of vented gauges enables insight into real water table fluctuations caused by both effect of increased loading from the surface or tides, which both cause increase of the aquifer pore water content internal pressure and the changes of barometric pressure, the available datasets can be used to determine BE. The obtained values varied in a minor to negligible way in between the three analyzed piezometers, ranging from 0.24 for D1 piezometer to 0.23 found at D2 and D4 piezometers. The values were determined for the largest significant correlation obtained for variable time lag between the barometric pressure and piezometric head residuals. Time lag for piezometer D1 corresponded to 1 h, the same as for D2, whereas no time lag was observed for D4. Following the definition of loading efficiency and its functional relationship with $BE$ [7,70], after $BE$ was known, $LE$ values were equal to 0.76 (D1) and 0.77 (D2 and D4).

The definition of $LE$ enabled the piezometric head time series definition. Hereby, we started from the known $LE$ value. A laboratory test to determine the hydraulic conductivity value of the sandy material found at the seabed corresponding to the outlet capping material was previously conducted and elaborated within Section 2. To finalize the conceptual model with its parameter values, several assumptions were made, as follows: (i) extended roof length was unknown, but existence was identified with available in situ documentation [65], and (ii) due to the evidence of the clay layer presence all over the domain with an on-average layer thickness of 22 m and specific spectrum features found in Figure 5, leakance was identified as negligible, which corresponded to aquifer confinement conditions. Considering assumptions made, Equation (13) from [26] was generalized to be used when four dominant tidal constituents were incorporated to finally determine outlet capping length, exact leakance value, hydraulic diffusivity, and distance from the shoreline. This was done by applying different values of the abovementioned parameters unless the lowest RMSE value was observed when the generated signal was compared to the tidally induced piezometric signal. In between analyzed piezometers, D2 showed the lowest RMSE value, corresponding to 0.003541, obtained for D = 1500 m² s⁻¹, $L$ = 1400 m, and $L_S/S$ = 1 day⁻¹, whereas $x$ was equal to 4080 m, representing the real distance of the piezometer inland from the shoreline. The determined leakance value could be assessed when the storativity value was known. Hereby, we relied on previous findings stemming from the central part of the study area. Below the compact clay layer, one finds a 10–15 m thick coarse gravel layer whose storativity ranges from 10⁻⁶ to 10⁻⁴ according to the literature [1]. This implies a leakance value of the same order of magnitude as the storativity value, which can be classified as negligible and furthermore confirms the aquifer’s confinement.

The RMSE value obtained for D1 was equal to 0.004579 and corresponded to $D$ = 470 m² s⁻¹, $L$ = 1400 m, and $L_S/S$ = 1 day⁻¹ obtained for x equal to 305 m. The lower hydraulic diffusivity value in comparison to D2 was caused by two factors: (i) smallest hydraulic conductivity of the aquifer layer piezometer penetrates to, and (ii) lower aquifer thickness (only 1.20 m) compared to more than 60 m at D2 [62,63,64]. For the sake of completeness, D4 was represented by the same methodology but with one difference in the conceptual model compared to D1 and D2, relying on the fact that the karstic bedrock located between the sea and the D4 location is identified as significantly fractured and with the presence of high conduits [65] representing a highly permeable border. In this manner, outlet capping was still used with the parameters $K_P$ and $m$, as defined previously. The lowest observed RMSE was equal to 0.004258 and corresponded to $D$ = 235 m² s⁻¹, $L$ = 0 m, and $L_S/S$ = 0 day⁻¹, which implies absolute confinement of the aquifer in this part of the domain and emphasizes the validity of the assumption regarding extended roof non-presence. The $x$ value obtained for the lowest RMSE value was equal to 2500 m, which represents the actual distance between D4 and the inner karstic area border southwest of the piezometer location. Good agreements and consequently low RMSE values of tidal induced piezometric head time series are demonstrated in Figure 8a–c.

By the use of solely tidal components of the sea level in tidal methods, we performed a validation of obtained parameters when they were applied to the period of 1 July to 31 July 2015. The same hydrogeological parameter values, as determined for the period of August 2015, were used to generate tidally induced piezometric head at piezometers D1, D2, and D4. In the same manner as for August 2015, the observed sea level was analyzed in the frequency domain and the main constituents were identified. After identification of tidal constituents, we generated tidally induced piezometric head by using modified Equation (13) from [26] within the piezometers D1, D2, and D4. The results shown in Figure 9a–c emphasize good matching between the modeled time series and the time series obtained from measured signals by identification of only tidal constituents. The RMSE values determined were as follows: D1 RMSE = 0.005724, D2 RMSE = 0.007939, and D4 RMSE = 0.00458. Although a slight increase in RMSE values can be seen, still good matching was evidenced by visual inspection of Figure 9a–c.

4. Discussion

To further elaborate our findings and their generality, in order to emphasize the tidally induced method’s applicability, we selected several issues, discussed below: (i) selection of appropriate aquifer conceptual model; (ii) the use of piezometric head residuals and barometric pressure detrending; (iii) additional interpretation of hydrogeological parameters $D$, $L$, and $L_S/S$; (iv) potential influence of surface water to piezometric head found within the confined aquifer; and (v) result sensitivity to different sampling time and temporal averaging scales.

Throughout the literature there is a wide range of available coastal aquifer conceptual models with appropriate solutions to define tidally induced transient piezometric head inland [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29], covering the possibility of taking into consideration a large number of parameters, namely, leakance effect, effects of elastic storage and leakage of the submarine outlet capping, extension of the confinement layer under the sea, permeability and thickness of the outlet capping, variable sea water density effects, and existence of multilayered systems. To select an appropriate conceptual model, we used available results obtained from in situ investigation works, which yielded clear stratigraphy of the aquifer. Following the results conducted in the past, we upgraded aquifer characterization by combination of ERT and geoelectrical sounding [2], which further confirmed and clarified the geological settings. Insight into bathymetric features and available borehole from the sea, offshore from Diga, undoubtedly confirmed the existence of a clay layer below the sea bottom, being covered by a sandy outlet capping layer.

The use of piezometric head residuals instead of full piezometric signals is in accordance with the general definition of observed signal defined by Equation (1) of this paper. Because the observed piezometric signal after elimination of long-term trends was assumed to be constituted solely of tidally induced sea level changes and barometric pressure, the residual represented the effect of the barometric changes to the observed piezometric signal. During long-term trend elimination, typical low frequency filters are not used, and thus all signals covered by this study were not subjected to same filtering procedure in a way of frequency cut-off definition. However, these signals were used within the method of Rahi [41,47,72] to determine $BE$. By deeper inspection of two additional settings made by Rahi, compared to the $BE$ determination method by Clark [42], it was obvious all nonphysical data pairs were excluded from the $BE$ determination process, thus enabling the use of detrended barometric pressure and piezometric time series residuals, as explained in his paper. Results by Turnadge et al. [50] demonstrated good agreement in $BE$ value determination by different methods (including Rahi) when time series were used at the duration interval of up to 1 month (Figure 4 of [50]).

Although only the hydraulic diffusivity was determined within the Results section, hereby we extend the range of hydrogeological parameters and their assessed range values. The expression to define aquifers specific storage is

$$S_S=\frac{{\rho }_W \times g \times \varnothing \times C_W}{1-TE} \left({\mathrm{m}}^{-1}\right),$$

(15)

where $S_S ({\mathrm{m}}^{-1})$ is aquifers specific storage, $g (\mathrm{m} {\mathrm{s}}^{-2})$ is gravitational constant, $\varnothing (-)$ stands for porosity, $C_W \left({\mathrm{Pa}}^{-1}\right)$ presents pore water compressibility, $TE (-)$ is tidal efficiency, and ${\rho }_W$ represents water density, as defined by [73].

Table 3. Hydrogeological parameter estimation.

Piezometer	ρW (kg m−3)	g (m s−2)	∅ (−)	CW (Pa−1)	BE (−)	SS (m−1)	k (m s−1)
D1	1023.00	9.81	0.15	4.58 × 10−10	0.24	2.8727 × 10−6	0.0014
D2	1019.00	9.81	0.15	4.58 × 10−10	0.23	2.9859 × 10−6	0.0045
D4	1016.00	9.81	0.15	4.58 × 10−10	0.23	2.9771 × 10−6	0.0007
D1	1023.00	9.81	0.25	4.58 × 10−10	0.24	4.7878 × 10−6	0.0023
D2	1019.00	9.81	0.25	4.58 × 10−10	0.23	4.9765 × 10−6	0.0075
D4	1016.00	9.81	0.25	4.58 × 10−10	0.23	4.9618 × 10−6	0.0012

summarizes the results when Equation (15) was used to calculate $S_S$, as well as the assessment of hydraulic conductivity using Equation (14). As can be seen, for the range of expected porosity values (0.15–0.25) [1,74], the obtained specific storage ranged from 2.87 × 10⁻⁶ to 4.98 × 10⁻⁶ (m⁻¹), whereas hydraulic conductivity ranged from 7.0 × 10⁻⁴ to 7.5 × 10⁻³ (m s⁻¹). The obtained values corresponded to fine-to-medium gravel, as found during in situ investigation [2,61,66]. Because the perforated screen of the monitoring piezometers were vertically positioned to similar absolute depth, and by insight into available boreholes, this layer did not change its depth significantly, and narrow intervals of specific storage and hydraulic conductivity values were expected. Further in situ investigation works should be focused in order to at least determine transmissivity and storativity values to increase the accuracy of the assessed values.

Obtained leakance values expressed by $L_S/S$ implied very low leakance values, thus increasing the level of aquifer confinement and dominancy of tidal influence. On the basis of the specific storage and local aquifer thickness values, the highest storativity value was assessed to 6.45 × 10⁻⁴, which led to $L_S$ = 6.45 × 10⁻⁴ (day⁻¹). Even in cases of high $L_S$, as shown, leakance through the clay layer and its effect upon piezometric signal characteristics was minor to negligible, thus confirming aquifer confinement [17,69]. Clay layer as a roof extension length obtained for D1 and D2 was equal to 1400 m away from the coastline. From the bathymetric features at the area off from the shore, depth of approximately 25–30 m was identified at the distance of 1400–1500 m away from the shore. This can be used as qualitative confirmation of roof extension length if the roof of the gravel layer is found at approximately 25–30 m below the mean sea level [58,59,60,61,62,63,64,65,66,67].

Although piezometers D2 and D4 are located in the vicinity of surface water bodies, preliminary analysis on the influence of surface water oscillations to the confined aquifer piezometric head time series did not show its relevance. Because piezometer D4 is located 80 m away from Mala Neretva river, we checked for similarity in features of (i) D4 piezometric head time series and (ii) Mala Neretva surface elevation time series. Inspection of these time series implied that Mala Neretva surface elevation did not follow typical tidal features. This can be explained as the consequence of downstream and upstream gate closure during dry season to mitigate salt water cline penetration upstream, thus strengthening the hypothesis that confined aquifer piezometric head is dominantly caused by direct connection to the sea. Similar findings hold for piezometer D2, whereas D1 is located closest to the shoreline.

Because the study relied on available datasets observed at the same sampling frequency, we performed the inspection of different sampling and averaging time scales. From 30 April to 4 May 2019, a 4 day dataset was performed with 5 min sampling time to conduct the convergence test of the sampling and averaging time analysis and its influence upon major piezometric head signal parameters $TE$ and $\Delta T$. Firstly, sea level was used as observed with a 5 min sampling scale, whereas piezometric head sampling time was variable (5, 10, 15, 30, and 60 min). Secondly, sea level was used as defined on a 60 min sampling scale, whereas piezometric head sampling time was variable (5, 10, 15, 30, and 60 min). Parameters $TE$ and $\Delta T$ were defined in a same manner as explained within Section 2.4 and shown in Figure 10a–c. Independently of the sea level signal sampling time, tidal efficiency showed less than 0.4% relative difference, even when 60 min sampling time was used. Apart from $TE$, $\Delta T$ was found to be more sensitive to sampling time. Figure 10b implies a 10 min sampling time to reduce relative difference in $\Delta T$ below 2%, with hourly sampling time relative error equal to 10% for D1, whereas D2 and D4 showed a satisfactory difference (below 3%). If, instead of sampling time, an averaging time was analyzed, it would be demonstrated that the obtained difference is double the one obtained for the equivalent sampling time.

5. Conclusions

In this paper, we presented an approach of determination of hydrogeological parameters of coastal aquifer based on tidally induced and barometrically induced groundwater fluctuations when applied to the River Neretva aquifer located in the southeastern part of the Adriatic Sea. To our best knowledge, this is the first paper demonstrating the applicability of tidal methods for hydrogeological characterization of Neretva Valley coastal aquifer. Besides the application of a wide range of well-established methods, we firstly enabled them to be used for polychromatic tidal signals and then systematically combined them, starting from raw observed signals and determination of barometric pressure influence to the final assessment of the aquifer’s hydrogeological parameters: hydraulic diffusivity, hydraulic conductivity, specific storage, storativity, leakance, clay layer extending roof length, and outlet capping.

Below, we summarize the highlights of the study:

• The conceptual model presented by Equation (13) of [26] can be successfully applied to the River Neretva Valley aquifer to capture transient and periodic features of piezometric heads observed inland from the shore, thus increasing the potential of the method and its capability.
• The solution of the nature of the transient groundwater fluctuations can be determined even in cases of the presence of multiple constituents, thus enabling the use of polychromatic signals in real cases.
• Even when the observed signal was composed of long-term trends, barometric pressure, and tidally induced sea level oscillations, as was the case in this study, by detrending procedure and using the Rahi method [41] a signal can be eliminated until only tidal components exist, thus enabling the use of tidal methods to assess hydrogeological parameters.
• Both spectral component analysis and leakance values were in agreement, confirming aquifer confinement assumed due to the detection of a compact clay layer with, on average, 18–20 m thickness, detected by previous in situ investigation work [2,61,62,63,64,65,66,67].
• Geological structure of the aquifer can be defined from top to bottom as stratigraphically divided in several layers: fine sand with, on average, 8–10 m thickness; a compact clay layer that is, on average, 18–20 m thick; and a gravel layer with a variable thickness ranging from 1 to 130 m until the bedrock. The clay layer extends under the sea, towards the offshore area, with an overall detected length of 1400 m.
• The determined hydraulic conductivity and specific storage did not converge significantly in between different piezometric locations. At D1, specific storage varied between 2.87 × 10⁻⁶ and 4.78 × 10⁻⁶ (m⁻¹) for the porosity range 0.15–0.25. For the same porosity interval, D2 was characterized by 2.98 × 10⁻⁶ to 4.97 × 10⁻⁶ (m⁻¹), whereas the specific storage value of D4 was found between 2.97 × 10⁻⁶ and 4.96 × 10⁻⁶ (m⁻¹). In the same manner, hydraulic conductivity was identified, ranging from 1.4 × 10⁻³ to 2.3 × 10⁻³ (m s⁻¹) for D1, 4.5 × 10⁻³ to 7.5 × 10⁻³ (m s⁻¹) for D2, and 7.0 × 10⁻⁴ to 1.2 × 10⁻³ (m s⁻¹) for D4, indicating horizontal layer geometry.

Findings arising from this paper offer valuable input for further analysis of sea water intrusion dynamics within the River Neretva Valley and enable the detection of the effectiveness of related mitigation techniques in the near future.