Temporal Variations in Chemical Proprieties of Waterbodies within Coastal Polders: Forecast Modeling for Optimizing Water Management Decisions

: In polder-type land, water dynamics are heavily inﬂuenced by the artiﬁcial maintenance of water levels. Polders are low-lying areas of land that have been reclaimed from the sea or from freshwater bodies and are protected from ﬂooding by dikes or other types of ﬂood-protection structures. The water regime in polders is typically managed using a system of canals, pumps, and sluices to control the ﬂow of water in and out of the area. In this study, the temporal changes in water salinity in the polder-type agricultural ﬂoodplain within the Neretva River Delta (NRD), Croatia, were analyzed by applying multivariate statistics and forecast modelling. The main aim of the study was to test the model that can be used in practice to forecast, primarily, water suitability for irrigation in a coastal low-lying agricultural catchment. The speciﬁc aim of this study was to use hydrochemistry data series to explain processes in water salinity dynamics and to test the model which may provide accurate salinity prediction, or ﬁnally select the conditions in which the model can be applied. We considered the accuracy of the model, and it was validated using independent data sets. To describe different patterns of chemical changes in different water classes due to their complex hydrological connectivity, multivariate statistics (PCA) were coupled with time-series analysis and Vector Autoregression (VAR) model forecasting. The multivariate statistics applied here did not indicate a clear connection between water salinity of the surface-water bodies and groundwater. The lack of correlation lies in the complex hydrological dynamics and interconnectivity of the water bodies highly affected by the artiﬁcial maintenance of the groundwater level within the polder area, as well as interventions in the temporal release of freshwater into the drainage canal network. Not all individual water classes contributed equally to the dominant patterns of ionic species identiﬁed by PCA. Apparently, land use and agricultural management practices in the different polders lead to uneven water chemistry and the predominant contributions of speciﬁc ions, especially nutrients. After applying the Granger causality test to reveal the causal information and explain hidden relationships among the variables, only two surface-water and two groundwater monitoring locations displayed a strong causal relationship between water electrical conductivity (EC w ) as an effect and sea level as a possible cause. The developed models can be used to evaluate and emphasize the unique characteristics and phenomena of low-lying land and to communicate their importance and inﬂuence to management authorities and agricultural producers in managing and planning irrigation management in the wider Mediterranean area.


Introduction
The fertile soils found in delta regions globally are critical for food production, and the loss of this land could have serious implications for global food security. The increase in sea level and the intrusion of saltwater into aquifers is already affecting the quality of the soil and reducing its productivity [1,2]. In addition, climate change alters the marine environment, primarily due to the rise in sea level. According to the IPCC (Intergovernmental Panel on Climate Change), the mean global rate of sea-level rise between 1901 and 2010 was 1.7 mm per year, with a higher rate of 3.2 mm per year between 1993 and 2010 and between 1920 and 1950 [3]. Climate change is also having a significant impact on the study area, particularly over the Adriatic Sea, leading to changes in water dynamics and salinity. Studies of long-term sea level  from a range of tide gauges in the Adriatic indicate sea-level rise between 2 and 3.4 mm per year when atmospheric and steric effects are considered, with an error of 1 mm per year that could be due to land subsidence [4]. Similar results [5] were obtained for both absolute (altimetry) and relative (tide gauges) sea-level changes from 1993 to 2008, both showing positive trends, with values of 3.2 mm per year and 1.9 mm per year, respectively. A recent review on historic and future trends in sea-level rise in the North Adriatic Sea, specifically Venice, gave an estimate of 1.23 ± 0.13 mm/year, based on an analysis of the observational period from 1872 to 2019 [6]. The same study emphasized Atlantic hydrographic boundary conditions as a major source of uncertainty for future projection of Mediterranean Sea levels. Land subsidence in the low-lying coastal areas, such as the Neretva River delta (NRD) and similar locations on both sides of the Adriatic Sea, makes the estimation of local sea-level variations even more difficult. Additionally, the continuing relative sea-level rise (RSLR) is already affecting groundwater in terms of freshwater quality and quantity due to the saline wedge into the coastal aquifers, as reported by Scardino et al. [7] for the Southern Adriatic.
Other challenges related to climate change are a decrease in water availability and quality in Eastern Mediterranean countries [8][9][10][11], an increase in soil aridity and salinity, land desertification, and exhaustion of water resources [12][13][14], which ultimately affect crop yield [15]. Along the Eastern Adriatic coast, particularly in its southern part, the shortage of water, severe droughts, and anomalous weather events in recent decades have become the main constraints on efficient agricultural productivity [16,17]. As in other seas and coastal areas, seawater intrusion (SWI) is a growing problem exacerbated by sea-level rise, especially in river deltas, influenced by both natural and anthropogenic causes [18]. Natural triggers generally include downstream impacts, such as sea-level change, tidal fluctuations, waves [19], and land subsidence [20,21], as well as impacts from the upstream catchment, such as change in freshwater inflow. In lowland coastal areas, even a low rate of sea-level rise causes an inland shift of saline groundwater. This problem has been addressed for many European low-lying coastal areas, both on a large scale [22,23] and for small watersheds [24,25]. Water issues, such as overallocation, environmental flows, and river salinity, are all influenced by the degree and nature of the connectivity between rivers and aquifers [26]. The temporal hydrological connectivity of different water bodies induced by the polder operating system enhances not only fluvial dynamics but also water chemistry. The NRD has been converted to agricultural land through numerous and extensive reclamation projects that began in the late 1950s and resulted in the creation of polder-type agricultural land. Past river floodplain hydroengineering measures reduced the side branches and meanders of the NRD from twelve to just three [15]. Generally, the surface-water levels within polder-type land are artificially maintained and controlled by a network of drainage canals, pumping stations, and sluice-gate complexes, which are used to release water into the sea [16,27]. In addition to natural impacts, freshwater inflow into the lower part of the Neretva River is severely affected by the operation of hydropower plants upstream in Bosnia and Herzegovina, especially during dry periods.
The actions to integrate spatial and temporal data on water-quality issues into sustainable water management require the consideration of additional information, primarily on the timing, intensity, and duration of water quality in specific environmental assets.
Long-term monitoring programs commonly produce abundant time-series data that allow, among other things, the evaluation of water-quality trends [28]. However, when dealing with a highly modified deltaic floodplain, designing and implementing an efficient monitoring program is a highly complex task [29]. As defined by Alfonso et al. [30], polders drain excess water independently into higher areas or to the sea, forming a system that requires observation of water levels in every polder so that proper management of the hydraulic structures that connect them is achieved. Optimization of the monitoring network in such an aspect is performed to include all water bodies that have different functions within the polder. The segmentation of the polder-type land by active land use with agricultural dikes, ditches, and roads strongly influences the dynamics of water circulation in the floodplains [31]. The rapid migration of saltwater into fresh aquifers by surface and subsurface flow paths over short time scales currently receives less attention than lateral seawater intrusion and long-term vertical SWI [32]. Accurate short-term predictions on water salinity dynamics can play a critical role in irrigation management for farmers. Understanding the changes in water salinity levels can help farmers make informed decisions on the timing, amount, and type of irrigation to use, which can ultimately impact crop yield and quality [33]. In the NRD, the Croatian Waters authority manages and operates the polder system to ensure feasible management and planning. As part of these activities, a soil and water monitoring program has been conducted regularly since 2009 [16], providing a comprehensive database of surface and groundwater-quality parameters related to salinization processes and intensive agricultural production, including soil-quality, land-use, and meteorological data, and adequate inputs for geostatistical modelling.
Many previous studies examining the temporal variations in water properties in coastal polders, or their details, highlight a practical application in improving flood forecasting, water-quality management, or sustainable land-use planning to implement a climate adaptation strategy. Karrasch et al. [34], in their attempt to develop an "actor-based" scenario in the low-lying coastal areas of Northwest Germany, proposed the application of freshwater retention areas (polders) to prevent potential uncontrolled flooding of the hinterland. In such cases, one should provide decision makers with clearly drafted and spatially explicit scenarios. The specific aspects of such empirical and practical research depend on the characteristics of the specific resource unit, the spatio-temporal scale, and/or the social-ecological setting. In the context of our study, the shortage of freshwater available for irrigation in a small-scale coastal area occurs seasonally due to climate change and seawater intrusion. Zuurbier et al. [35] studied the effects of aquifer storage and recovery (ASR) on freshwater management in coastal areas of Nootdorp, the Netherlands, using wells. In this approach, freshwater surpluses are stored for later use in times of demand, creating a self-sufficient freshwater supply which makes external freshwater supply (including infrastructure) and/or costly desalination superfluous. The final success of the application and operational optimization of ASR systems depends on both geochemical reactions and processes in the aquifer and the management strategies to combat potential water-quality issues.
Yu et al. [36] used the data on 25 different variables on 144 polders to investigate impact of groundwater on surface-water quality and nutrient loads in the Amsterdam area, the Netherlands. Polders with high chloride concentrations were located along the coastline (Zuiderzee margin) and in the Upconing area, which is known to experience upconing of brackish water. The seepage of brackish water is driven by the low surfacewater levels in addition to drainage via pumping stations. Similar results regarding chloride concentrations were observed in surface waters. The authors found similar spatial patterns in solute concentrations between groundwater and surface waters, with much higher concentrations in groundwater. Analyzing geological structures and different hydrological and water-quality variables in a deep Duch polder, de Louw et al. [37] found that groundwater seepage is not a uniform process and that different types of seepage differ in regard to flux, chloride concentration, and the location within the polder. Therefore, their contributions to surface-water salinization also differ.
In areas where both water quantity and quality are influenced by a complex array of natural and anthropogenic factors, multiple lines of variables and statistical methods should be applied to understand and manage the processes [38]. Various geostatistical models are being applied as well to investigate surface-water connectivity dynamics [39,40]. Romić et al. [16] tested two linear mixed-effects models for predicting electrical conductivity and nitrate concentration in surface waters and groundwater, which were found to be efficient in adequately reflecting the heterogeneity and complexity of the NRD. Principal Component Analysis (PCA), a statistical tool of feature extraction, in which fewer independent variables are created from a combination of the original parameters, has been applied to several environmental issues, including water quality and water-quality monitoring [41,42]. This process facilitates the visualization of obtained data and understanding of the interrelationship of various variables, conditions, and locations. Additionally, the use of predictive models that consider the artificial maintenance of the water regime is useful in polder-type land as it could help in managing the water level and ensure that the land is used efficiently and sustainably.
Thus, to improve the understanding of the complex interactions of ecological, hydrological, and water-quality factors, several numerical or mathematical models can successfully be applied. The potential of coupling specific statistical tools, such as PCA and Granger causality (GC) analysis, has been tested in the studies of the relationships among water-quality indicators and other environmental variables. By exploring how PCA and GC analysis can explain relationships between water quality, hydroclimatic variables, and watershed characteristics, Zavareh et al. [43] found a large variation within watersheds, despite their similarity. The Granger causality test describes the dynamic correlation among variables and does not indicate true causal relationships [44], but the recognized patterns present in the water-quality data enable an effective forecast. Classic prediction models employed in a wide array of water-quality studies include the Autoregressive Moving Average (ARMA) and the Box-Jenkins Autoregressive Integrated Moving Average (ARIMA) [45]. Vector Autoregression (VAR) is the other statistical model used in time-series forecasting. The main difference between ARIMA and VAR is their scope and complexity. ARIMA is a univariate model, which means that it considers only one time series at a time. ARIMA models the relationships between one observation and a set of lagged observations and residual errors, making it suitable for modelling stationary time-series data. VAR, on the other hand, is a multivariate model that considers multiple time series at the same time. VAR models the relationships between multiple variables, making it more suitable for modeling non-stationary time-series data and data sets with multiple variables; in research by Kaur et al. [46], a VAR model was compared to a Vector Error Correction Model (VECM) and a more complex machine-learning Long Short-Term Memory (LSTM) model for determining the changes in water levels and the water flow of different water bodies across Italy. The results showed the significance of time-series models (VAR and VECM) over LSTM, where the VAR model gave reliable results for water bodies, such as aquifers, rivers, and lakes, while the best results for springs were obtained by the VECM model. Keng et al. [47] coupled Granger causality with VAR in forecasting groundwater levels in Malaysia. A VAR model with eight lags was the most suitable model for representing the relationship between groundwater level and rainfall. While this study showed that VAR can be used in predicting groundwater level in the study area, other variables should have been considered, since VAR allows multiple variables to be included in modeling. Similarly, Ramli et al. [48] used monthly data of rainfall and discharge (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) to predict rainfall for the next five years using a VAR model in Peghasin District, Indonesia. The predicted rainfall value had an R 2 value of 0.66. The results indicated that from the available data set (2008-2015) reliable prediction can be performed for only two years (2016-2017). As in the research of [47], other variables influencing rainfall may have been considered.
In the context of water-quality prediction, a VAR model was found to be more appropriate if multiple factors affecting water quality were considered, such as salinity by SWI and nutrient inputs from nearby agricultural areas. For instance, Veerendra et al. [49] used machine-learning techniques that included ARIMA and VAR to predict the results of water tests for a few water components and found that the VAR model was more reliable and that its values were nearer to the original values than those of the ARIMA model. In a poldertype flood plain, waterbodies are hydrologically connected in various dimensions. Such interactions between waterbodies result in chemical heterogeneity, including heterogeneity in salt and nutrient contents [50].
Rising sea levels will obviously aggravate the issue of SWI so that a finely determined spatial and temporal scale is required to capture the essential dynamics of water salinity associated with mixed processes and various conditions in the case of the NRD. Accurate forecasting of water salinity offers useful insights for near-future polder management and crop production, including management options for water and salinization caused by SWI polder operation and irrigation practices, which was an expected contribution of this study.
For this purpose, multivariate statistics were combined with time-series analysis and VAR model forecasting to propose an indicator-based monitoring and assessment program. In this study, we used monthly water-quality data for the period 2010 to 2020 to develop a model for short-term forecasting of water salinity. The general objective was to use data on water quality in the NRD to explain processes related to water salinity dynamics and to test the model which may provide accurate salinity prediction. The specific aims were to (i) identify common features in the chemical composition of various classes of water bodies by applying multivariate statistics, (ii) establish the relationship between different classes of surface waters and groundwater that governs the magnitude and extent of water salinity, (iii) test the stationarity of time series and investigate the causal relationship between sea level and water salinity, and (iv) recognize patterns in water salinity by developing a 12-month-ahead salinity forecast using the VAR model.

Study Area Description
The research area is in the NRD, on the southeastern coast of the Adriatic Sea (43 • 00 N, 17 • 30 E) ( Figure 1). The Neretva River is the largest river of the Adriatic catchment area and, according to the annual water discharge, is one of the largest rivers of the Mediterranean Sea [51]. The Neretva River, with its numerous tributaries in its lower course forms an extensive low-lying delta, which influences the main river either directly at the surface or indirectly in the groundwater. The entire delta is under the strong influence of the Adriatic Sea, both through SWI through the main riverbed and groundwater. In the delta, the Neretva River overlies Quaternary deposits. Carbonate rocks, mainly limestone, form the edge of the valley, its flanks, and smaller isolated hill hummocks within the valley [52,53]. These rocks are intensely fractured and deeply karstified. Surface Quaternary sediments consist mainly of peat and clay (organic marshy deposits), which are underlain by clayey sands, sands of variable grain size, sandy clay, gravelly sands, sandy gravel, Holocene gravels (alluvial sediments), and Pleistocene conglomerates. Geomorphology, soil characteristics, and reclamation in the study area are described in detail by Romić et al. [54] and Kralj et al. [55]. Elevations above sea level of individual polders were determined using a digital elevation model (DEM) of the study area. Since the elevation of the polders within the ameliorated area is lower than in the surrounding area, continuous drainage of excess water is required through a dense network of drainage canals, gates, and pumping stations [16]. Polder-like agricultural land is intensively used for the cultivation of citrus fruits and vegetables. The characterization of the Quaternary aquifer was made by a group of authors [25,56] who compiled the results of a series of hydrogeological, geotechnical, and geophysical studies carried out in the area during the last 60 years ( Figure S1). The aquifer system is composed of three main layers. The uppermost layer consists mainly of fine sand with a locally present shallow sublayer of clay and silt. This layer is present throughout the NRD area and represents the unconfined aquifer with a total thickness of 1-10 m. Beneath the unconfined layer a compact clayey-silt layer with a thickness of 10-25 m, which rises slightly towards the sea and forms an aquitard, has been determined. This layer continues below the seabed offshore. Finally, beneath the clay layer is a confined aquifer of mostly coarse gravel with a total depth of 130 m unless bedrock is encountered. The bedrock in this area is karstified limestone, whose depth varies from the ground surface at the edge of the valley to a maximum depth of 165 m below ground level in the middle of the valley [56].
The area is semiarid with a Mediterranean climate, with dry, hot summers and humid, mild winters. The mean annual rainfall for the main weather station of the area (Ploče) is 1077 mm (1988-2020), with most of the rainfall occurring in the period from October to April. The minimum rainfall occurs in July (28 mm) and the maximum occurs in November (153 mm). The mean annual air temperature is 15.9 °C, with January being the coldest (7.0 °C) and July the warmest month (25.7 °C). The annual Penman-Montheith reference evapotranspiration is 1217 mm, with the highest value of 197 mm occurring in July.

Monitoring Network Design and Data Collection
In the study area, the surface-water-and groundwater-quality monitoring program, established in 2009, has been focused on salinization to identify long-term changes that have both natural and anthropogenic causes and to provide farmers and decision makers with effective data from which to derive information for management and restoration. The monitoring network covers an area of 5815 ha and includes four different polders with surface and groundwater monitoring locations. In addition, surface waters were The characterization of the Quaternary aquifer was made by a group of authors [25,56] who compiled the results of a series of hydrogeological, geotechnical, and geophysical studies carried out in the area during the last 60 years ( Figure S1). The aquifer system is composed of three main layers. The uppermost layer consists mainly of fine sand with a locally present shallow sublayer of clay and silt. This layer is present throughout the NRD area and represents the unconfined aquifer with a total thickness of 1-10 m. Beneath the unconfined layer a compact clayey-silt layer with a thickness of 10-25 m, which rises slightly towards the sea and forms an aquitard, has been determined. This layer continues below the seabed offshore. Finally, beneath the clay layer is a confined aquifer of mostly coarse gravel with a total depth of 130 m unless bedrock is encountered. The bedrock in this area is karstified limestone, whose depth varies from the ground surface at the edge of the valley to a maximum depth of 165 m below ground level in the middle of the valley [56].
The area is semiarid with a Mediterranean climate, with dry, hot summers and humid, mild winters. The mean annual rainfall for the main weather station of the area (Ploče) is 1077 mm (1988-2020), with most of the rainfall occurring in the period from October to April. The minimum rainfall occurs in July (28 mm) and the maximum occurs in November (153 mm). The mean annual air temperature is 15.9 • C, with January being the coldest (7.0 • C) and July the warmest month (25.7 • C). The annual Penman-Montheith reference evapotranspiration is 1217 mm, with the highest value of 197 mm occurring in July.

Monitoring Network Design and Data Collection
In the study area, the surface-water-and groundwater-quality monitoring program, established in 2009, has been focused on salinization to identify long-term changes that have both natural and anthropogenic causes and to provide farmers and decision makers with effective data from which to derive information for management and restoration. The monitoring network covers an area of 5815 ha and includes four different polders with surface and groundwater monitoring locations. In addition, surface waters were grouped into classes based on their common characteristics and functions within the polders. Shallow piezometers made of 110 mm diameter plastic pipes with a porous section of 3 to 4 m in length for groundwater monitoring were installed in the central part of each polder. Water samples were collected at several locations within each water class (see Table S1). Water-quality data from the 1 January 2010 to the 31 December 2020 from 15 surface-water bodies (three classes: Class 1 RIVER/LATERAL CANAL, Class 2 PUMPING STATION, and Class 3 DRAINAGE CANAL) and 7 groundwater locations (Class 4 PIEZOMETER) were processed in this study. The water-quality data from 2021 were used to validate the developed forecasting models. The geographical setting of the study, including four polders within the floodplain and the water sampling locations, is shown in Figure 1.
Water samples were collected monthly by direct filling with a hand-held bottle sampler (surface water) and by pipe sampling from piezometers specifically designed for groundwater sampling. Sampling, sample transport, and handling followed a standardized procedure [57,58]. The samples used for the determination of orthophosphate (o-PO 4 ) were filtered through a membrane filter (0.45 µm mesh size) immediately after sampling. Samples used for the determination of ammonium nitrogen (NH 4 -N) were preserved by adding sulfuric acid to adjust the sample to a pH of approximately 2. Water electrical conductivity (EC w ) was measured using a Mettler Toledo MPC conductivity meter with a cell constant control and temperature compensation device (25 • C), according to [59]. pH was determined with a Schott Lab 870 pH meter [60]. Analysis of Cl − and SO 4 2− , NO 3 -N, NH 4 -N, and o-PO 4 was performed with a segmented flow system (SFA) and spectrometric detection using a Skalar San+ Analyzer following standard procedures [61][62][63][64], respectively. Determination of sodium and potassium was performed using an AAS PerkinElmer 3110 atomic emission spectrometer following standard procedures [65]. Determination of Ca 2+ and Mg 2+ was performed with a Varian VISTA-MPX CCD Simultaneous ICP-OES [66]. The determination of bicarbonates was performed by acid-base titration with H 2 SO 4 . Data quality-control assurance for water analysis was implemented using a quality system accredited to [67] standard, by using field blank samples, field duplicates, participation in the international proficiency testing program, and internal reference samples.
Monthly mean sea levels were recorded at the tide gauge station located in the Adriatic Sea at the Mala Neretva River mouth, and the data were provided by Croatian Waters (the public water management company in Croatia).

Statistics
The univariate summary statistics of the variables that describe water chemistry were calculated first. A total of twelve variables were selected (pH, EC w , Ca 2+ , Mg 2+ , K + , Na + , HCO 3 − , Cl − , SO 4 2− , o-PO 4 , NO 3 -N, and NH 4 -N). In addition to summary statistics, oneway analysis of variance (ANOVA) was performed for all variables between the different water classes along with Tukey's multiple comparison test. Univariate statistical analysis and ANOVA were performed using the statistical software package XLSTAT [68].

Principal Component Analysis
PCA was applied to determine correlations among the variables and to summarize the variation within the data. PCA is a technique for reducing the dimensionality of large data sets, increasing interpretability but at the same time minimizing information loss [69]. PCA involves a mathematical procedure that reduces the variables into major controlling components, thereby enabling the interpretation of variable distributions (common sources of ions and hydrochemical processes). Here, PCA was performed in the R programming language using the FactoMineR package [70].

Time Series and Testing for Stationarity
A time series is a sequence of observations taken sequentially in time [71] and there are several possible objectives in analyzing time series, such as description, explanation, control, and prediction [72]. Stationarity is an important concept in time-series analysis that has a large impact on how data are perceived and predicted. When forecasting or predicting future events, most time-series models assume that each data point is independent of the others (stationary). A time series is said to be strictly stationary if its properties (or the process that generates them) are unaffected by changes in its temporal origin [73]. There are several tests for testing the stationarity properties of time series, and one of the most commonly used is the unit root test proposed by Dickey and Fuller [74]. The null hypothesis of the Dickey-Fuller (DF) test is that there is a unit root in the first-order Autoregressive (AR) model, implying that the data series is not stationary [75]. In the case of higher-order autocorrelation, the Augmented Dickey-Fuller (ADF) test is used to detect the higher-order autocorrelation. If the calculated p-value is below the set significance level, the null hypothesis is rejected, indicating that the series is stationary. If non-stationarity is present in the data, the data must be transformed, often by differencing, i.e., by applying the difference operator to the original time series to obtain a new one [73]. In first order differencing, the difference between consecutive data points is calculated [76] (Equation (1)).

Granger Causality Test
The GC test is one of the most robust methods available to detect the causal interactions between two time series and was initially introduced in the field of economics [77,78]. The GC test is a statistical hypothesis test for determining whether one time series is useful for forecasting another. If the probability value is less than any level, then the hypothesis would be rejected at that level. GC addresses the issue with prediction rather than correlation to identify causation between time-series variables [79]. One of the main assumptions of the GC test is the stationarity of the time series. To conduct the analysis, Granger [78] proposed a bivariate model between two time series (Equations (2) and (3)), where X and Y are two stationary time series; a, b, c and d are coefficients; and ε t η t are two uncorrelated white-noise series.

Vector Autoregressive Model
A Univariate Autoregressive (AR) model is a linear model with a single equation and a single variable, in which the current value of a variable is explained by its lagged (previous) values [80]. In the VAR model, on the other hand, there are multiple variables, and each of them is modelled as a linear combination of the lagged values of itself and the lagged values of the other variables in the system. It represents an extension of the univariate AR model to dynamic multivariate time series and was proposed by Sims [81]. For a set of K time-series variables y t = (y 1t , . . . , y Kt )', a VAR model captures their dynamic interaction. The basic p-lag VAR model is of the following form (Equation (4)).
where A i 's (K × K) are coefficient matrices and ε t is an unobservable error term [82]. One of the advantages of the VAR model is that it does not require as much knowledge about a variable affecting the factors but rather a list of variables that are assumed to interact intertemporally [83]. The lag length for the VAR(p) model may be determined using model selection criteria. The general approach is to fit VAR(p) models with orders p = 0, ..., pmax and to choose the value of p that minimizes some model selection criteria. The three most common information criteria are the Akaike (AIC), Schwarz-Bayesian (BIC), and Hannan-Quinn (HQ) criteria [84].

Model Performance Assessment
Many different criteria, such as the coefficient of determination (R 2 ), mean square error (MSE), sum of square errors (SSE), root mean square error (RMSE), mean absolute error (MAE), and many others, are used by authors to evaluate the performance of developed statistical models [85]. The coefficient of determination (R 2 ) is a statistical tool that helps in identifying the relationship between variables in a regression model and represents a measure of the proportion of explained variance present in the data [86]. Low values of R 2 (<0.5) indicate a weak relation between the predictor variables and the response variables. The calculation of R 2 is as follows: The root mean square error (RMSE) is used to calculate the standard deviation of the residuals between the observed and predicted data values [87]. Whereas R 2 is a relative measure of fit, RMSE is an absolute measure of fit, and lower values indicate better fit [88]. It is calculated as follows: Mean absolute error (MAE) represents the average of the absolute difference between observed and predicted values and is often used to measure the prediction error in timeseries analysis [87]. It is calculated as follows: Sum of square error (SSE) represents the difference between the observed and predicted value (Equation (8)) while the mean square error (MSE) represents the average of the squared difference between the observed and predicted data (Equation (9)).
PCA, ADF tests, GC tests, and VAR models were performed using R statistical software [89].

Surface-Water and Groundwater Chemistry
In this study, three classes of surface-water bodies were defined in terms of their common characteristics and using two categories of indicators: (1) natural state or/and extent of modification and (2) hydrological functions of the water body within the poldertype agricultural catchment. Piezometers representing the groundwater body are installed in each polder in the NRD to (1) assess the impact of agricultural practices on groundwater and (2) assess the impact of saline groundwater on the salinization of soil and consequently on agricultural crops. Descriptive statistics on the chemistry of the water bodies divided into classes, as shown in Table 1, clearly indicate different environmental influences.
Different water types can be distinguished based on their ion composition.     Table 1.
Water bodies belonging to the Class 1 RIVER/LATERAL CANAL (R/LC) were supposed to reflect natural conditions and were mostly unaffected by anthropogenic pressures, although they were considered pristine reference waterbodies. The Class 1 R/LC water bodies are fed freshwater from the catchment. They are least affected by saline water intrusion and are low in nutrients ( Table 1). The order of the main cation mean concentrations was Na + > Ca 2+ > Mg 2+ > K + , and their concentrations were significantly lower than those in the other water-body classes. Maximum concentrations of all cations and  Table 1.
Water bodies belonging to the Class 1 RIVER/LATERAL CANAL (R/LC) were supposed to reflect natural conditions and were mostly unaffected by anthropogenic pressures, although they were considered pristine reference waterbodies. The Class 1 R/LC water bodies are fed freshwater from the catchment. They are least affected by saline water intrusion and are low in nutrients ( Table 1). The order of the main cation mean concentrations was Na + > Ca 2+ > Mg 2+ > K + , and their concentrations were significantly lower than those in the other water-body classes. Maximum concentrations of all cations and anions except nutrients (N species) indicate the periodical intrusion of the saline water into the river courses, which has been obviously modifying the natural state common for the river water composition, particularly in Class 1 LC Vidrice and R Mala Neretva, which are located closest to the sea (Figure 1). Class 2 and Class 3 water bodies exhibit similar physicochemical characteristics, which are highly affected by hydrological changes due to technical maintenance of the polder-type floodplain. The mean electrical conductivity values (EC w ) at Class 2 PUMPING STATION (PS) and in Class 3 DRAINAGE CANAL (DC) were higher than those of the Class 1 water bodies. Although the mean EC w values of Class 2 PS and Class 3 DC water bodies did not differ substantially, the range of EC w was much wider in the DC, reaching 36 dS m −1 in certain periods of the year. The order of the main cation mean concentrations did not differ compared to all classes in the water of Classes 2 and 3, being Na + > Ca 2+ > Mg 2+ > K + . Mean nutrient concentrations were low, and the highest maximum concentrations of nitrate and ammonium were recorded in drainage canals.
Class 4 PIEZOMETER (Pi) exhibited distinctively different hydrochemistry compared to the classes of the surface-water bodies. As a result of periodic intense SWI into coastal aquifers, the mean EC w value was the highest (10.7 dS m −1 ) and significantly different from the surface-water classes. The range of EC w values in groundwater was similar to the range for the same parameter in Class 3 DC, but the mean was significantly higher (11 dS m −1 ). Class 3 DC water bodies showed similar EC w value ranges, but with distinctively lower mean values, which indicated hydrological connectivity between surface water and groundwater. The mean concentrations of the main cations were significantly higher in Class 4 Pi than in the surface-water classes, but the order of the main cation mean concentrations was still the same: Na + > Ca 2+ > Mg 2+ > K + . The maximum EC w value was 38.9 dS m −1 , which was close to the seawater salinity level. Obviously, all water bodies in the study area have been impacted by increasing EC w .
ANOVA showed that pH differed significantly among all classes of water bodies ( Table 1). The main cation concentrations in Class 1 R/LC water bodies differed significantly from those in Class 2 PS and Class 3 DC, which did not differ. The nutrient concentrations were not significantly different in the surface-water bodies but were significantly different from the nutrient concentrations in the groundwater. The significant difference in ion composition between all surface-water bodies and groundwater was confirmed by the analysis of variance (Table 1).

Multivariate Statistics
As shown in Figures S2-S4, PCA was performed on hydrochemical variables for all classes of surface water as well as for every class of surface water separately. Most of the variation for all classes is explained by the first two principal components, which cumulatively explained more than 70% of the variation in surface-water bodies, and more than 76% in groundwater ( Table 2). Even though the water bodies were grouped into classes based on their common hydrological features, not all of them contributed equally to the dominant patterns of ion species identified by PCA (Figure 3).  As shown in Figure 3A, which presents a PCA biplot for surface-water bodies, the highest contribution to water salinity belongs to the DC Luke, located in the polder that receives a large amount of groundwater seepage, while the drainage canals and pumping stations collect and transfer most of the nutrients leached from the agricultural fields. The PCA variable contributions for each surface-water class are given separately in Table S2, and related PCA biplots are shown in Figures S2-S4. Visual comparison of the biplots and analysis of clusters of the same observations in different classes of surface-water bodies shows distinct patterns for the contribution of ions, which are confirmed by the variable contribution indices shown in Table S2.
In the case of Class 1 R/LC, the locations closest to the sea, which are R Mala Neretva and LC Vidrice, contribute most to the water salinity, explaining 46.5% of the variance by The PCA results for all surface-water bodies analyzed jointly indicate that the most relevant variables in all water classes (i.e., the ones that explain most of the variance) are EC w , Na + , Ca 2+ , Cl − , Mg 2+ , K + , and SO 4 2− , which are related to the SWI. In this context, water salinization may be assigned as a dominant process, explaining 57.6% of the total variance (Table 2). Nutrient loadings to the surface waters (NO 3 -N, NH 4 -N, and o-PO 4 ) were evident as well ( Figure 3A), and they explain 12.5% of the variance in the second component, reflecting direct anthropogenic pressure deriving from agricultural activities. Nevertheless, due to differences in agricultural management practices in the polders, not all surface-water bodies contributed equally to both principal components.
As shown in Figure 3A, which presents a PCA biplot for surface-water bodies, the highest contribution to water salinity belongs to the DC Luke, located in the polder that receives a large amount of groundwater seepage, while the drainage canals and pumping stations collect and transfer most of the nutrients leached from the agricultural fields. The PCA variable contributions for each surface-water class are given separately in Table S2, and related PCA biplots are shown in Figures S2-S4. Visual comparison of the biplots and analysis of clusters of the same observations in different classes of surface-water bodies shows distinct patterns for the contribution of ions, which are confirmed by the variable contribution indices shown in Table S2.
In the case of Class 1 R/LC, the locations closest to the sea, which are R Mala Neretva and LC Vidrice, contribute most to the water salinity, explaining 46.5% of the variance by the first component. Major variable contributions in the first PCA dimension come from EC w, Mg 2+ , K + , Na + , Cl − , and SO 4 2− . The dominant process affecting water chemistry identified in the second component is the freshwater inflow through the main river course and lateral canals. Major variable contributions in the second PCA dimension come from pH, EC w , Ca 2+ , K + , and HCO 3 − . The contribution of nutrient loadings is rather low, scattered unevenly within the third, fourth, and fifth dimensions, going in the opposite direction of the two main driving forces ( Figure S2). This phenomenon is related to the intrinsic character of the streams grouped in Class 1 R/LC. Numerous tributaries influence the main river stream, either directly on the surface or indirectly in groundwater, and drain mountains of Bosnia and Hercegovina built of magmatic, metamorphic, carbonate, and clastic sedimentary rocks.
In contrast, surface-water bodies of Class 2 PS and Class 3 DC contribute the most to the nutrient loadings. The contribution of nutrient concentrations in surface-water classes is scattered unevenly among the second to the fifth components (Table S3). However, the share of contribution to the nutrient loadings among the surface-water classes suggests several sources and possibly different intensities of NO 3 -N, NH 4 -N, and o-PO 4 fluxes within the floodplain. The highest loading of nutrients is attributed to the water bodies at the pumping stations. Among Class 3 DC, Luke exhibits the highest salinity for drainage water, while DC Jasenska, DC Vidrice, and DC Vrbovci drain mostly nutrients. Obviously, differences in the agricultural land use of polders and dominant crop production result in unequal water chemistry and the prevailing contributions of specific ions, particularly nutrients. Class 2 PS cluster water bodies are supplied by pumping stations within the polders in the process of maintaining the groundwater table and the water level in the canals. Visual inspection of the biplot presented in Figure S3 shows an even contribution to the nutrient loadings of all monitoring locations belonging to Class 2 PS but a more pronounced impact of PS Luke on the water salinity. For all surface-water classes, PCA recorded a similar contribution of nutrient loadings to the water chemical composition, although the specific variable contributions in the third, fourth, and fifth dimensions indicate different nutrient sources and/or flow paths.
The highest contribution of DC Luke to the salinity of Class 3 DC is evident from the biplot shown in Figure S4. Along with DC Vidrice, the same location highly contributes to the nutrient loadings. The surface elevation of the polder Luke located on the right side of the Neretva River varies between −3.04 and 0.63 m above sea level. Land surface elevation is the highest on the river embankments and fluvial terraces along the main water courses, decreasing toward the former swamp area. In this polder, land subsidence is the most pronounced, which further accelerates groundwater seepage and impacts surface-water quality. Unlike other polders, where saline water intrudes into groundwater through the karstified aquifer, water salinization in Luke derives from the salt wedge in the Neretva River. However, the flow pattern and intensity of salinization of the water in this polder depend on the dynamics of the pumping station operation and obviously increase with the pumping of excess water.
The PCA biplot of Class 4 Pi ( Figure 3B) shows a clear distinction in groundwater chemistry among monitoring locations. The first principal component that explains 61% of the variance is evidently related to the water salinity, most pronounced in Pi Jasenska, as shown in Table S4. This polder was developed by the reclamation of the salt lake, and the location obviously communicates with the sea even at present. The only partial overlapping with the Pi Jasenska was observed for Pi Komin, located at the edge of the former lake.

Augmented Dickey-Fuller Test
An ADF test is performed before developing time-series models, and the p-values for each site and water-quality parameter are given in Table S5. In all surface-water classes (R/LC, PS, and DC), the time series for variables related to seawater, i.e., EC w and all major ions, are stationary, except for Cl − at DC Vidrice (p-value = 0.07) and SO 4 2− at R Neretva  Table S5, are more complex and are unevenly scattered both among locations and among parameters at individual locations. According to the ADF test, only the time series for pH at each location within the Pi water class are stationary. At Pi Luke, the time series for all variables except NH 4 -N are non-stationary, in addition to the low p-values for the other two nutrients, o-PO 4 (0.06) and NO 3 -N (0.09). In addition, the parameters related to salinity are non-stationary, due to the salinization processes deriving from the salt wedge in the main riverbed of the Neretva River, which is seasonally affected by both river inflow and SWI. Only at two sites, Pi Jasenska and Pi Vidrice, are the EC w time series stationary (Table S5). In addition to the EC w time series being stationary, the time series for ions closely related to water salinity are also either stationary at both locations (Mg 2+ , Cl − , and HCO 3 − ) or have p-values very close to alpha of 0.05 (Na + and K + at Pi Jasenska). Polders Vidrice (where Pi Vidrice was installed) and polder Opuzen (where Pi Jasenska was installed) are located closest to the Adriatic Sea and have low mean elevations of −0.44 m a.s.l. and −1.13 m a.s.l., respectively. Since the stationarity of the time series is one of the prerequisites for the development of forecast models, all non-stationary time series were transformed to stationary ones by differencing.

Granger Causality
As we studied the hydrological response of sea-level fluctuation on water salinity within the polders, the GC test was applied to assess whether one time series is useful in forecasting another, i.e., whether one variable at time t-lag causes another variable at time t. Table 3 shows p-values for the GC test when water salinity (EC w ) is considered as the effect and sea levels are considered as possible causes. We presumed that Class 3 DC and Class 2 PS are under the dominant influence of the polder maintenance system (i.e., pumping station operation dynamics), so we excluded these classes from the GC testing. The GC test was performed on groundwater and water bodies of the Class 1 R/LC to identify and describe the hydrological effect of sea-level fluctuation on groundwater and surface-water salinity caused by direct SWI into aquifers. As already mentioned, in polder Luke, water salinization derives from the salt wedge in the River Neretva, and therefore this location was excluded from the GC test. The null hypothesis was that there is no GC between the cause and the effect and that the lower the p-value is, the stronger the causality. The results presented in Table 3 show that there is a strong causal relationship between sea level and water salinity only at two locations of class 1 R/LC (R Mala Neretva and LC Vidrice) and at two locations of Class 4 Pi (Pi Jasenska and Pi Vidrice). R Mala Neretva and Pi Jasenska are located in polder Opuzen, which has the lowest individual point elevation of −4.05 m a.s.l. and an average elevation of −1.13 m a.s.l. The polder Vidrice, where LC Vidrice and Pi Vidrice are located, has a mean elevation of −0.44 m a.s.l., and the lowest point is −2.61 m a.s.l. Moreover, these two polders are closest to the Adriatic Sea, which may explain the causal relationship between the water salinity and the sea level.

Forecasting Based on the Vector Autoregressive Model
The VAR model is a multivariate time-series model used to examine the dynamic relationship between two or more variables that interact with each other, where each variable is modelled as a linear combination of the lagged value of both itself and the lagged values of other variables in the system. VAR models were developed for four locations that showed causal relationships with the GC test, Class 1 R/LC R Mala Neretva and LC Vidirce as surface-water bodies and Class 4 Pi Jasenska and Pi Vidrce as groundwater bodies. Nine endogenous variables associated with SWI were selected: pH, EC w , Ca 2+ , Mg 2+ , K + , Na + , HCO 3 − , Cl − , and SO 4 2− . The maximum number of lags was determined using the AIC. R 2 was used to evaluate how well the developed models fit the observed data.
The results show that the VAR models developed for both Class 1 R/LC locations for pH, EC w , Ca 2+ , and HCO 3 − and for K + for LC Vidrice have R 2 values greater than 0.9, indicating that they fit the observed data well ( Table 4). The models for Mg 2+ , Na + , and Cl − for both locations and for K + for R Mala Neretva have R 2 values between 0.6 and 0.8, and only the model for SO 4 2− at LC Vidrice has a value below 0.6 and explains only 50% of the variability of the fitted data. In general, the VAR models developed for two Class 1 R/LC locations explain more than 60% of the variability in the observed variables for all parameters except SO 4 2− at LC Vidrice and can therefore be used for reliable forecasting. The results of the VAR models for the two Class 4 Pi locations Pi Vidrice and Pi Jasenska differ significantly from those for Class 1 R/LC. None of the developed VAR models for forecasting any of the nine variables have R 2 values higher than 0.53 (Table 4). At Pi Jasenska, all models have a very low R 2 value, and only the models for Ca 2+ and Mg 2+ have R 2 values higher than 0.2. At Pi Vidrice, the VAR models for Mg 2+ and Ca 2+ have the highest R 2 values, 0.50 and 0.53, respectively (Table 4). Consequently, none of the developed VAR models could provide reliable predictions for Class 4 Pi locations Pi Vidrice and Pi Jasenska; therefore, only the surface-water Class 1 locations LC Vidrice and R Mala Neretva were selected for forecasting.
Forecasts for EC w based on the VAR model with monthly granularity for one year (2021) ahead were made. Figure 4A shows the forecast for EC w for 2021 at Class 1 R/LC location R Mala Neretva. The forecasted EC w values for 2021 varied from 0.48 dS m −1 in April to 2.09 dS m −1 in October, with a mean value of 1.36 dS m −1 ( Table 5). The average value of EC w for the eleven-year period (2010-2020) was 1.18 dS m −1 . Figure 4B shows the forecast for EC w for 2021 at the Class 1 R/LC location LC Vidrice. The forecasted minimum EC w for 2021 was 1.39 dS m −1 (January) and the maximum of 2.63 dS m −1 was for May (Table 5). The forecasted annual mean was 1.92 dS m −1 , which corresponds to the eleven-year average at this location (1.75 dS m −1 ). The forecasted values of EC w and confidence intervals are given in Table 5. The results show that it is possible to develop reliable short-term (1-year-ahead) forecasts of water salinity at two Class 1 R/LC locations within the NRD that are very important not only for managing the system but also for planning agricultural production. Indeed, these two locations, the watercourses R Mala Neretva and LC Vidrice, whose main purpose is to protect the area from flooding, also have an important function as major sources of freshwater needed for irrigation of agricultural crops in polders Opuzen (where R Mala Neretva is located) and Vidrice (where LC Vidrice is located), which together cover an area of over 3000 ha, representing more than half of the agricultural area of the NRD.  Figure 4B shows the forecast for ECw for 2021 at the Class 1 R/LC location LC Vidrice. The forecasted minimum ECw for 2021 was 1.39 dS m −1 (January) and the maximum of 2.63 dS m −1 was for May ( Table 5). The forecasted annual mean was 1.92 dS m −1 , which corresponds to the elevenyear average at this location (1.75 dS m −1 ). The forecasted values of ECw and confidence intervals are given in Table 5. The results show that it is possible to develop reliable shortterm (1-year-ahead) forecasts of water salinity at two Class 1 R/LC locations within the NRD that are very important not only for managing the system but also for planning agricultural production. Indeed, these two locations, the watercourses R Mala Neretva and LC Vidrice, whose main purpose is to protect the area from flooding, also have an important function as major sources of freshwater needed for irrigation of agricultural crops in polders Opuzen (where R Mala Neretva is located) and Vidrice (where LC Vidrice is located), which together cover an area of over 3000 ha, representing more than half of the agricultural area of the NRD.    Model validation is an important part of multivariate data analysis. Model validation is used to compare model predictions with the observed and unknown data set to assess model accuracy and predictive ability [90,91]. The observed EC w data from 2021 were plotted with forecasted data for 2021, including confidence intervals at the Class 1 R/LC locations R Mala Neretva ( Figure 5A) and LC Vidrice ( Figure 5B). Visually, all observed and predicted data points at both locations fall between the lower and upper bounds of the confidence interval. To check the goodness of fit of the model, RMSE values were calculated. The calculated RMSE was 0.8 dS m −1 (Table 5)

Discussion
In this study, time series of hydrochemical data of the polder-type agricultural floodplain within the NRD were used to interpret the temporal changes in water salinity and to develop and test the short-term forecasting model needed to assist in optimizing water management decisions. The most specific features of the terrain in the study area comprise the lower-lying parcels of predominantly polder-type land formed in the past by intensive land reclamation and hydro-technical interventions. The complex land reclamation projects in the study area have been intended to ensure flood protection in the rainy season

Discussion
In this study, time series of hydrochemical data of the polder-type agricultural floodplain within the NRD were used to interpret the temporal changes in water salinity and to develop and test the short-term forecasting model needed to assist in optimizing water management decisions. The most specific features of the terrain in the study area comprise the lower-lying parcels of predominantly polder-type land formed in the past by intensive land reclamation and hydro-technical interventions. The complex land reclamation projects in the study area have been intended to ensure flood protection in the rainy season with dikes along the Neretva River. On the Mala Neretva River and at the "mouth" into the sea, sluice gates have been constructed to regulate the flow of high waters and prevent SWI into the local waters of the Mala Neretva during the irrigation season. Maintaining the appropriate regime of flow in drainage canals is directly related to the groundwater level. Low land elevations (partly even below sea level) and shallow groundwater levels require the application of permanent pumping to prevent flooding polders. Such hydrological interventions foster water and solute flow, and hydrologic models provide valuable information for identifying problems and supporting decisions [92]. Due to disparate characteristics, stressor combinations, scales, and management practices, no conceptual methodologies for analyzing linked groundwater-surface-water interactions in a lowland catchment have been developed to date [93]. Groundwater and surface water are not separate components of the hydrological system. Rather, they are linked and interact across a wide range of physiographic and catchment settings, and the exchange of water between streams and shallow groundwater may influence the hydrochemistry in both directions.

Water Salinization Patterns
Human alterations to flow paths, such as ditches, drainage canals, and similar features in the polder-type floodplain, are aimed at increasing drainage efficiency. However, the natural paths of water and solute flow through the floodplain can be modified, as confirmed by the PCA results. Similar to the results of Li et al. [39], which clearly showed surface hydrological connectivity across the floodplain that resulted in an enhanced spatial similarity in the water quality of the seasonal floodplain lakes, we also found such similarity in terms of a range of water-quality parameters, primarily those related to salinity. This arises from the polders' flooding control and water pumping system.
Although the region has been extensively studied, none of the previous studies have attempted to analyze the nature of spatial and temporal variations in hydrochemical characteristics in a connected system of water bodies in an artificially managed poldertype floodplain of the NRD. As an ultimate 280 km 2 segment of a more than 10,000 km 2 catchment of the Neretva River strikingly different from the rest of the catchment in many features, geological and hydrogeological, the Neretva River itself has been modified through the construction of in-stream hydropower plants and reservoirs. The study of Kralj et al. [55] explained the nature of depositional conditions recorded in Neretva River sediment cores, which are certainly influenced by long-lasting human activity in the whole catchment area. Hydroelectric dams and water-storage structures have changed the historical discharge regime and contributed to sediment deficit in the deltaic part of the river, causing changes in river morphology and ecology.
Romić et al. [16] applied two linear mixed-effects models for the prediction of electrical conductivity and nitrate concentration in the surface waters and groundwater in the NRD. The model was proven to be useful to alert for water-salinization risk due to natural processes. In contrast, the prediction of nitrate concentration in the water was not equally satisfactory, apparently stemming from the heterogeneity of the land management and agricultural practices. In our research, we were further focused on the temporal patterns of surface-water and groundwater salinization induced by both natural processes and human activities.
In terms of food security, ensuring a safe water supply for agriculture in terms of both quantity and quality is the foremost objective in integrated water management. The estuary mouth empties into the Adriatic Sea in the Neretva Channel. Water inflow from the wide catchment spread of over 10,000 km 2 reaching this point, along with the small streams and springs on the karst edges, is fresh and constant but seasonally variant. In the summer months, when the freshwater inflow is low, seawater enters the riverbed and forms a seawater wedge that can reach more than 20 km upstream from the river mouth (at the town of Metković on the Croatian-Bosnian border). On the other hand, if the freshwater inflow from the catchment area is higher, the salt wedge is pushed downstream. Changes in river flow control salt-wedge dynamics. The spatiotemporal patterns of the extent of the salt wedge and stratification studied by Ljubenkov and Vranješ [94] showed that saline water reaches the border 20 km upstream at a freshwater inflow of less than 180 m 3 s −1 , while it is completely displaced from the riverbed at an inflow of more than 500 m 3 s −1 . In other words, the Neretva River is highly stratified, i.e., the freshwater flows above the seawater wedge with a very narrow transition zone (approximately 0.5 m). The balance of input from tidally driven seawater and river inflows that sets up stratification patterns becomes highly disturbed by the water pumping within the polders.

Water Hydrochemistry and PCA Analysis
Different water bodies or classes of surface waters and groundwater tend to develop distinct hydrochemical characteristics depending on the dynamic mixing of water from different sources [16]. The SWI intensity and seasonal dynamics as well as the water chemistry of surface-water bodies are related predominantly to the inflows from the river catchment. As a narrow delta-shaped riverine estuary, the area is morphologically complex due to the highly transient patterns of river inputs and hydrological connectivity within the different types of water bodies. Hydrochemical changes in various surface-water bodies, such as different types of drainage canals, are induced by more complex processes activated primarily by the pumping-station operation regime that maintains the water table within the polder area. Exfiltration of shallow saline groundwater to prevent the water table from accumulating water from drainage canals may adversely affect the quality of surface water used for irrigation and threaten crop production [95].
The hydrochemistry of surface-water bodies analyzed together by PCA is dominantly affected by the SWI. However, hydrological processes differ among both different surfacewater classes and distinct locations within each class. Patterns of water salinization in low-elevation land are highly dynamic, with large variations over not only seasonal but also diurnal timescales due to pumping-station operations [96]. According to the classification of Geyer et al. [97], the NRD is a typical salt-wedge estuary due to its microtidal environment and high freshwater flow rates [98]. As a result of the low tidal dynamics of the Adriatic Sea, where mean tidal amplitudes do not exceed 30 cm [99], the Neretva River is highly stratified, which means that freshwater flows above the seawater wedge. Groundwater salinization in the Quaternary aquifer occurs from several sources. In the west, it is in direct contact with the sea, so SWI is clearly visible in both the unconfined and the confined aquifers (very high EC w values recorded). In the north is the Neretva River, which for most of the year contains a wedge of seawater in the lower part of the water column. In the dry summer months, the thickness of the seawater wedge increases. Given that the groundwater level in the unconfined aquifer is kept lower than the water level of the Neretva River [25], by means of pumping stations and drainage canals, saline water from the Neretva infiltrates into the aquifer. According to the available measurements [100], the EC w of the groundwater at the location of the piezometer located south of the Neretva River and near Pi Jasenska generally exceeds 2 dS m −1 . A third source of groundwater salinity is the seepage of saline water from a confined aquifer into an unconfined aquifer, whose piezometric level is higher than the water level in the unconfined aquifer, through an aquitard. However, considering the significantly lower permeability of the aquitard, the share of salinization from this source is significantly lower compared to the intrusion of saline water from the sea and the Neretva River. The groundwater recharge by freshwater takes place through the infiltration of precipitation and inflow from the marginal parts of the valley where the groundwater accumulated in the karst aquifer flows out at the springs.
In the southern part of the research area is the river Mala Neretva, with occasionally increased EC w . Previous research [25] did not show a clear connection between the Mala Neretva River and groundwater. Regarding the groundwater from the piezometer which is located in the immediate vicinity of the Mala Neretva River, in the first 60 cm of the water column, occasionally the EC w reaches less than 1 dS m −1 , which could indicate the existence of a connection with the river and the inflow of freshwater. However, this is not always the case. Sometimes the EC w levels of the Mala Neretva River are lower than in the groundwater. Given that PS Vidrice is located nearby, it can be assumed that the operation regime of this pumping station affects the EC w distribution. This issue will definitely need to be investigated in the future. Therefore, the study area, although it represents a rather small flood plain, is an area of high geomorphological and hydrological complexity, including intensive hydrodynamics caused by the operation of pumping stations.
Polders in the NRD are maintained by constant pumping at pumping stations located at the lowest point at each polder, and excess polder water is discharged to the sea through the surface-water system at higher elevations. Due to the complex geomorphology of the wider catchment area, especially regarding the karst formations surrounding the deltaic floodplain, the hydrochemistry of the various surface-water bodies, including polder open channels, is highly erratic. In the Netherlands, due to sandy and intensively drained low-land catchments, most of the streamflow originates from groundwater [93]. In contrast to the findings of Louw et al. [37] that upwards saline groundwater seepage leads to surface-water salinization of deep lying polders in the Netherlands, in a densely drained lowland catchment of the Neretva River, salt fluxes are affected mainly by polder pumping operations and salt-wedge intrusion into the main river course. Nonetheless, the water monitoring network set up in the NRD proved to be efficient in identifying the hydrochemical heterogeneity of different water bodies. Alfonso et al. [30] addressed the problem of network design and evaluation of monitoring water levels and quality in polders due to a considerable number of different target water bodies that need to be monitored in relatively small, hydraulically disconnected areas. The water management system in the NRD is operationally optimized to control the water level and salinity level in the polders. A key objective in monitoring the water quality of the water courses in the polder network is to control the downstream salinity level in rivers, interconnected channels, and groundwater. A specific aspect of this study was to capture the essential dynamics of water salinity within the polders associated with the seawater level that allow short-term predictions.
Among Class 1 R/LC, R Mala Neretva and LC Vidrice are the locations most influenced by seawater due to their proximity to the Adriatic Sea. For Class 2 PS, most of the saltwater is collected at PS Luke and PS Vidrice. However, different processes are responsible for such a picture. Polder Luke is under the strong influence of the stratified flow of the Neretva River, and the saline water from the river that enters the polder is discharged by continuous pumping towards the monitoring location PS Luke. Although the operation of the pumping station diverts and collects water from the entire polder, the contribution of PS Luke to nutrients is low, which is due to poor land use, many abandoned parcels, and the low intensity of agricultural production in this polder. In contrast to PS Luke, the nutrient contribution of PS Vrbovci was higher, which was due to the intensive use of agricultural land. At PS Vidrice, there was a mixed influence of fresh and saltwater due to its position in the polder. PS Vidrice is located on the left bank of the Mala Neretva River and pumps the collected water from the Vidrice polder into the river. As a result of proximity and the fact that the Mala Neretva River is at a higher elevation than the pumping station, PS Vidrice is influenced by freshwater from the Mala Neretva River (contributions of Ca 2+ and HCO 3 − ). Last, PS Opuzen, although closest to the sea, showed very little to no contribution to salinity, which is due to the dominant freshwater inflow from the entire polder area rather than SWI at the pumping station. The intensive agricultural production of fruits and vegetables in this area contributes to some extent to nutrient loadings.
Unlike PS Luke, DC Luke has a much greater contribution of saltwater but also of nutrients. The contribution to saltwater comes from the Neretva River and the salt wedge, and the higher contribution of DC Luke in nutrients is local in nature, as the DC monitoring location is located on a drainage canal between intensively used agricultural parcels. Other Class 3 DC locations have a minor effect on saltwater transport but have high nutrient contents due to agricultural production. In Class 4 Pi, clear distinctions in groundwater chemistry were found among the monitoring locations. The first principal component, apparently related to water salinity, was most pronounced in Pi Jasenska. The high EC w values observed at Pi Jasenska may be the result of saline water "pockets" somehow related to the remnants of the saline Modrič Lake, which occupied this area before the reclamation projects in the 20th century [15,101]. Remnant saltwater might have been present in the low permeable formations, which now diffuses into the local aquifers. In some coastal areas, ancient seawater might have been trapped in low permeable formations, and as the relative sea level drops and fresh groundwater infiltrates, the saltwater can diffuse into the local aquifers, causing salinization [102]. The understanding of past hydraulic conditions, coastal geology, and relative sea-level changes can serve as a basis for further research and investigations into the salinization of coastal aquifers in the study area.

Time-Series Analysis and Forecasting
To develop and test a model for short-term forecasting of water salinity, monthly data on water quality for the period 1 January 2010-31 December 2020 were used. The model was eventually tested using the 2021 observed data. Primarily, the forecasting of salinity enables management decisions in relation to salinity levels in the irrigation season. The results suggest that SWI is the primary cause of salinization in all polders except the Luke polder. To investigate the causal relationship between sea level measured at the tide gauge and water salinity, a GC test was conducted.
Stationarity as a prerequisite for time-series modelling and forecasting was observed for most parameters in surface waters, in contrast to groundwater, where results were unevenly distributed across locations and parameters. Non-stationarity in hydrological data is common and can be attributed to environmental and/or anthropogenic climate variability and change [103], especially for groundwater time series [104,105]. Although the non-stationary time series were transformed by differencing, the developed VAR models for the selected Class 4 Pi locations Pi Vidrice and Pi Jasenska had low R 2 values for all variables and could not be used for reliable predictions. This could be due to the complex combination of different parameters affecting groundwater in the NRD. Studies by Lovrinović et al. [25] in the NRD have shown that the shallower unconfined aquifer is subject to a more complex influence caused by many factors simultaneously, including SWI, inflow from larger watercourses, recharge from precipitation, and, of course, maintenance of polders (PS regimes and drainage canals) and agricultural production (irrigation).
The developed VAR models for two Class 1 R/LC locations had a high R 2 value for most parameters, particularly EC w , indicating a good fit with the observed data and implying that a reliable short-term forecast (1 year ahead) could be made. Comparison of the forecasted EC w for 2021 with the EC w actually observed in 2021 showed relatively low RMSE values at both locations.

Management Implications of Forecasting and Indicator Development
In the past, there have been several projects in the NRD aimed at improving water management through the use of pumping and drainage systems in polders, which altered the pattern of outflow. The study of Kralj et al. [55] identified reduced sediment deposition in the river channel system, floodplain, and estuary; the decrease in terrain elevation; and the progressive lowering of the delta surface due to reduction in sediment supply. In addition, a number of worrying side effects occurred as well, such as streambank erosion and main-channel deepening in the lower flow due to gradual erosion caused by the let-off of the excess water accumulated in the dams. Additionally, the recent efforts to enable navigation of the river, the construction of the flood-prevention system, and the closing of one of the armlets to reduce SWI have contributed significantly to the changes in sediment fluxes in the delta. Developing reliable short-term forecasting models for EC w is important not only for managing the system, but also for the planning of agricultural production, especially in the polders directly affected by SWI. In addition to flood protection, the monitoring locations R Mala Neretva and LC Vidrice have an important function as main sources of freshwater needed for irrigation of agricultural crops in the Opuzen polder (where R Mala Neretva is located) and the Vidrice polder (where LC Vidrice is located), which together cover an area of more than 3000 ha, representing more than half of the agricultural area of the NRD. To mitigate salinization in coastal polders, such as NRD, forecasting models can provide insights into future salinity levels and help in implementing preventive measures. By considering factors such as tidal dynamics, groundwater intrusion, and irrigation practices, these models can simulate and predict salinity variations within the polder. Based on these predictions, farmers and water managers can implement strategies such as controlled drainage, which involves managing water levels to prevent the rise of saline groundwater and subsequent salt accumulation in soils. Additionally, the models developed can assist in the design and operation of water-control structures, such as sluices and pumps, to regulate water levels and prevent flooding or waterlogging. Specifically for the NRD study region, models can aid in identifying suitable locations for freshwater reservoirs or constructing barriers to prevent seawater intrusion. We should not avoid mentioning the benefits of forecasting models in improving water quality in coastal polders by predicting nutrient concentrations and identifying potential pollution hotspots. The models can assist in identifying critical periods when nutrient loading is high, allowing for the timely implementation of measures to minimize environmental impacts.
In order to accurately predict water salinity in the study area's surface-water monitoring locations, the developed models have proven to be effective. However, it is important to acknowledge that these models do not adequately represent groundwater conditions. Therefore, future studies should consider the inclusion of additional variables to enhance the accuracy of salinity prediction models for various locations and water bodies. In addition to incorporating additional variables, exploring non-linear relationships between variables can also contribute to more accurate predictions. Machine-learning models, which have the capacity to analyze complex and non-linear relationships, offer a promising avenue for future research.

Conclusions
Sea-level rise in the context of climate change is one of the most important factors triggering SWI in Mediterranean coastal aquifers. Other factors were identified as well in the study area's polder-type land that have an impact on the movement of dissolved substances. These factors include flood-protection structures, the artificial maintenance of water levels, and the control of water flow into and out of the area. According to the defined goals, the following conclusions can be drawn:

•
The study revealed a mutual influence between stream water and shallow groundwater in terms of their chemical composition. The results of the PCA performed in the study showed that the exchange of water between these two water bodies is a significant factor affecting their hydrochemistry. However, not all water classes had the same impact on the dominant patterns of ionic species. Differences in land use and agricultural practices in the different polders resulted in uneven water chemistry and a greater influence of certain ions, particularly nutrients.

•
The multivariate statistics applied here did not indicate a clear connection between water salinity of the surface-water bodies and groundwater. The area is influenced by seawater intrusion from the Adriatic Sea, which can occur directly into the aquifer or through the Neretva River and its stratified flow, particularly in combination with low freshwater inflow during summer when irrigation demands are high.
• Granger causality (GC) test results showed a strong causal relationship between electrical conductivity (EC w ) and sea level at only two surface-water and two groundwater monitoring sites. These sites were used to make 1-year-ahead predictions for EC w using a Vector Autoregression (VAR) model.

•
The VAR models were successful in making accurate predictions for the two surfacewater locations but not for the two groundwater locations. Two prediction models were validated using data from 2021 and showed a high degree of accuracy. As a result, the developed VAR model may be useful for authorities in watershed management and agricultural producers in two polders with over 3000 hectares of agricultural land in the NRD to make future production and irrigation plans.
The successful application of the VAR model applied in this study is limited to the two surface-water monitoring locations assessed in this study. Further research and model development are needed to improve the accuracy of predictions for groundwater conditions and expand the applicability of the model to other areas.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agriculture13061162/s1, Figure S1: Hydrogeological map (a) and cross section of the Neretva River valley (b); Table S1: Surface-and groundwater-quality monitoring locations (NRD, Croatia) grouped into four water classes and important characteristics of each class; Figure S2: PCA Biplot Class 1 RIVER/LATERAL CANAL (R/LC); Table S2: Results for the PCA variable contributions in the first five dimensions for each water class; Figure S3: PCA Biplot Class PUMPING STATION (PS); Table S3: PCA results-eigenvalues, variance (%), and cumulative variance (%) of surface-water-body and associated water classes; Figure S4: PCA Biplot Class 3 DRAINAGE CANAL (DC); Table S4: PCA results-eigenvalues, variance (%), and cumulative variance (%) of groundwater body/class piezometer; Table S5: Augmented Dickey-Fuller stationarity test results for all location and water-quality parameters with lag orders and calculated p-values.