Empirical relationships of groundwater head–salinity response to variations of sea level and vertical recharge in coastal confined karst aquifers

Coastal aquifers are characterized by a mixing zone with freshwater–saltwater interactions, which have a strong relationship with hydrological forcings such as astronomical and storm tides, aquifer recharge and pumping effects. These forcings govern the aquifer hydraulic head, the spatial distribution of groundwater salinity and the saline interface position. This work is an empirical evaluation through time-series analysis between aquifer head and groundwater salinity associated with the sea-level dynamics and the aquifer recharge. Groundwater pressure, temperature and salinity were measured in a confined aquifer in the northwest coast of Yucatan (México) during May 2017–May 2018, along with precipitation. Cross-correlation and linear Pearson correlation (r) analyses were performed with the data time series, separating astronomical and meteorological tides and vertical recharge effects. The results show that the astronomical and meteorological tides are directly correlated with the aquifer head response (0.71 < r < 0.99). Salinity has a direct and strong relationship with the astronomical tide (0.76 < r < 0.98), while the meteorological tide does not (r < 0.5). The vertical recharge showed a moderate correlation with the aquifer head (0.5 < r < 0.7) and a nonsignificant correlation with the groundwater salinity (r < 0.5). In this study, the sea level (r > 0.7) is a more important forcing than the vertical recharge (with 0.5 < r < 0.7). Empirical relationships through time-series analysis and the separation of individual hydrological forcings in the analysis are powerful tools to study, define and validate the conceptual model of the aquifer.


Introduction
Coastal aquifers, characterized by freshwater-saltwater interactions, are crucial for maintaining regional biodiversity in various ecosystems, supporting socio-economic development, and providing a reliable source of water supply in coastal zones. These systems are vulnerable to natural and anthropogenic hazards, both chronic stressors (e.g., sea-level rise, aquifer pumping) and acute shocks (e.g., storm surges, extraordinary runoffs, tsunamis), and processes that may cause pollution, saltwater intrusion and freshwater depletion (Werner et al. 2013;Ketabchi et al. 2016). Overall, sea-level position, regional aquifer discharge/abstraction and recharge are the main natural processes and forcing agents that ultimately control coastal aquifer dynamics (Levanon et al. 2017;Dessu et al. 2018).
The saline interface in coastal aquifers has been studied using a wide variety of techniques such as geochemical and geophysical methods, hydrodynamic techniques, simple-multivariate statistical evaluation, geo-spatial Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10040-020-02151-9) contains supplementary material, which is available to authorized users. analysis, and environmental tracers, among others. Each approach has advantages and disadvantages due to their applicability, conditions of use, and uncertainty in the results (Werner et al. 2013;Kumar et al. 2015;Post and Werner 2017;Rachid et al. 2017;Bakker and Schaars 2019;Fernández-Ayuso et al. 2019).
Statistical approaches through time series analysis are effective methods to study the interaction between hydrological forcings and their associated aquifer response. The approaches include empirical and analytical approaches (Bakker and Schaars 2019). Time series analysis has been applied in aquifer head data to study the change in the regime of the aquifer, to obtain hydraulic characteristics such as storage and hydraulic conductivity, and to understand multiple interactions between hydrological forcings and the water table (Trglavcnik et al. 2018;Bakker and Schaars 2019;Fernández-Ayuso et al. 2019). These studies have focused on the effects of hydrologic forcings on the aquifer head, but not on the groundwater salinity or on the dynamic of the saline interface. This gap is explored in the research reported in this paper.
Other studies explore the qualitative and quantitative relationships of groundwater level and aquifer salinity due to temporal variations in hydrological forcings (Sun and Koch 2001;Tularam and Keeler 2006;Kovacs et al. 2017;Levanon et al. 2017;Fiorillo et al. 2018), but these approaches do not analyse each forcing separately. Dessu et al. (2018) study the interaction of hydrological forcings separately with multivariate regression analysis, but the work does not separate astronomical and meteorological tides in marine forcings. Trglavcnik et al. (2018) report that this separation is important to understand the interaction between offshore forcings and coastal aquifers, especially in heterogeneous cases, but their study was not implemented for groundwater salinity. It means that hydrological forcings in coastal aquifers must be separated (especially the tide components) and later analysed for better conceptualization. This separation has not been explored for study of the interaction between forcings and groundwater salinity.
In the Yucatan Peninsula, México, similar studies have been carried out using fractures, sinkholes, and/or offshore submarine springs for monitoring (Vera et al. 2012;Parra et al. 2016;Coutino et al. 2017;Kovacs et al. 2017Kovacs et al. , 2018Young et al. 2018). However, these methodologies are valid for studies on the turbulent mixing within a cave-related domain (e.g., free water-column mixing) in the absence of porous or fractured/karstic-rock matrix . In contrast, studies on the N-NW Yucatán peninsula aquifer show evidence of confinement in the coastal zone (Perry et al. 1989), suggesting that the aquifer behaviour here may differ from those presented in the previous studies cited in the preceding. In fact, Levanon et al. (2017) suggest that the pressure wave propagates from the seafloor boundary, faster than the fluctuation of the freshwater/saline-water interface, and it is controlled by heterogeneities in both the saturated and unsaturated permeability and the specific storage coefficients. This scenario is different compared with the turbulent mixing in a free water column observed in karstic manifestations as sinkholes or offshore submarine springs.
Besides, Coutino et al. (2017) and Kovacs et al. (2017) found a positive correlation between extreme rainfall events (caused by hurricanes), water salinity and groundwater level measured in sinkholes in a karstified environment. Results suggest that sea-level position has a linear correlation with aquifer head variation, while storm surges, wind conditions, droughts and rainy periods are other important stochastic factors that control these phenomena. With respect to offshore spring discharges, the principal forcing agents are tide and the hydraulic gradient between the groundwater-level and sealevel positions (Vera et al. 2012;Parra et al. 2016;Young et al. 2018).
In summary, previous studies offer only an explanation for the relationship between hydrological forcings (precipitation and astronomical and meteorological tides) and the aquifer head, and they do not focus on the salinity. In Yucatan, the studies were only for karstic manifestations (sinkholes and offshore springs) where interactions are different than in the matrix/porous media with confined aquifer conditions. This research analyses the specific influence of each hydrological forcing in the aquifer head and salinity, using linear and crosscorrelation methods, with the aim to assess the empirical relationships between the aquifer head and the groundwater salinity associated with sea level, astronomical and meteorological tidal effects, and vertical recharge.
Thus, this study builds and improves upon methodological precedents by: (1) obtaining a quantitative analysis of the hydrological forcings in the aquifer using reproducible techniques that can be applied to other coastal aquifers in the world with similar conditions, (2) assessing the effects of each hydrological forcing on the aquifer, (3) understanding the salinity mixing mechanism in the aquifer using empirical evaluation through time series analysis and linear and crosscorrelation relationships between hydrological forcings and the aquifer response, and (4) evaluating the differences in the salinity behaviour measured in previous works on karstic sinkholes and in groundwater wells. This approach was tested in the northwest coast of the Yucatan Peninsula, which is a karstic confined aquifer with the presence of a shallow water table and mixing in the fresh-saline water interface.

Study zone
The study zone is located in the northwest coast of Yucatan, Mexico, in a transect between the towns of Hunucma and Sisal, located between 20.80°and 21.20°N latitude and 89.80°and 90.20°W longitude (Fig. 1). The study area has a humid tropical climate with an average annual rainfall of 600 mm and average potential evaporation of 1,600 mm/year; the rainy season is from May to November. The average annual temperature is 25.6°C (CONAGUA 2012). The surface geology ( Fig. 1) shows the existence of Quaternary deposits at the coast (Holocene and recent), and Tertiary deposits consisting of the Carrillo Puerto Formation (Pliocene-Miocene). Regional surficial structures are located in the socalled cenote (a local name referring to karstic sinkholes) ring zone, while the study area is located 10 km to the east of this fractured region (Butterlin and Bonet 1963;Mexican Geological Service 2005;Villasuso and Méndez 2000).
Within the study area, the regional system is an unconfined karstic aquifer, except near to the coastline, where a thin superficial impermeable layer of limestone, locally known as caliche, confines the regional aquifer. Overall, groundwater flows from SE to NW (to the coast). The aquifer has a freshwater lens that overlies a saltwater body whose thickness decreases towards the coast; the saltwater-freshwater mixing zone defines a saline interface (Perry et al. 1989;Villasuso-Pino et al. 2011). The vertical stratigraphy defines two aquifers: (1) towards the coast a local unconfined aquifer located in the dune-based deposits, and (2) the regional aquifer located in karstified limestones; these aquifers are separated by the caliche confining layer. Towards the continent, there is a low-permeability layer of sandy-limestone allowing the formation of another local aquifer (Villasuso-Pino et al. 2011). This conceptual model was validated with physical evidence collected from well drilling samples and video camera surveys, which also show a low degree of heterogeneity, compared to other highly karstified zones in Yucatan Peninsula (e. g., Cenote Ring and Holbox Fault). Figure 2 shows a detailed stratigraphy associated with the geology in the transect P7b-P7a wells (Fig. 1).

Field measurements
Continuous measurements of groundwater pressure, temperature and salinity were collected in the study zone during May 2017-May 2018. Pressure transducer loggers (Onset HOBO, model U-20L-04) were installed in six monitoring wells, to monitor hydrogeological parameters every 30 min, considering a submergence of 1.5 m from the aquifer static water table (Fig. 1). Atmospheric pressure was collected with a BARO-DIVER installed in the P9 well, for barometric compensation. . Section A-A′ represents a transversal geological section between P7a and P7b wells. Black arrows show the regional groundwater flow path in the study zone Electrical conductivity loggers (Onset HOBO, model U-24-002-C) were installed in the P7a (coastal zone) and the P8 (transition inland zone) wells, considering two fixed depths: upper saline interface (USI) and bottom saline interface (BSI) as shown in Fig. 2. The USI was defined using a fixed salinity halocline greater than 1.5 parts per trillion (ppt; first meter in halocline), while the BSI was defined by the halocline zone nearest to saltwater concentration (35 ppt).
Monthly electrical conductivity and temperature profiles in monitoring wells were run using a YSI EXO 2 multiparametric device. These measurements were used to calibrate the electrical conductivity loggers. The salinity (Sal) was calculated using practical salinity scale criteria as reported in Lewis (1980). Precipitation and atmospheric pressure data were measured at meteorological stations installed at Sisal (SSL), P9 well and Merida (MID) (Fig. 1). Data were collected at the 1-min interval and downsampled to 30-min time intervals. Precipitation data were used in the analysis of vertical recharge in the aquifer. Atmospheric pressure was used to audit the BARO-DIVER field measurements and complement data loss due to maintenance. Tide measurements were collected from a tide gauge installed at Sisal, as part of the National Mareographic Network (UNAM 2000). The tidal gauge data were also downsampled to 30-min time intervals.

Forcing analysis
This study analyses the influence of sea level and precipitation, as forcings, on the head, salinity and discharge of the aquifer. This influence was analyzed using the following techniques: (1) tidal harmonic using T_Tide (Pawlowicz et al. 2002) and the power spectrum (Rao and Swamy 2018) for the identification of the main tidal components and the separation of the astronomical and meteorological signals of the tide, (2) a sample cross-correlation (Box et al. 1994, Eq. 1) to calculate the time-lag between forcings and aquifer response, and (3) linear correlation evaluated with the Pearson correlation coefficient described by Eq. (2) (Fisher 1958;Box et al. 1994). These methodologies are described in the following.
where, y1 t and y2 t are time series with sample means y1 and y2, and lags k = 0, ± 1, ± 2…. For data pair (y1 1 , y1 2 ), (y2 1 , Fig. 2 Geological cross-section A-A′ of the study zone (see Fig. 1). Vertical stratigraphy has been validated with drilling wells. Blue and red dots represent the location of pressure and electrical conductivity loggers, respectively. The saline interface is located between the upper saline interface (USI) zone and the bottom saline interface (BSI) zone y2 2 ), …, (y1 t , y2 t ), n is the number of data points in the time series. Finally, Cy1y2 is the covariance of y1 and y2, and sy1 and sy2 are the sample standard deviations of the time series.

Analysis of the sea-level position and its influence on the aquifer
Tides are composed by astronomical and residual components. Astronomical tides are the sum of several tidal constituents, in particular seven fundamental constituents with periods between 11.97 to 25.8 h, and other tidal constituents with higher and lower frequencies (Masselink and Hughes 2003). Long-period tidal constituents have not been analysed in this research because the time series is only a year of data. Residual tides are divided into radiational tides and storm surge. The former is the product of atmospheric variations with fluctuation periods between 12 to 24 h. Storm surges (also called meteorological tides) are a product of noncyclic phenomena such as wind and storm events (Masselink and Hughes 2003). The residual tide in this research is a product of only storm surges because atmospheric variations were compensated using a barometric logger. The principal tidal components were obtained through tidal harmonic analysis using T_Tide (Pawlowicz et al. 2002) and the power spectrum (Rao and Swamy 2018). The difference between the tidal measurements (T, full signal) and the astronomical tide (AT) defines the meteorological tide (MT), as MT = T − AT. Similarly, the well pressure and salinity time series were processed through harmonic analysis as well, to separate the astronomical and meteorological signals in the aquifer head (AH and MH) and the aquifer salinity (AS and MS).
The main tidal components were adjusted using the analytical model of tide propagation in a confined aquifer (Ferris 1952), in which the behaviour is modelled by means of an exponential decay of the sinusoidal tidal wave as where h is the variation in the aquifer's hydraulic head (measured in m), x (m) is the distance, t (s) is time, H o is the tidal amplitude (or sum of tide components), τ is the tidal period, and D (m 2 /s) is the hydraulic diffusivity of the aquifer, which is the ratio of the aquifer transmissivity and storativity. The amplitudes in Eq.
(3) were adjusted using the harmonic analysis for the main tidal components. The goodness of fit was evaluated using the root mean square error (RMSE) calculated between the sum of the tidal components obtained from the harmonic analysis and the sum of the components fitted by the Ferris (1952) model. The meteorological head (MH) was filtered to remove the effects induced by local recharge. For this purpose, MH signals in the aquifer were derived ( Δh = Δt ) and data for which Δh = Δt presented values one or two orders of magnitude larger than the gradient associated with events without a local recharge were removed.
The time lag between AT, MT and the aquifer response (AH, MH, AS and MS) were performed using a sample cross-correlation (Box et al. 1994). After adjusting the lag time between signals in the aquifer and AT/MT (i.e., the astronomical aquifer response f(x i + lag ) is a function of the astronomical tide x i ), a linear regression was adjusted for sea-level signals (AT and MT) and the aquifer response (AH, MH, AS and MS). Similar to other studies and hydrologic models (Legates and McCabe 1999;Kovacs et al. 2017;Dessu et al. 2018), the Pearson correlation coefficient (r) was used as an indicator of the relationship between hydrologic variables. In this study r > 0.70 was accepted as a significant correlation, 0.50 ≤ r ≤ 0.70 as a moderate correlation, and r < 0.5 as a nonsignificant correlation.

Analysis of the aquifer discharge response to forcings
The flux rate is defined by Darcy's law for porous media as Q = K × A × i; while the cubic law defines the flux rate for a fractured matrix as Q = 2K f × A × i; where Q is the flux rate (m 3 /s), K is the hydraulic conductivity (m/s) in porous media, K f is the hydraulic conductivity (m/s) in the fracture, A is the cross-section area (m 2 ), and i is the hydraulic gradient (Bear 1972;Domenico and Schwartz 1997;Cook 2003). This means that the discharge Q is a direct function of the hydraulic gradient i, assuming constant geometry and aquifer characteristics; thus, the hydraulic gradient i can be considered a qualitative indicator of the aquifer discharge variation. Therefore, the hydraulic gradient was computed using a time series of MH and MT. Eqn (4) defines the hydraulic gradient i between two points as the ratio of the difference between the piezometric head in a well (h well ) and the sea level elevation corrected to the equivalent freshwater head (h sea ) and the distance between both points (d). Young et al. (2018) proposed a similar analysis using a hydraulic index, that expresses the pressure gradient between inland springs and submarine springs at the coast, and their results suggest that the discharge flux are strongly related to the pressure gradient.
Precipitation and its influence on the aquifer A total of three rainfall recharge events were isolated to analyse the short-term influence in the aquifer. These events were selected due to precipitation data availability. The time lag between precipitation events and the aquifer response (water level and salinity) and relationship significance were obtained using a sample cross-correlation analysis (Box et al. 1994).
The seasonal influence of rainfall in the aquifer was analyzed using the monthly average of the hydraulic head (M-MH), salinity (Sal), and the monthly accumulated precipitation (MP), values previously corrected to account for meteorological tides effects. The variables MP, Sal and M-MH were plotted in a scatter plot, adjusted with the time lag (the meteorological aquifer response f(x i + lag ) is a function of the monthly accumulated precipitation x i ). A linear regression was adjusted for MP vs M-MH and MP vs Sal. Pearson correlation coefficient (r) was used as an indicator of significant correlation.

Results and discussion
Monthly salinity profiles in wells show a shallower and thinner saline interface towards the coast, until there is an abrupt change (close to 1 m thickness). The freshwater (salinity <1 ppt) thickness is 25 m in the continental aquifer (in well P4), and reduces towards the coast where there is brackish groundwater with a 2 ppt concentration (Fig. 3a). Salinity profiles show a stable fresh/brackish water thickness without discontinuities or stratification. Figure 3b shows that in rainy months (May-November), the aquitard stores freshwater derived from rainfall, attenuating the recharge in the regional aquifer. The salinity in the aquitard decreases by the freshwater storage.
Time series of the aquifer head, salinity and hydraulic gradient in the study area are shown in Fig. 4. Hurricanes were not registered during the study period, although the signal of northern cold fronts can be seen from October 2017 to March 2018. Figure 4a shows the water surface astronomical variation in the sea level and the aquifer, confirming that the astronomical tide propagates significantly into the latter. Figure 4b is a plot of the meteorological head (MH), composed by (1) the effect due to the aquifer recharge from precipitation and (2) the nonastronomical tide propagation into the aquifer. Precipitation from rainy months (June-October, Fig. 4d) recharges the aquifer, as can be seen from the MH signal (1) during specific short-events in early September and late October, and (2) for the rainy period during which MH slowly increases from May to October (Fig. 4b). This aquifer response is clear at the wells farther from the coast (P4, P7b, P8 and P9), but not at P5 and P7a, the wells closer to the coast, where the ocean influence predominates and the aquifer is confined. Zavala-Hidalgo et al. (2010) report that the monthly mean sea level was the lowest in July and highest in October, which is consistent with the monthly mean MH from well P5 and P7a, as shown in Fig. 4b.
The rainy season commonly is related to increments in the hydraulic gradient and an increment in the groundwater flux  (Fig. 4c). These increments in the discharge flux rate related to precipitation can be appreciated in Fig. 4c,d. In fact, given the monthly average sea-level variation mentioned in the preceding, the hydraulic gradient in wells is lower during the higher sea-level period (September-November, Fig. 4c), but also coincides with the largest freshwater recharge period. From mid-December to the end of April the hydraulic gradient in the wells is consistently high (0.5 m/km in P7a), which may be due to a low monthlyaveraged sea level, even if the precipitation is low in that period. Finally, the effects of the cold fronts from December to April, associated with strong northerly winds which produce 20-30-cm storm surges lasting 1-4 days (Torres-Freyermuth et al. 2017), can be seen in the MH signal (Fig.   4b), in particular at the stations closer to the coast and a lesser extent in station P4 and P7b.
The salinity records (Fig. 4e,f) at the upper and bottom saline interface show a variable diurnal response related to the astronomical tide, except at BSI in P8 well, where the ocean influence is minimal. At P7a-BSI the signal decreases during some periods, coincident with the behaviour of the monthly salinity profiles in P7a shown in Fig. 3c.
Throughout the study period, salinity at P7a BSI shows a seasonal dynamic behaviour, ranging from a maximum of 27.7 ppt on July 10th to 13.5 ppt on April 12th (Fig. 4f); this suggests that the BSI moves seasonally up and down, as registered in the electricalconductivity logger located in a fixed position where the diurnal variability is small. The high values in June are associated with (1) the end of the dry season, (2) a low monthly-averaged sea level and (3) a relative increase in the hydraulic gradient. In contrast, the low salinity values occur in three periods: (1) a decrease in October (lowest value of 13.5 ppt on October 20, 2017), which coincides with both the period of the year with highest meteorological head and the largest freshwater input to the aquifer, (2) gradual quasi-monotonic decreases of salinity from July to September 2017 and (3) from February to April 2018, possibly associated with the decrease of the sea level and apparently an increase in freshwater storage, i.e., thickening and deepening of the freshwater layer near the coast.

Sea level forcing
The tidal forcing, arising from planetary motions influenced by the rotation of the earth and the orbits of the moon around the earth and the earth around the sun, can be defined by a set of spectral lines, i.e., the sum of a finite set of monochromatic sinusoids with specific frequencies and phases (tidal constituents, Pawlowicz et al. 2002). Power spectra (PS) of the head data in the monitoring wells correspond to the PS of the astronomical tide (Fig. 5a). The peaks' frequencies, in cycles per hour (cph), coincide with the main lunar diurnal constituent O1 (0.0387 cph), the diurnal lunisolar K1 (0.0417 cph), and main lunar semidiurnal M2 (0.0805 cph). These values are coincident with values reported in local studies (Vera et al. 2012;Torres-Mota et al. 2014) and express the most important forcing wave components of the astronomical tide in the study zone.

Freshwater aquifer head
The tidal influence is observed in wells with distances up tõ 12 km from the coast. The most distant well from the coast (P4 well) did not register the evidence of tidal influence; in the case of P7b well, harmonic analysis suggests a highly attenuated O1 and K1 signal, but the range of aquifer head variation (8 mm) is equivalent to the pressure logger accuracy; thus, this result must be taken with caution. The power spectra in these wells are four orders of magnitude smaller than sea level, suggesting that tidal effects are negligible (Fig. 5a).
The analytical model by Ferris (1952) reveals that the hydraulic diffusivity ranges from 20.5 m 2 /s (P7a) to 2,450 m 2 /s (P9; see Fig. 6), adjusting with the astronomical tidal constituents O1, K1, M2; the percent of AT propagated into the aquifer (linear correlation slope between AT and AH) ranges from 1.4 to 79%- Table 1, Fig. S1 of the electronic supplementary material (ESM)-with an expected Pearson correlation coefficient (r) from 0.77 to 0.99. This model indicates that the astronomical tide has a direct relationship with the elevation of the aquifer inland, associated with the hydraulic diffusivity. The highest correlation values are located in wells near to the coast (P7a and P5 wells), while correlation decreases inland (P8, P9 and P7b wells). The time lag between AT and Fig. 5 Power spectra for the a sea-level elevation and aquifer head (Patm is atmospheric pressure), and b salinity at the wells in the study area Fig. 6 Fit of the astronomical tide propagation in the aquifer using the Ferris model for O1, K1 and M2 tidal components. AH is the hydraulic head for astronomical tide effects and D is the hydraulic diffusivity AH ranged between 0 (P7a well) and 7.0 h (P7b well, see Table 1). Cross-correlation for meteorological effects (MH and MT) shows a direct positive correlation between the meteorological tides and the aquifer head in P7a, P5, P8 and P9 wells; crosscorrelation in P4 and P7b wells does not show significant correlation. The time lag between MT and MH was estimated between 0 and 3.0 h (Table 1), similar to the astronomical tide propagation time. MT and MH data are plotted in Fig. 7. The coastal zone (P7a, P5, and P9 wells) shows higher Pearson correlation (r between 0.82 and 0.71) and a moderate correlation for P8; correlation is nonsignificant in wells farther from the coast (P7b and P4, r < 0.40), suggesting that other variables, such as precipitation or other aquifer oscillations, are related to the aquifer head in this zone.
Meteorological effects propagate into the aquifer with MT percent (linear regression slope) between 40% (P8) and 61% (P7a; see Table 1), increasing near the coastline. Linear regression results indicate that the sea level rise by meteorological effects increases the aquifer head in the continent, being more noticeable near the coast (Fig. 7).
Regarding the tide propagation in coastal aquifers, previous studies have reported a direct relationship between the sea level and the aquifer head (Vera et al. 2012;Levanon et al. 2017;Dessu et al. 2018), but these studies have not considered the individual influences of both the astronomical and meteorological signals. In terms of the decay of sea level propagation into the aquifer, Fig. 8 shows that the astronomical signal (AT) experiences a faster decay than the meteorological component (MT), suggesting that the nonastronomical signal of the ocean tide, which typically has a longer period than the astronomical components analysed here (Torres-Freyermuth et al. 2017) can, theoretically, travel farther inland, as reported by Trglavcnik et al. (2018). The differences between the propagation of the forcings in the aquifer can only be observed when the astronomical and meteorological signals are analysed separately.

Aquifer salinity
Power spectra (PS) in salinity present similar features as the PS in the astronomical tide: the spectra reveal the same frequencies at the main tidal components: O1 (0.0387 cph), K1 (0.0417 cph) and M2 (0.0805 cph). PS in salinity shows a similar order of magnitude to the PS in sea-level and aquifer response. This behaviour suggests that a harmonic analysis can also be an option to separate astronomical and meteorological signals in terms of aquifer salinity. P8-BSI spectra do not show similitude with the main astronomical components in the sea, suggesting that the salinity variations in the aquifer associated with the astronomical tide are strongly attenuated in the fixed logger position in this specific well (Fig. 5b).
The correlation for astronomical effects (AS and AT) shows a strong relationship between the astronomical tide and the astronomical signal in the aquifer salinity. Time lags are between −2.5 and 7 h; the time-lag in P7a BSI is abnormal because the salinity peak occurs before the high tide (−2.5 h), while in stations P7a USI and P8 USI, the time lag is estimated between 3 (P7a) and 7 h (P8). These distortions have been reported in boreholes within the coastal zones (Levanon et al. 2017), as an effect of the absence of a capillary zone. Levanon et al. 2017 suggest that burying the sensors in the aquifer is an option for filtering these distortions. In this study, sensors were installed in boreholes only, and these effects may be present in the results as observed in the PS of P7a-BSI (Fig.  5b). The astronomical signals' responses in the study zone are similar to those in the conceptual model of Levanon et al. (2017) which was defined for unconfined aquifers, suggesting that this conceptualization can be applied in confined systems.
AT and AS are plotted in Fig. 9, showing a direct strong correlation of the astronomical signal for salinity with the astronomical tide (r between 0.76 and 0.95, Table 1). These correlations were expected because the PS of the sea level is coincident with the PS of the salinity in the monitoring wells (Fig. 5b). Results suggest that during a rising tide, the salinity in the aquifer groundwater increases, and inversely decreases during ebb. Furthermore, Fig. 9 shows that (1) the diurnal salinity variation at P7a is only 0.1 ppt at USI, and 0.8 ppt at the BSI; (2) the correlation at P7a is larger at BSI than at USI, and (3) the tidal effect in salinity occurs 5.5 h before at BSI compared to USI. These observations suggest that the astronomical tide effect in the salinity propagates in the vertical direction from the bottom to the upper saline interface (BSI to USI) and then towards the aquifer inland, as previously described by Levanon et al. (2017).
Previous studies report a direct response of the aquifer salinity to the meteorological tides; but the monitoring was implemented only in sinkholes and submarine springs (Vera et al. 2012;Parra et al. 2016), where the salinity mixing mechanism is turbulent in a free water column in the vicinities of these structures, as described by Coutino et al. (2017). In this research, the results do not show a significant relationship between MT and MS, with r < 0.50 (Table 1 and Fig. S2 of the ESM); these differences can be associated with the effect of assessing the hydrological forcings separately (especially the tide components), and the effect of different mixing mechanisms on groundwater salinity compared to highly karstified systems. In this respect, there are no previous studies reporting a similar methodology and/or similar results. Besides, the studies previously mentioned (Vera et al. 2012;Parra et al. 2016;Coutino et al. 2017;Kovacs et al. 2017) were carried out with extreme climatological events (e.g. hurricanes), which were not present in this investigation; therefore, it is not known whether extreme climatological events could generate these differences.

Precipitation forcing
Aquifer head Figure 10 shows three selected rainfall events to compare precipitation, aquifer head and salinity (events in June, July and September 2017). Meteorological stations at MID, P9 and SSL yielded data to show the rainfall distribution throughout the study period.
Cross-correlation analysis was evaluated only for MID and P9 (transition and continent) because the conceptual model suggests a coastal confined aquifer (SSL) with poor response from precipitation; therefore, P7a and P5 were not considered in the analysis. In the case of salinity correlations, only P7a USI and P8 USI were correlated with precipitation in P9, because the recharge effects could only be identified in the upper limit of the saline interface.
Results show that short-term precipitation events at MID and the aquifer head in the study zone present a time lag between 26 and 26.5 h (moderate correlation in September 2017 with 0.50 < r < 0.7, see Table S1 of the ESM), while in station P9 the time lag ranged from 0.5 to 3.0 h, and are poorly correlated (r ≤ 0.5; see Table S1 of the ESM). This suggests that the local precipitation is not a significant forcing in terms of the vertical recharge response of the aquifer, probably due to the low permeability of the vadose zone and the aquitard layer (Fig. 3b).
The seasonal precipitation response in the aquifer was evaluated using monthly accumulated precipitation, which reflects the seasonal trend in the aquifer head due to precipitation and represents seasonal changes over the analysed period. Figure 11 shows that during rainy months (June-November 2017), the aquifer head increases by recharge, and decreases in dry months. Besides, precipitation data at P9 and MID suggest a 3-months delay in the aquifer response. Using this time lag, monthly accumulated precipitation and aquifer level (previously corrected by MT influence, M-MH) were correlated by means of a linear Pearson correlation model as shown in Table 2. Table 2 reveals that the seasonal precipitation presents a moderate direct relationship with the mean aquifer head, similar to previous studies Kovacs et al. 2017). The values of r for P7b, P8 and P9 range between 0.50 and 0.56 with respect to precipitation at MID, and between 0.68 and 0.70 with respect to precipitation at P9. The aquifer head in P4 does not show a correlation with the precipitation data (See Fig. S3 of the ESM). These outcomes suggest that the vertical recharge is not as strongly correlated with the aquifer head as the sea level.

Aquifer salinity
Cross-correlation and Pearson correlation coefficients for precipitation and salinity are nonsignificant (r < 0.50), suggesting  . Previous studies carried out in sinkholes show an inverse relationship between precipitation and salinity due to the direct connection of the aquifer and the surface terrain Kovacs et al. 2017). In this research, the monthly salinity profiles in wells reveal that the vadose zone stores the water of precipitation during rainy months (Fig. 3), thus the saline interface does not show response to this stress. Rocha et al. (2015) report similar results, suggesting that recharge in the aquifer is not related to the variation in the saline interface.

Summary and conclusions
In this study, empirical relationships through time-series analysis and linear and cross-correlation methods were used as a useful tool to understand the dynamics of coastal aquifers and to assess the aquifer response due to regional hydrologic forcings.
In the study area, the astronomical tide is the most correlated (Pearson r-value) hydrological forcing to the aquifer head, followed by the meteorological tides, and the precipitation (vertical recharge). Therefore, the hydraulic gradient and the regional discharge is more sensitive to the sea-level elevation than the vertical recharge. The most correlated forcing associated with the groundwater salinity is the astronomical tide, while meteorological tides and precipitation show a nonsignificant correlation.
High-frequency salinity variations in the study zone are correlated with the diurnal tide, while the meteorological tides are stochastic local waves and do not have a significant relationship with the salinity, suggesting that the groundwater salinity is more sensitive to astronomical waves than stochastic waves, similar to the tide pressure propagation.
In the study region, the vadose zone and the aquitard layer are important factors to consider. The effects of precipitation are attenuated by the confining layer, meaning that the vertical recharge does not affect the saline interface position, as opposed to the highly karstified aquifer in which turbulent salinity mixing occurs in the free water surface (e.g., sinkholes) during rainfall. In fact, the salinity mixing mechanisms in confined porous and/or poorly karstified aquifers are better explained by the propagation of the marine wave head pressure (astronomical and meteorological) towards the aquifer.
This study shows that hydrological forcings in karstic coastal aquifers can be better understood when they, as well as their effects, are assessed separately; for instance, this separation allows one to identify (1) the asymmetry between the propagation of the meteorological and astronomical tides towards the continent, (2) the fact that the aquifer salinity does not show a clear relationship with the meteorological tides as reported in other highly karstified systems, and (3) the relative correlation of each hydrological forcing in terms of the aquifer head/salinity response.
The results of this research offer an opportunity to identify the most correlated hydrological forcings and simplify the variables in the aquifer conceptualization, which is important, for instance, when a numerical model is implemented. Therefore, empirical relationships through time series analysis can be considered as a first target in the study of the aquifer and for the development of conceptual models, and later, to simplify the problem for numerical modelling purposes.
It is important to implement this analysis in the future with longer-term series data using a multiple variable regression approach, which can explain better the interaction of coupled hydrological forcings and their relative importance. A limitation of this research is that empirical relationships are valid only in the portion of aquifer studied, and explains only local effects. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.