*Article* **Participatory and Integrated Modelling under Contentious Water Use in Semiarid Basins**

**Rodrigo Rojas 1,\* , Juan Castilla-Rho 2,† , Gabriella Bennison <sup>3</sup> , Robert Bridgart <sup>4</sup> , Camilo Prats <sup>5</sup> and Edmundo Claro <sup>3</sup>**


**Abstract:** Addressing modern water management challenges requires the integration of physical, environmental and socio-economic aspects, including diverse stakeholders' values, interests and goals. Early stakeholder involvement increases the likelihood of acceptance and legitimacy of potential solutions to these challenges. Participatory modelling allows stakeholders to co-design solutions, thus facilitating knowledge co-construction/social learning. In this work, we combine integrated modelling and participatory modelling to develop and deploy a digital platform supporting decisionmaking for water management in a semiarid basin under contentious water use. The purpose of this tool is exploring "on-the-fly" alternative water management strategies and potential policy pathways with stakeholders. We first co-designed specific water management strategies/impact indicators and collected local knowledge about farmers' behaviour regarding groundwater regulation. Second, we coupled a node–link water balance model, a groundwater model and an agent-based model in a digital platform (SimCopiapo) for scenario exploration. This was done with constant input from key stakeholders through a participatory process. Our results suggest that reductions of groundwater demand (40%) alone are not sufficient to capture stakeholders' interests and steer the system towards sustainable water use, and thus a portfolio of management strategies including exchanges of water rights, improvements to hydraulic infrastructure and robust enforcement policies is required. The establishment of an efficient enforcement policy to monitor compliance on caps imposed on groundwater use and sanction those breaching this regulation is required to trigger the minimum momentum for policy acceptance. Finally, the participatory modelling process led to the definition of a diverse collection of strategies/impact indicators, which are reflections of the stakeholders' interests. This indicates that not only the final product—i.e., SimCopiapo—is of value but also the process leading to its creation.

**Keywords:** stakeholder participation; surface water-groundwater interaction; scenario modelling; integrated water management; agent-based modelling; SimCopiapo

### **1. Introduction**

Water resources are fundamental for supporting livelihoods, food production, energy generation and ecosystem services across the globe. Despite their relevance, water systems are under continuous threats, thus undermining water security [1] and promoting water stress [2]. Interdependencies between water, ecological and social systems across multiple

**Citation:** Rojas, R.; Castilla-Rho, J.; Bennison, G.; Bridgart, R.; Prats, C.; Claro, E. Participatory and Integrated Modelling under Contentious Water Use in Semiarid Basins. *Hydrology* **2022**, *9*, 49. https://doi.org/10.3390/ hydrology9030049

Academic Editors: Il-Moon Chung, Sun Woo Chang, Yeonsang Hwang and Yeonjoo Kim

Received: 14 February 2022 Accepted: 11 March 2022 Published: 16 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

scales and dimensions (e.g., water–energy–food–environment nexus [3,4]) continuously challenge the way water resources have been managed [5]. In this regard, Hoff [6] states that " . . . water management and governance have not yet adapted to these cross-scale and cross-sectoral interdependencies and their dynamics and associated uncertainties".

Water management challenges are no longer addressed solely as technical problems but rather have become part of complex policy and decision-making processes, where multiple stakeholders and institutions reflecting an array of diverse values and interests are involved [7–9]. "Integrated" approaches to account for the array of drivers that help to constrain/condition these water management challenges have therefore received a surge of attention in recent years [10–12].

Kelly et al. [13] discuss the term "integration" in the context of integrated assessment and define five levels with multiple loci in the modelling process: (a) integrated treatment of issues, (b) integration with stakeholders, (c) integration of disciplines, (d) integration of processes and models, and (e) integration of scales of consideration. Integration of biophysical and socioeconomic aspects [10,14] and integration across processes/models (e.g., surface water and groundwater interactions [15,16]), as well as integration with stakeholders [17], have all been documented in the water management-related literature. In the context of surface water–groundwater interactions, Barthel and Barnhaz [16] suggest that "integrated modelling" should explore aspects beyond the purely physical coupling process between surface water and groundwater systems and cover multiple scientific domains and disciplines, thus aligning with the level "integration of processes and models" proposed by Kelly et al. [13].

Jakeman et al. [18] suggest that the development and application of integrated modelling stands on several building blocks, with participatory modelling [19,20] and the development of modelling tools and software/hardware technologies considered as key pillars. In the context of policy analysis, more specifically, participatory processes are essential for linking science and policy [21] and to achieve the legitimacy of processes [22], with stakeholder participation and computer-based models regarded as key components of the participatory and collaborative modelling [17]. This society–science–policy interface [23] is usually moulded by different contextual pressures and communication protocols, thus rendering early stakeholder participation critical for successful outcomes in policy making [24].

There is no doubt that stakeholder participation in water resource management has received substantial attention in the last years [25–27]. The popularity of participatory modelling in particular has seen a substantial growth due to its compatibility with environmental paradigms such as Integrated Water Resources Management (IWRM) and Adaptive Management (AM) [28]. An advantage of participatory modelling resides in the potential to integrate meaningful input from decision makers and stakeholders into the modelling process [29]. Based on case studies from Africa, Asia, Europe and Oceania, Penny and Goddard [30] noted however that experimentation and learning beyond the "expert" group (to include non-expert participation) was mostly absent from discussions around model development.

To enable a participatory involvement, Basco-Carrera et al. [17] suggest that developed tools and models in the context of participatory modelling should be built using open source or freeware software where possible to facilitate distribution and use by stakeholders. Similarly, Carmona et al. [31] suggest that decision-making tools for successful stakeholder participation in natural resources management should be transparent, flexible and designed to elicit knowledge from different groups. Transparency and flexibility in the process of model development are also advocated by Bots et al. [21], with the aim of increasing stakeholders' trust by making the usually perceived "black box" model transparent.

Despite the clear need for stakeholder participation in the modelling development process, van Bruggen et al. [32] suggest that limited attention has been given to the modelbased exploration and design of policy pathways with stakeholders. They argue that disciplinary fragmentation and the "not-invented-here" academic syndrome ("a negative

attitude to knowledge that originates from a source outside the own domain" [33]) are factors hindering the development of modelling with stakeholders.

In arid and semiarid regions, collaborative processes and water governance are usually a major challenge [31,34–36] driven by contentious water use and competing stakeholder interests and values. This situation impacts the successful materialisation of the integration levels described by Kelly et al. [13] and poses challenges to the participatory modelling process as highlighted by Carmona et al. [31]. In particular, in these regions, the interaction between surface water and groundwater plays an important role [16,36,37]. This interaction is often complicated by agricultural and/or mining activities, as they will potentially alter the fragile flow regimes of the coupled water system [37]. As highlighted by Gorelick and Zheng [38], groundwater plays an important role, and its relevance will continue in the coming years, more importantly in arid and semi-arid regions.

This article describes the implementation of a participatory modelling process to develop and deploy an integrated modelling tool and digital platform (SimCopiapo) to support decision making in water management in a semiarid basin under contentious water use. The purpose of this digital platform is to explore alternative water management strategies to support scenario analysis and potential policy pathways with stakeholders, thus contributing to addressing the research need identified by van Bruggen et al. [32].

We build upon the work of Galvez et al. [34] to set up a participatory process in the Copiapó River Basin (CRB), northern Chile. We follow the integrated modelling levels suggested by Kelly et al. [13], and as such we include in the proposed integrated modelling tool surface water–groundwater interactions (*integration of processes and models*); local knowledge and expertise in water operational rules (*integration with stakeholders*); short-, mid- and long-term outputs, as well as sub-daily reservoir operations, daily water balance in irrigation districts and monthly time steps in groundwater assessment and different spatial scales for aquifer sectors and irrigation districts (*integration of scales of consideration*); and an agent-based model (ABM) to account for farmers' compliance against imposed caps on groundwater allocations (*integrated treatment of issues*). As suggested by the literature, we develop the integrated modelling tool and software platform in open source code with constant input from different stakeholder groups (water users, regulators, civil society, academy) [34] for transparency and flexibility [17,31] and to promote the legitimacy of the process [22] and ownership of results. The novelty of this work lies in advancing previous modelling efforts in the CRB [39–41] by improving on the operational rules of critical infrastructure in the CRB and co-designing water management strategies and impact indicators, all of which are designed with continuous input from key stakeholders by employing formal participatory and stakeholder engagement processes. A major feature of the proposed digital platform (SimCopiapo), compared to previous modelling efforts in the CRB, is the ability for users to run "on-the-fly" a loosely coupled [16] node–link water balance model and fully distributed groundwater model during interactive participatory sessions, thus facilitating social learning and knowledge co-creation. This was done in order to address research needs identified in the specialised literature [17,31,42]. Finally, the proposed digital platform (and integrated model) can be seen as a boundary object [43] bridging stakeholders and facilitating mutual understanding and cooperation—a practical exercise that has not been implemented before in the CRB [34].

The remainder of this article is arranged as follows. Section 2 describes the case study and the two-step methodological framework implemented in this work. Results of the integrated modelling process are analysed in Section 3 for a series of water management strategies and a base scenario. Section 4 presents a discussion of these results, and concluding remarks are offered in Section 5.

#### **2. Materials and Methods**

#### *2.1. Case Study: Copiapó River Basin*

The Copiapó River Basin (CRB) covers an area of 18,700 km<sup>2</sup> and is located in Northern Chile at the southern boundary of the Atacama Desert (Figure 1). The discharge contributions of the main tributaries Pulido, Jorquera and Manflas rivers in the headwater basins are regulated by the Lautaro Reservoir (26 Mm<sup>3</sup> ). The gauging station "Río Copiapó en La Puerta" shows an average discharge of 2.6 m<sup>3</sup> s −1 , whereas average annual precipitation in Copiapó city is 19 mm, reaching up to 500 mm y−<sup>1</sup> for altitudes over 5000 m above sea level (asl) [44]. The CRB is a clear example of a semiarid basin under sustained water stress originating from both natural and anthropogenic causes, where water management can be regarded as inadequate [40,41,45,46].

**Figure 1.** Location of the Copiapó River basin, main groundwater (aquifer) sectors, surface water irrigation districts (I: irrigation district D1, II: irrigation district D2, III: irrigation district D3, IV: irrigation district D4, V: irrigation district D5, VI: irrigation district D6, VII: irrigation district D7, VIII: irrigation district D8, IX: irrigation district D9. Most downstream districts (VIII and IX) are combined into irrigation district D89), and Lautaro Reservoir at the confluence of main tributaries (Jorquera, Manflas and Pulido rivers) to the Copiapó River. Urban areas in black color, coloured circles represent gauging stations (after [34]).

Currently, available surface water in nine irrigation districts (see Figure 1) is fully allocated for consumption, whereas the overexploitation of groundwater has been previously well documented in the literature [34,40,41] and is manifested by deteriorating groundwater quality and persistent deepening of groundwater levels. Administration of groundwater rights/licenses takes place in six groundwater/aquifer sectors (see Figure 1). Around 60% of groundwater demand is used for highly technified irrigation, whereas mining activities account for 30% and drinking water for 10% of the demand [47]. On an average water year, the Copiapó River dries up halfway to the outlet at the Pacific Ocean due to upstream water consumption and zero contributions from lateral intermediate sub-basins [48].

Groundwater rights for consumption total ca. 19 m<sup>3</sup> s −1 in the CRB, which contrasts with the estimated average recharge rates of 3 to 4.8 m<sup>3</sup> s <sup>−</sup><sup>1</sup> and effective groundwater demands of 6 to 14 m<sup>3</sup> s −1 [47,49]. Therefore, permanent conflicts between water users at different levels (upstream vs downstream users, surface water vs groundwater users) are detrimental factors for effective water resources management in the basin. These conflicts can be typified as Type 2, river basin conflicts, and Type 3, overexploited groundwater systems, by Bauer [50]. Rinaudo and Donoso [45] identified five factors leading to the current over-exploitation of Copiapo's groundwater resources: (i) limited knowledge of groundwater, (ii) legal complexity and political pressure, (iii) poorly-defined water permits, (iv) compliance and enforcement problems and (v) inconsistencies between management of surface water and groundwater.

Despite the peculiarities of the water resource management model in the Copiapó Basin [45,50], this management landscape is likely to reflect similar operational conditions as in other semiarid water-stressed basins around the world, thus providing generality to our findings.

#### *2.2. Methodological Framework*

We applied an intertwined two-step methodological framework in this work (Figure 2). First, we implemented a participatory process with existing key stakeholders to define and explore potential water management strategies and impact indicators of interest to stakeholders and to collect data on farmers' behaviours regarding groundwater regulation and operational rules for critical hydraulic infrastructure (c.f. Lautaro Reservoir). This step builds upon the work by Galvez et al. [34], who identified key stakeholders and barriers to collaborative water governance in the CRB and feeds into step 2. The second step consisted of designing and implementing a digital platform (termed SimCopiapo), which hosts a Graphical User Interface (GUI) (see Figure A1 in the Appendix A) with capabilities for stakeholders/modelers to run "on-the-fly" a node–link water balance model, a groundwater model and an agent-based model (ABM) [51,52]. This second step implemented the aspects identified through the participatory process in step 1. During participatory workshops, SimCopiapo was mainly run by stakeholders organised into groups with guidance provided by the research team. The purpose of this digital platform was to collaboratively explore different water management strategies as support for exploratory scenario analysis during participatory decision-making sessions with stakeholders. In the following sections, details for both steps are described.

**Figure 2.** General methodological framework highlighting development of the digital platform SimCopiapo, participatory processes and the integration of the surface water and groundwater models (red dashed line).

#### 2.2.1. Participatory Process

Galvez et al. [34] provide an overview of the stakeholders involved in the participatory process. We implemented 4 plenary workshops with 31 institutions and 5 working sessions with specific groups of stakeholders for data/local knowledge collection. In addition, we implemented two surveys with regional organisations in the study area: Water Resources Regional Advisory Committee (CARRH) (on-line) and Copiapó Exporters and Producers Association (APECO) (on-line). Surveys were used to collect information about farmers' behaviours regarding tolerance towards groundwater regulation (e.g., follow groundwater allocation rules) and the propensity to breach these rules following the social sub-model proposed by Castilla-Rho et al. [51,52]. This information was used to parameterise an agentbased model (ABM) to assess compliance against caps on groundwater use as explained in Section 2.2.3.

#### 2.2.2. Water Management Strategies and Impact Indicators

From the participatory process, we co-designed with stakeholders 13 water management strategies grouped in 4 domains: (a) exchanges of water uses/rights among users, (b) improvement to current hydraulic infrastructure, (c) management of groundwater recharge and (d) management of water demand. These water management strategies are shown and described in detail in Table 1 and were implemented in the SimCopiapo platform.

The participatory process also allowed the definition of a series of key impact indicators of interest to stakeholders. To this end, we followed a similar approach as that proposed by Santos Coelho [53]. The main impact indicators identified were (a) river flows through the Copiapó city (termed as urban flows), (b) environmental flows at the Copiapó basin outlet, (c) storage at Lautaro Reservoir (headwater basin), (d) percentage change in aquifer storage in groundwater sectors 2 to 6 after implementing water management strategies and (e) compliance with the cap on groundwater use. A description of the main impact indicators is presented in Table 2. These and other impact indicators (e.g., water security for individual irrigation districts) are automatically generated and exported by the SimCopiapo platform.






46



#### 2.2.3. SimCopiapo Platform

SimCopiapo is a digital platform built mostly in Python, including an API web and front-end HTML/JavaScript (Figure A1 in Appendix A). It includes capabilities to select pre-defined water management strategies devised by stakeholders, run "on-thefly" a loosely coupled node-link water balance and groundwater models [16,55], display graphs/plots/maps for rapid assessment and export reports summarising the results of management strategies selected by the user. It also contains an ABM associated with water management strategy 4.1 (groundwater demand management) to assess farmers' compliance against caps in groundwater use imposed by the regulator.

Figure 3 shows the interaction between the components of the SimCopiapo platform, including node–link water balance, groundwater and agent-based models. SimCopiapo uses geospatial information (irrigation areas, aquifer sectors, channel network, production wells, etc.), historical hydrological timeseries at the headwater basins for the period 1991–2016 and a series of alternative water management strategies (see Table 1) selected by the user to set up a specific scenario run. SimCopiapo users also have the opportunity to select/input pre-defined alternative hydrological time series driving the simulation (historical 1991–2016, 50% historical, etc.) or to include new time series if required. SimCopiapo is built as an open-source tool to allow continuous, replicable, reproducible and transparent research and improvements by other users [56–58], thus improving on previous efforts developed under proprietary hydrological software [39–41,49].

**Figure 3.** SimCopiapo digital platform and interactions of the integrated surface water and groundwater models and different modules and platform components.

It is worth noting that the objective of the digital platform is to facilitate the interaction among stakeholders in the contentious Copiapó river basin, where competing and conflicting water uses exist and the level of collaboration for water management is limited and generally perceived as inadequate [34]. Therefore, the focus is on providing a research tool able to run basin-scale assessments in order to support rapid appraisal of water management strategies enabling stakeholder discussion, collaboration and decisionmaking. Under these premises, SimCopiapo is aligned more closely to what Oxley et al. [24] define as a policy-oriented model, where *accurate* process representation is traded for *adequate* process representation and emphasis is focused on addressing practical policy issues. For this purpose, we illustrate how the digital platform can be populated with surface water and groundwater models previously documented and validated by stakeholders (e.g., [39–41,49]).

DGA-HIDROMAS [39] provides the latest surface water model available for the Copiapó basin implemented in the AQUATOOL software—a proprietary generic decision support system for water resources planning [59]. We translated the topology of the AQUA-TOOL model for the Copiapó basin into a node–link model coded in Python, accounting for the daily mass balance of the surface water system. The conceptual representation of this topology and the operation of the irrigation districts in the Copiapó River Basin is presented in Figure 4. This figure shows the upstream Lautaro reservoir as the main hydraulic infrastructure regulating surface flows to supply irrigation water to nine districts (D1–D9, irrigation districts D8 and D9 are combined into a single district, D89. See Figure 1 for locations of irrigation districts). Available surface water regulated from Lautaro reservoir is equally allocated between districts D1 to D7 (12% each), whereas districts D8 and D9 are allocated 8% each (i.e., 16% combined for D89). Figure 4 also shows the crop sectors (e.g., R2a-XX) belonging to each irrigation district, with some of the irrigation districts including more than one crop sector (e.g., R2a-13 and R2a-14 belong to irrigation district D6). It is worth emphasising that D8 and D9 are the most downstream irrigation districts located at the outskirts of Copiapó city, thus experiencing changes in land use patterns from rural to urban. Upstream of the gauging station "Copiapo @ Mal Paso", most of the available surface water is conveyed through a channel (1000 L/s maximum capacity), leaving just excess water flowing through the natural river course.

**Figure 4.** Conceptual model for the (**a**) operation of the Lautaro Reservoir and the (**b**) irrigation districts in the Copiapó River Basin. Values for parameters of the conceptual model are obtained from [39,40] and included in the node–link water balance model. Water sources for each irrigation district are distinguished between C: channel, W: wells and M: mixed source and are based on crop area supplied by that source.

For each irrigation district, information on crop types/surfaces, irrigation and water conveyance efficiencies, water sources, etc. is used to perform an internal supply–demand balance per irrigation area, considering alternative water sources (surface water, groundwater or combination). The water source indicates the water volume supplying crop areas for each irrigation district. The volumes demanded for each irrigation district are obtained from the crop surveys by DGA-DICTUC [49] and the water licenses. In this way, irrigation demands are supplied first by channel source, then mixed source and finally groundwater (wells). Results of the water balance for each irrigation district are spatially coupled to the aquifer sectors implemented in the groundwater model, which together with the infiltration through the riverbed define the main recharge rates for the aquifer sectors (see Figure 1).

This node–link mass balance model was fully coupled with the latest available groundwater model for the Copiapó aquifer developed in MODFLOW-2005 [60] by DGA-HIDROMAS [39]. Using the FloPy Python package [61], we translated this MODFLOW-2005 model for operation in Python and coupled it with the node–link model of the surface water system in the SimCopiapo platform. The daily node–link water balance model was aggregated to a monthly time-step for consistency with the MODFLOW model for the Copiapó aquifer. For full details on the coupling process we refer the reader to [62].

As shown in Figure 3, SimCopiapo also includes an agent-based model (ABM) to assess farmers' compliance against caps imposed on groundwater use in the Copiapó aquifer. This ABM is based on the social sub-model developed by Castilla-Rho et al. [51,52], which represents a social utility function, S, that follows a Cobb–Douglas functional form:

$$\mathbf{S} = \text{grid}^{\text{m}} \text{ (1 -- group)}^{\text{n}} \tag{1}$$

where m = number of times a farmer reports a neighbour taking groundwater illegally, n = number of times a farmer is seen taken groundwater illegally, and grid-group are categories of the Cultural Theory proposed by Douglas [63]. S (social utility function) represents the loss of social reputation and the social costs to groundwater users when reporting noncompliant behaviour. Using survey data collected from farmers in the Copiapó basin [62] and the four grid-group categories (Egalitarian–Hierarchist–Individualist–Fatalist) proposed by Douglas [63] in Cultural Theory, we were able to parametrise equation 1 and thus farmers' decision-making processes. The user of SimCopiapo can adjust two parameters associated with the ABM model: (a) the percentage of groundwater users monitored by the regulator to check compliance and (b) the severity of the fines (as a percentage of the total farm revenue) if the farmer is caught taking groundwater illegally.

Equation (1) quantifies the loss of social reputation and the social costs to farmers when reporting non-compliant neighbours engaged in illegal extraction of groundwater in the Copiapó basin, and thus impact farmers' future decisions of engaging in non-complaint behaviour (i.e., taking groundwater illegally). Other factors impacting this decision relate to farmers' probability of being monitored by the regulator and the severity of fines if farmers are caught in non-compliant behaviour [48,49]. This social metric is combined with an economic (gross margins from crop enterprise) and institutional (monitoring/monetary fines) score into each farmers' objective function for decision making; i.e., whether to take groundwater illegally or not. For details on this implementation, the reader is referred to [51,52].

#### 2.2.4. Improvements on Previous Integrated Modelling Tools

Suarez et al. [40] presented an integrated model for the CRB using the SIMGEN module of AQUATOOL software [59] (surface water only). More recently, Hunter et al. [41] presented an integrated model for the CRB based on [40] coupling the WEAP (Water Evaluation and Planning system) model [64] and the MODFLOW [60] groundwater model described by DGA-HIDROMAS [39]. Although these works claim the advantages of their corresponding integrated modelling frameworks, both tools rely on proprietary software and are therefore not amenable for rapid modifications by interested stakeholders, have been developed with limited input from key stakeholders in terms of potential

water management strategies as well as impact indicators of interest to stakeholders, and concentrate only on demand management strategies. In this work, we improved on several aspects on the integrated modelling framework for the CRB: (1) the groundwater model developed by DGA-HIDROMAS [39] has been checked for spatial and temporal consistency of aquifer contributions to surface water at La Puerta and Angostura gauging stations; (2) evapotranspiration and groundwater demands were revised and activated; for the surface water model, (3) Lautaro Reservoir operational rule for water allocation was completely re-designed and implemented thanks to the advice of the Vigilance Board of the Copiapó River and its Tributaries; and (4) the operation of the irrigation districts was revised and deployed in Python considering (a) controlled discharge from the Lautaro Reservoir and (b) a supply–demand model for the irrigation districts considering allocation volumes, irrigation demands, irrigation losses, conveyance losses and gross water demand. For details on other improvements regarding updates on mining groundwater demands, downstream irrigation districts, drinking water demands and losses in the potable water network, the reader is referred to [62].

#### **3. Results**

#### *3.1. Results Participatory Process*

A total of 31 organisations representing the civil society (6), regional state agencies (16) and private/productive (9) sectors were engaged in the participatory process. On average and across all participatory workshops and working sessions, stakeholders from the civil society were the least involved (35% of organisations engaged), whereas stakeholders of the private sector were the most engaged, with 52% of the institutions of this sector taking part in the participatory sessions.

Regarding the on-line surveys for parameterising the ABM for groundwater regulation, the CARRH survey showed that civil society stakeholders were the most engaged, with 83% of the institutions of this sector providing responses, whereas only 31% of the institutions of the public sector were engaged in this process. Forty-four percent of private sector stakeholders participated in this survey. The APECO survey targeted 25 farmers of the CRB, with more than 83% concentrated in the upstream aquifer sectors 1 to 3 and the remaining 17% concentrated in downstream aquifer sectors 4 and 5. No responses were obtained from farmers in aquifer sector 6. Although not shown here, results of both surveys indicate a clear trend towards validating the importance of regulating groundwater resources for sustainable use and minimizing impacts to ecosystems and third parties, and enforcing this regulation in practice. Discrepancies among the CARRH stakeholders were observed on justifying the illegal extraction of groundwater on economic (profits) grounds and allocating importance to social costs (loss of reputation) if caught breaching caps on groundwater use imposed by the regulator. This discrepancies can be attributed to the heterogeneity of the stakeholders composing the CARRH [34]. On the contrary, the APECO survey indicated that farmers attributed a higher importance to the loss of social reputation (individual) if caught breaching the imposed cap on groundwater use and allocated a higher importance to collective enforcement policy such as effective monitoring of groundwater use. For details about the outcomes of the surveys and the implementation in the ABM model, the reader is referred to [50,51,60].

#### *3.2. Validation of Node–Link Water Balance Approach to Surface Water Modelling in SimCopiapo*

As the main driver controlling surface water flows and the recharge to aquifers is the operation of the Lautaro Reservoir, we focused on reproducing the observed discharges measured at the gauging station immediately downstream of the Lautaro Reservoir (Figure 5). For the simulated period implemented in SimCopiapo (1991–2016), we observe a much better correspondence between observed and simulated discharges compared to the original surface water models implemented by [40,41]. In general, peak releases are properly simulated with the exceptions of hydrological years 1997/1998 and 2001/2002,

where both the original surface water model and the node–link model implemented in SimCopiapo show difficulties in capturing the peak behaviour.

**Figure 5.** Release discharge from Lautaro Reservoir for the (**A**) original model developed by Hunter et al. [41] and Suarez et al. [40] and (**B**) adapted operational rule implemented in SimCopiapo.

Similarly, Figure 6 shows the simulated and observed discharges for the period 1991– 2016 in La Puerta and Copiapó City gauging stations (see Figure 1 for locations). This figure shows that the operation of the Lautaro Reservoir and the improved supply–demand model for the irrigation districts implemented in the node–link water balance model in the SimCopiapo platform can reproduce the observed discharges in a reasonable manner, preserving the long-term trend and capturing relevant peaks in 1998/1999 and 2003/2004. It is worth noting that "Rio Copiapó en Ciudad" (Figure 6b) is a gauging station located downstream of all irrigation districts, and as such reflects excess volumes after irrigation use located upstream of Copiapó city.

**Figure 6.** Simulated versus observed discharges in (**a**) "Rio Copiapó en La Puerta" and (**b**) "Rio Copiapó en Ciudad" gauging stations. Horizontal red line represents the average controlled discharge from Lautaro Reservoir (3020 L/s).

Closing the water balance obtained from loosely-coupled surface water and groundwater models has been identified as a drawback of the integration process [16,55]. La Puerta gauging station in the CRB is the main control point to verify that the coupling

process of both the daily node–link mass balance model (implemented in Python) and the MODFLOW model (implemented in FloPy) is correct. In this sector, the Copiapó valley shows an important reduction in its cross section, and basement rocks are uplifted, thus substantially constraining the aquifer cross section, resulting in substantial contributions from the aquifer to the Copiapó River [39,49]. In general, discharges at La Puerta gauging station are influenced by the groundwater throughflows from groundwater sector 2 and the infiltration rates from the Lautaro Reservoir (see Figure 1). The infiltration losses in the Lautaro Reservoir were therefore adjusted between 500 L/s and 3500 L/s (as a function of the stored volume) until a reasonable match between the daily node–link mass balance model and groundwater model outputs at La Puerta was obtained. Figure 7 shows the match between both the node–link water balance and MODFLOW models at La Puerta gauging station. In general, we observe a good match between river gains from groundwater simulated through the drain package of MODFLOW and the node balance at La Puerta gauging station. We observe a good fit when simulating the temporality and magnitude of the time series, with a range between 1000 and 2500 L/s at "Rio Copiapó en La Puerta". Few discrepancies are observed at the start of the simulation period, most likely attributed to the stabilisation of parameters in the node–link model (warm-up period) and due to the aggregation of the monthly time steps into 6 month stress periods in MODFLOW.

**Figure 7.** Validation of coupling process for daily surface water mass balance models and groundwater model at "Rio Copiapó en La Puerta" gauging station.

#### *3.3. Results for Individual Water Management Strategies*

To assess the results of individual (and combined) water management strategies we defined a base scenario reflecting hydraulic infrastructure, water demands, crop types and land uses corresponding to year 2018. This base scenario reflects a business-as-usual (BAU) approach. Both the base scenario and water management strategies are assessed for a 25-year period (2018–2042) in order to isolate the marginal impacts of implementing such strategies, using the observed hydrology for the period 1991–2016 as forcing data.

Table A1 in the Appendix A shows the individual results for each water management strategy described in Section 2.2.2. Individual water management strategies show spatially bounded impacts and marginal cumulative impacts at the end of the simulation period. Strategy 3.1, for example, shows an increase in infiltration flows in the Copiapó river of less than 3% the potential recharge volume due to constraints in the size of the recharge ponds and the available surface flows for infiltration. In the next sections, we analyse a selected group of results for individual strategies.

53

#### 3.3.1. Water Uses/Rights Exchanges

For strategies promoting water uses/rights exchanges among users, the most attractive corresponds to strategy 1.3-c in Table 1, which promotes urban flows increases up to 263 L/s on average for the 25-year simulated period. Figure 8 shows the monthly frequency of the occurrence of urban flows for the base scenario and strategy 1.3-c. Stakeholders have defined the occurrence of urban flows through Copiapó city as an important indicator of quality of life, and results show a substantial increase in the occurrence of urban flows from 19% to 96% by implementing strategy 1.3-c.

**Figure 8.** Monthly frequency of occurrence of urban flows at the Copiapó City gauging station for (**a**) base scenario and (**b**) strategy 1.3-c, for dry, acceptable, and flooded thresholds. Blue cells record the occurrence of average monthly discharges greater than 0 L/s and less than 5000 L/s at Copiapó city.

#### 3.3.2. Improvement to Current Hydraulic Infrastructure

In terms of water management strategies promoting improvements to current hydraulic infrastructure, strategy 2.1 (impermeabilisation of Lautaro Reservoir) shows substantial impacts (positive and negative). Figure 9 shows the water balance for the Lautaro Reservoir for both the base scenario and strategy 2.1. For this figure, we observe that releases from the Lautaro Reservoir become regular and over 2000 L/s, with 18 out of 20 water years using the spillway to regulate the reservoir's capacity. After fully implementing the impermeabilisation of the inundated surface by year 5, infiltration losses become 0, and the reservoir volume is above 50% its capacity for 14 out of the 20 remaining years. Although not shown here, this results in increases in water security for irrigation districts no. 6, 7, 8 and 9, which are most closely located downstream of the reservoir. In addition, by implementing strategy 2.1, an increase in the frequency of occurrence of urban flows through Copiapó city from 18% to 30% of the months in the simulation period 2018–2042 is observed.

Despite the positive impacts from the surface water perspective, a negative impact is observed for the groundwater sector 2, located immediately downstream of the Lautaro Reservoir. Groundwater sector 2 is a narrow tube-like aquifer recharged mainly through upstream groundwater throughflows originating from upstream aquifers and, most importantly, infiltration losses from the Lautaro Reservoir. Figure 10 shows the groundwater levels of representative observation wells located in groundwater sector 2. After implementing the impermeabilisation of the Lautaro Reservoir, a sustained decreasing trend is observed in the mid- (MT) and long-term (LT), reaching average values of −0.8 m/y. It is worth noting that these decreasing trends concentrate in the upstream half of groundwater

sector 2, whereas groundwater levels around La Puerta remain stable or increase given the constriction of the aquifer section explained in sections above.

**Figure 9.** Water balance for the Lautaro Reservoir for (**a**) base scenario and (**b**) strategy 2.1, reflecting the impermeabilisation of 100% inundated surface after a 5 year period (vertical red dotted line).

**Figure 10.** Time series of groundwater levels in representative observation wells located in the groundwater sector 2 for (**a**) base scenario, and (**b**) Strategy 2.1. ST: short-term, MT: mid-term, LT: long-term. Increasing average trend for each period for all observation wells in blue. Decreasing average trend for each period for all observation wells in red.

#### 3.3.3. Demand Management

For the strategies promoting management of the groundwater demand, strategy 4.1 (a proportional reduction of 40% of groundwater use across groundwater sectors 3, 4 and 5) shows substantial impacts. This strategy was implemented through the ABM and assessed the level of compliance of the imposed cap on groundwater use and the impact on groundwater level/balance. Based on our previous experience, both monitoring of groundwater users and fine levels are strong deterrents when dealing with non-compliant behaviour in groundwater management [48,49]. Figure 11 shows two levels of enforcement tested in this strategy: (a) monitoring of 90% of users in groundwater sectors 3, 4 and 5 and substantial fine levels (90% gross profit from farm enterprise) if farmers are caught breaching the cap, thus defining a strong enforcement policy; and (b) lax monitoring (20% of users in groundwater sectors 3, 4 and 5) and fine levels (20% gross profit from farm enterprise), thus defining a weak enforcement policy. Figure 11a shows that for groundwater sectors 3, 4, 5 and 6, implementing the cap on groundwater use together with a strong enforcement policy brings storage volumes in these aquifer sectors back to values that are better than the initial state of the base scenario. On the contrary, when the enforcement policy is weak (Figure 11b), the impacts are limited to aquifer sector 5 and to a lesser extent in aquifer sector 6 given the proportional volume of groundwater for irrigation in these sectors. This indicates that at least 20% of the groundwater users of aquifers sectors 3, 4 and 5 need to be monitored if a cap on groundwater extraction is imposed by the regulator. For other impact indicators, implementing strategy 4.1 has a limited impact.

**Figure 11.** Impact indicators for strategy 4.1 (cap on groundwater use implemented through ABM) for two enforcement strategies: (**a**) high level of monitoring (90% of groundwater users) and fines (90% of revenue if farmer caught breaching the cap) and (**b**) low level of monitoring (20% of groundwater users) and fines (20% revenue if farmer caught breaching the cap).

#### *3.4. Results for Combined Water Management Strategies*

Any single water management strategy cannot address the basin-scale water management challenges identified in the CRB by different authors (see e.g., [34,40,41,45]). Impacts of individual strategies are sectoral and, in some cases, spatially and temporally constrained. It then seems appropriate to combine alternative water management strategies to assess multiple stakeholders' interests/perspectives. Based on the participatory process, a prioritised combination of water management strategies attractive to stakeholders of the CRB was defined by simultaneously implementing strategies 1.3-c, 2.1, 2.3, and 4.1a (see Table 1). Strategy 2.3, which corresponds to operating a desalination plant to supply 90 L/s (first

5 years), 450 L/s (next 5 years) and 930 L/s (last 15 years), has been included as this strategy is already in early operation in the CRB.

Figure 12 shows the results for impact indicators when combining these strategies. We found a substantial increase in the occurrence of urban flows from 20% to ca. 100% of the simulated period (25 years), thus enhancing the quality of life perceived by stakeholders of the CRB; marginal increases in the environmental flows at the basin outlet, thus promoting a healthy habitat for the wetland at the Copiapó River mouth; substantial increases in the storage volume of the Lautaro Reservoir from 20% to 72% of the simulated time with volumes greater than 50% its maximum storage capacity, thus impacting water security for irrigation districts; recoveries in heads and stored volumes in groundwater sectors 3, 4, 5 and 6, thus decreasing pumping costs to users and contributing to groundwater sustainability in the long-term; high levels of compliance (>80% groundwater users) to caps on groundwater use supported by a robust enforcement policy formulated around high monitoring rates and substantial fines.

**Figure 12.** Impact indicators for the combinatory of water management strategies prioritised through the participatory process.

All these positive impacts also carry a negative impact, which is the detrimental impact on groundwater heads and stored volumes in the upper section of the aquifer sector 2, immediately downstream of the Lautaro Reservoir. Figure 13 shows that after implementing the combination of water management strategies, decreases in groundwater heads in the upper section of aquifer sector 2 can reach up to 40 m compared to the base scenario. Long-term increases in groundwater heads in sectors 3 and 4 can reach between 30 m and 40 m in aquifers immediately downstream of La Puerta gauging station and around 10 m in aquifer sector 5 downstream of Copiapó City gauging station.

**Figure 13.** Groundwater head difference (m) between base scenario and combined water strategies. Inside panels show results of Lautaro Reservoir water balance (storage (% maximum volume) and infiltration losses (L/s)) and river flows (L/s) at the outlet of the CRB.

#### **4. Discussion**

The situation in the CRB indicates human and environmental vulnerabilities [38] stemming from the sustained exploitation of groundwater resources. Several authors have tried to explain the factors driving the water crisis in the CRB from economic, regulatory and management perspectives [34,45,46,48], whereas others have suggested technical solutions such as basin-scale or sectorial groundwater use restrictions [40,41,49]. Evidence by Wurl et al. [65] in a similar context (arid overexploited aquifer in Mexico) suggests that water management problems can no longer be addressed purely as technical problems and should consider a wide range of stakeholders' perspectives to strengthen the resilience of water resources. In this work, we have contributed towards this by devising potential water management strategies (i.e., technical solutions) and relevant impact indicators through a bottom-up participatory process driven by key stakeholders in the CRB.

Our approach directly addresses one of the tasks Rinaudo and Donoso [45] suggest the regulator should implement as part of a groundwater management model; i.e., implement an efficient enforcement strategy, e.g., by proposing minimum monitoring coverage and fine levels to achieve compliance on cap reductions based on ABM results. Although not explicitly addressed by our integrated modelling approach, the remaining tasks defined by Donoso and Rinuado (i.e., calculate sustainable groundwater abstraction limits, defining sharing rules, reallocation of water use rights and rules to adjust volume of water use rights) can be assessed implementing minor modifications to the SimCopiapo tool (e.g., test rules to adjust volumes of groundwater rights in different aquifer sectors).

SimCopiapo can be classified as a policy-oriented model [24] where adequate process representation, addressing practical policy issues and supporting decision-making with stakeholder participation are regarded as key features [42]. One of the purposes of the participatory modelling was social learning and acceptance of model improvements through direct participation in designing the conceptual model, water management strategies and impact indicators for discussion. Following the classification of Hare [42] the participatory modelling exercise implemented in SimCopiapo aligns with a Front and Back-End (FABE) category, where stakeholder involvement concentrates on early (conceptual model design, definition of operational rules, impact indicators and water management strategies) and later stages (assessment of water management strategies, discussion of potential policy pathways) of the modelling process. The effectiveness of the methods used in this participatory process is yet to be asserted through a follow-up process with decision and policy-makers of the CRB. A promising research direction is a post-hoc assessment using boundary objects attributes (credibility, salience, legitimacy) as suggested by Falconi and Palmer [43] to assess the success of SimCopiapo as a participatory model.

Our integrated approach and proposed digital platform (SimCopiapo) contribute to addressing the challenges identified in the use of models to operationalise IWRM by Badham et al. [66]. First, the bottom-up participatory process contributed to addressing a difficult problem in water policy by streamlining multiple pressures, conflicting values, competing goals and limited resources in a transparent way; and second, it helped in handling the human element in IWRM by reconciling potential conflictive agendas by stakeholders. The latter has been recognised as an important research avenue in water management [46,67,68].

Results show that no single strategy is able to provide definite long-term solutions to the water management challenges observed in the CRB nor to capture the multiplicity of stakeholders' interests expressed through the impact indicators identified. DGA-DICTUC [49], Suarez et al. [40] and Hunter et al. [41] proposed basin-wide or sectoral reductions in groundwater use (demand management) by values between 20% and 50% on the basis of cost analysis or a multi-dimensional measure of sustainability. While useful, assuming monetary motivations are central to water management and a key driver of behavioural change in groundwater users ignores the role that social, ecological and cultural values might have in this regard, thus constraining the assessment of water management strategies [69,70]. SimCopiapo contributes to equilibrating the assessment by transparently assessing physical, ecological and social aspects of the water management strategies devised, thus counter-balancing the bias towards exclusively cost-based assessments observed in the literature (e.g., [71–73]).

Our results indicate that management of groundwater use is one of the most critical water strategies to recover the aquifer sectors in the CRB. However, there needs to be a clear enforcement policy to trigger the minimum momentum required to achieve social acceptance of this policy. This is fully aligned to one of the drivers suggested by Rinaudo and Donoso [45] triggering the water crisis in the CRB. Results indicate that there seems to be a middle point between lax and strong enforcement policies to achieve this reduction in demand in a sustainable way. The definition of where this middle point lies is beyond the scope of this article, but our results bring a first approximation to this; i.e., between 20% and 90% monitoring coverage and between lax (fines accounting for 20% revenue) and strong (fines accounting for up to 90% revenue) fine levels for breaching the imposed cap on groundwater use. These results are fully aligned with findings by Castilla-Rho et al. [51,52] for other aquifers around the globe.

Future research avenues might consider including crop choices in the supply–demand and ABM models, implementation of other water management strategies and impact indicators, optimising the level of monitoring and fines to achieve a target compliance in different aquifer sectors or groundwater user groups (e.g., mining, industry) and testing multi-level ABM parametrisations for time-varying water management policies as in Du et al. [74].

#### **5. Conclusions**

In this work, we demonstrated the value of combining participatory modelling and integrated modelling to develop a digital platform tool (SimCopiapo) supporting social learning and knowledge co-construction for water management in a semiarid basin under contentious water use. We have contributed to transparently positioning a policy-based model as a boundary object to bridge a diverse group of stakeholders with individual and competing interests and perspectives.

Current water management in the Copiapó River Basin (CRB) is not providing the required solutions for resource sustainability, with previous research to date proposing (optimised) cost-based technical solutions focusing solely on groundwater demand reductions. This narrow perspective however does not seem to hold for complex water management problems with no single/simple solutions that act and depend on values and priorities by multiple stakeholders. Early stakeholder engagement and participation for social learning and knowledge co-construction are therefore essential steps in this process.

Our results suggest that management of groundwater demand in the CRB together with a portfolio of strategies including water rights exchanges, improvements to hydraulic infrastructure and robust enforcement policies are best suited to capture the diversity of stakeholders' interests and perspectives when addressing the water management challenges observed in the CRB. This diversity is expressed through a series of impact indicators, which, directly or indirectly, are a reflection of not only available groundwater resources for use but also ecological (e.g., basin outlets) and social (e.g., urban flows) aspects of relevance to stakeholders. An important aspect to manage groundwater demand is the establishment of an efficient enforcement policy to monitor caps imposed on groundwater use. In the absence of a clear policy and an institutional/legal framework to achieve this, water users' behaviours will continue to be non-cooperative, therefore leading to unsustainable groundwater use in the long-term.

Finally, we can conclude that including stakeholders in the participatory modelling process has led to the definition of a rich and diverse collection of water management strategies and ways to assess these strategies, which are a good reflection of stakeholders' interests and visions. This indicates that not only the final product—i.e., the SimCopiapo digital platform—is of value but also the process leading to its creation.

**Author Contributions:** Conceptualisation, R.R., J.C.-R. and G.B.; methodology, R.R. and J.C.-R.; software, R.B. and J.C.-R.; validation, R.R., J.C.-R., G.B. and C.P.; formal analysis, R.R. and J.C.-R.; investigation, R.R., J.C.-R., G.B. and C.P.; resources, E.C.; data curation, R.R., J.C.-R., G.B., R.B. and C.P.; writing—original draft preparation, R.R.; writing—review and editing, R.R., J.C.-R., G.B., R.B., E.C. and C.P.; visualisation, R.R. and J.C.-R.; supervision, R.R. and E.C.; project administration, E.C. and G.B.; funding acquisition, R.R., G.B., J.C.-R. and E.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Atacama Regional Government (Chile) through the "Fondo de Innovación para la Competitividad", FIC-R 2016, código BIP N◦ 30486475.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study, and clearance from the CSIRO Ethics Committee was obtained for this research (097/18).

**Data Availability Statement:** All data collected and necessary to reproduce this research are available from CSIRO Research Data Planner (rdp.csiro.au) and the corresponding project website https://research. csiro.au/gestion-copiapo/en/projects-simcopiapo-participative-modelling-for-water-management/, accessed on 14 March 2022.

**Acknowledgments:** We are grateful to the many regional/local stakeholders who took part in the participatory sessions throughout this research. Authors also thank Alexis Medina and Victor Galvez for comments made to improve earlier versions of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Graphical user interface for SimCopiapo platform v1.0 (in Spanish).

**Table A1.** Summary of impacts for water management strategies implemented in SimCopiapo.



**Table A1.** *Cont*.

(oo): no decreases/increases over the simulation period with respect to base scenario; (– –): decreases over the simulation period more than 20% with respect to base scenario; (o−): decreases over the simulation period less than 20% with respect to base scenario; (o+): increases over the simulation period less than 20% with respect to base scenario; (++): increases over the simulation period more than 20% with respect to base scenario.

#### **References**


### *Article* **Inferring Hydrological Information at the Regional Scale by Means of** δ **<sup>18</sup>O–**δ **<sup>2</sup>H Relationships: Insights from the Northern Italian Apennines**

**Federico Cervi 1,\* and Alberto Tazioli <sup>2</sup>**


**Abstract:** We compared five regression approaches, namely, ordinary least squares, major axis, reduced major axis, robust, and Prais–Winsten to estimate δ <sup>18</sup>O–δ <sup>2</sup>H relationships in four water types (precipitation, surface water, groundwater collected in wells from lowlands, and groundwater from low-yield springs) from the northern Italian Apennines. Differences in terms of slopes and intercepts of the different regressions were quantified and investigated by means of univariate, bivariate, and multivariate statistical analyses. We found that magnitudes of such differences were significant for water types surface water and groundwater (both in the case of wells and springs), and were related to robustness of regressions (i.e., standard deviations of the estimates and sensitiveness to outliers). With reference to surface water, we found the young water fraction was significant in inducing changes of slopes and intercepts, leading us to suppose a certain role of kinetic fractionation processes as well (i.e., modification of former water isotopes from both snow cover in the upper part of the catchments and precipitation linked to pre-infiltrative evaporation and evapotranspiration processes). As final remarks, due to the usefulness of δ <sup>18</sup>O–δ <sup>2</sup>H relationships in hydrological and hydrogeological studies, we provide some recommendations that should be followed when assessing the abovementioned water types from the northern Italian Apennines.

**Keywords:** stable water isotopes; young water fraction; global meteoric water line; northern Italian Apennines

#### **1. Introduction**

Oxygen (18O and <sup>16</sup>O) and hydrogen isotopes (2H and <sup>1</sup>H) of water are commonly used in surface and subsurface hydrology [1–3]. They are considered environmental tracers in the form of δ <sup>18</sup>O and δ <sup>2</sup>H, where δ(‰) = h *Rsample Rstandard* − 1 <sup>×</sup> <sup>1000</sup><sup>i</sup> and *R* is the corresponding isotopic ratio (18O/16O or <sup>2</sup>H/1H) in the water sample or in the standard (usually V-SMOW, i.e., the Vienna Standard Mean Oceanic Water). If a multitude of water samples is collected from the same source (a rain gauge, a river, a spring, or a well), the corresponding δ <sup>18</sup>O–δ <sup>2</sup>H pairs in a Cartesian graph will be aligned along a regression line in the form of y = mx + q, where y is δ <sup>2</sup>H, x is δ <sup>18</sup>O, m is the slope, and q the intercept. This fact was first noted by [4] when reporting several δ <sup>18</sup>O and δ <sup>2</sup>H values from precipitation waters worldwide, which allowed definition of the so-called "global meteorological water line" (GMWL) equal to δ <sup>2</sup>H(‰) = 8.0 · <sup>δ</sup> <sup>18</sup>O + 10.0. The authors of [5] have substantially confirmed this slope and intercept by processing more recent data. Although this relationship is valid everywhere, more accurate regression lines (with different slope and intercept) can be obtained by selecting δ <sup>18</sup>O and δ <sup>2</sup>H values from restricted areas. This is due to the specific isotopic fractionation processes (i.e., vapour pressure and temperature conditions) controlling

**Citation:** Cervi, F.; Tazioli, A. Inferring Hydrological Information at the Regional Scale by Means of δ <sup>18</sup>O–δ <sup>2</sup>H Relationships: Insights from the Northern Italian Apennines. *Hydrology* **2022**, *9*, 41. https:// doi.org/10.3390/hydrology9020041

Academic Editors: Il-Moon Chung, Sun Woo Chang, Yeonsang Hwang and Yeonjoo Kim

Received: 6 January 2022 Accepted: 16 February 2022 Published: 21 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

precipitation over each area. Moreover, since further post-condensation and temperaturedriven processes such as evaporation and evapotranspiration could act prior to infiltration and/or during runoff, δ <sup>18</sup>O–δ <sup>2</sup>H regressions from rivers (river water lines, RWLs) and groundwater (groundwater lines, GWLs) may also differ from that of the precipitation occurring in their recharge areas (meteoric water lines, MWLs). In fact, evaporation and evapotranspiration lead to a fractionation between the different isotopologues of water, with lighter water molecules (1H<sup>2</sup> <sup>16</sup>O) vaporising faster than heavier ones (2H<sup>2</sup> <sup>18</sup>O) and inducing an enrichment of the latter into the liquid residual. In this case, as a water parcel evaporates, its isotopic composition usually evolves with a δ <sup>18</sup>O–δ <sup>2</sup>H regression line whose slope is lower than those of MWLs.

It is evident that such changes in slopes from RWLs and GWLs concerning those of MWLs can be used to infer information on the hydrological processes occurring at the slope and the catchment scales. As an example, and without claiming to be exhaustive, the following studies highlighted that a change in slope from RWLs and GWLs concerning those of MWL can be used to:


Starting with the pioneering works by [4,12] regarding meteoric water and up to this day, δ <sup>18</sup>O–δ <sup>2</sup>H regression lines are usually carried out by means of the ordinary least squares (OLS) method, i.e., an approach that minimises the sum of the squared vertical distances between the y data values and the corresponding y values on the fitted line (the predictions). Thus, the OLS design assumes that there is no variation in the independent variable (x) and is considered as the simplest method among the several available linear regression models. By focusing on the aforementioned isotopes of water, we should also take care of the variations associated with variable (x) as the same or different isotopic fractionation processes, which may have developed even at different rates and may have affected both δ-values.

For this reason, [13] proposed a more complex linear regression approach for obtaining MWLs lines that consider associated errors with both dependent (y) and independent variables (x), such as the reduced major axis (RMA) and the major axis (MA). In the end, it is found that MA approaches usually led to the smallest discrepancies between the estimated and predicted values (a measure of goodness of fit usually described with the well-known coefficient of determination R<sup>2</sup> , i.e., smallest discrepancies are identified by higher values of R<sup>2</sup> ) and larger slopes in MWLs calculation than those obtained with RMA and OLS [14]. The recent attempt made by [15] on several water types (river water, groundwater, soil water) confirmed that RWLs and GWLs were also characterised by larger slopes in the case of MA regression while the lowest values were noted when using OLS. The same authors found that in all cases (MWLs, RWLs, GWLs), the higher the R<sup>2</sup> between δ <sup>18</sup>O and δ <sup>2</sup>H values, the smaller were the differences in the slopes obtained by MA, RMA, and OLS. This was due to the different sensitiveness of the linear regression approaches to outliers and large measurement errors rather than temperature-driven post-condensation processes.

This study aims to verify whether such discrepancies in slopes and intercepts from different regression methods are present (thus significant) or not in four water types (precipitation, surface water, groundwater collected in wells from the lowlands, groundwater from low-yield springs) from the northern Italian Apennines. For this reason, we exploited datasets already published in the literature, e.g., [7,11,16–18], and we carried out visual inspection (heat maps) and statistical comparison of results from the three aforementioned approaches (OLS, RMA, MA) already tested in [15]. With reference to OLS approach, we further verified whether preliminary weighting of the isotopic data to the corresponding values of discharge or precipitation may have induced changes to our results or not. In addition, we tested two methods, namely, Prais–Winsten (PW) and robust (R), to investigate possible influence on the final δ <sup>18</sup>O–δ <sup>2</sup>H alignments of nonstationary processes (here, we recall that the OLS, RMA, and MA approaches are based on the assumption that data are not serially correlated, thus a δ <sup>18</sup>O–δ <sup>2</sup>H pair from a determined time period is not correlated with the earlier one, while properties as mean, variance, and autocorrelation are constant over time, i.e., stationary, while recent papers in the literature highlighted the possible nonstationary behaviour of such series of isotopic data [19]) or even outliers (i.e., anomalous isotopic values) within the series of isotopes.

Furthermore, we provided a possible explanation for the geographic and climatic factors (i.e., catchment descriptors) influencing the several regressions and finally we made some considerations concerning their applicability in the context of mountainous areas such as the northern Italian Apennines.

#### **2. Overview of the Climatic, Geomorphological, and Hydrogeological Features of the Study Area**

The study area is located in the northern Italian Apennines and belongs to the Emilia– Romagna Region (Figure 1). It includes nine catchments between the Trebbia River and the Savio River, with river gauges (in which the samples were collected) located close to the foothills of the mountain chain. The area has an overall extension of 6261 km<sup>2</sup> . Maximum altitudes are in the southern sectors, where the main watershed divide lies (with main peaks showing elevations higher than 2000 m a.s.l., such as Mt. Cimone with its 2165 m a.s.l.) is the southern border of the Emilia–Romagna Region. Elevation decreases towards the NE direction to approximately 40 m a.s.l. at the Savio River gauge.

All the nine rivers originate from the main watershed divide and flow towards the NE. Six rivers (namely, Trebbia, Nure, Taro, Enza, Secchia, and Panaro) are tributaries of the Po River while the other three (Reno, Lamone, and Savio) enter the Adriatic Sea. Catchment areas are between 193 km<sup>2</sup> (Lamone) and 1300 km<sup>2</sup> (Secchia), while flow lengths range from 28 km (Enza) to 85.2 km (Secchia). Mean annual discharges during the period 2006–2016 are included between 8.4 m<sup>3</sup> s −1 (Savio) and 30.4 m<sup>3</sup> s −1 (Secchia).

From a hydrogeological point of view, a report [20] grouped bedrock outcrops in the aforementioned catchments into six main classes (or hydrogeological complexes, namely: clay, marl, flysch, foreland flysch, ophiolite, and limestone). Those composed of poorly permeable or impermeable materials (clay, marl, and flysch hydrogeological complexes; see Figure 1) are the most represented in terms of areal coverage, leading to a runoff response of rivers that closely follows the rainfall distribution during the year (pluvial discharge regime). Rivers originating from the most elevated parts of the main watershed divide (Secchia, Panaro) are characterised by a nival–pluvial discharge regime as they are influenced by the melting of snow cover accumulated during the winter months in the upper parts of their catchments [21].

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 4 of 20

**Figure 1.** Sketch map of the area (modified after [7]) with locations where sampling has been carried out by previous studies for precipitation (rain gauges with letters a to d), surficial water (rivers numbered 1 to 9), groundwater from springs (Greek letters α and β), and groundwater from wells (located in the alluvial fans with capital letters A to D). Hydrogeological complexes are reported following [20]; GC: clay; GM: marl; GF: flysch; GFF: foreland flysch; GL: limestone; GO: ophiolite. For further details, see Table 1. **Table 1.** Main features of the sampling points from which isotopic data were derived. For the **Figure 1.** Sketch map of the area (modified after [7]) with locations where sampling has been carried out by previous studies for precipitation (rain gauges with letters a to d), surficial water (rivers numbered 1 to 9), groundwater from springs (Greek letters α and β), and groundwater from wells (located in the alluvial fans with capital letters A to D). Hydrogeological complexes are reported following [20]; GC: clay; GM: marl; GF: flysch; GFF: foreland flysch; GL: limestone; GO: ophiolite. For further details, see Table 1.

**Location Type Code Number of Samples Timing of Sampling Length of Time Reference s Table 1.** Main features of the sampling points from which isotopic data were derived. For the corresponding map locations, readers are referred to Figure 1.


corresponding map locations, readers are referred to Figure 1.


**Table 1.** *Cont.*

In the vicinity of the foothills (therefore close to the corresponding river gauges), several wells drilled in the alluvial fans of the Trebbia, Taro, Enza, and Secchia rivers continuously pump groundwater for both agricultural and drinking purposes. As previously highlighted by [11], by exploiting water stable isotopes, wells pumping water from confined aquifers in Trebbia and Taro alluvial fans are also likely to be recharged by zenithal precipitation infiltrating through gravels and sands that outcrop at the foothills of the northern Apennines (i.e., apical part of the alluvial fans). On the contrary, an important quota of recharge also occurs from streambed dispersion (focused on the apical part of their alluvial fans, see [22]) seems to affect wells located in the alluvial fans of the Enza and Secchia rivers. Two groups of low-yield springs from the Secchia River catchment (namely, Pietra di Bismantova and Montecagno) were also considered. These springs should be considered as representative of the common ones in the northern Italian Apennines, whose discharges are closely related to the rainfall pattern while outflows are strongly reduced (often in the order of 1 L·s –1 or less) at the end of the summer periods (shallow groundwater flow paths; for more details see [20,23,24]).

From the climatic point of view, and as already reported in [25], the mean annual rainfall distribution during the period 1990–2015 exceeds 2200 mm/y near the main watershed divide and progressively decreases to about 900 mm/y in the foothills. The rainfall distribution during the year is characterised by a marked minimum in the summer season and two maxima during autumn (the main one) and spring. Close to the main watershed divide, the cumulative annual snow cover can reach 2–3 m at the end of the winter season [21]. Potential evapotranspiration ranges from about 500 up to 650 mm/y in the lowlands and is mainly active during the summer months [25].

#### **3. Methodology**

The methodology used in this study consists of five steps involving data inspection and statistical comparison on datasets, including rainfall (4 rain gauges, namely: Parma, Lodesana, Langhirano, Berceto), surface water (9 rivers, namely: Trebbia, Nure, Taro, Enza, Secchia, Panaro, Reno, Lamone, Savio) groundwater from wells (aggregated isotopic values from wells belonging to 4 alluvial fans, namely: Trebbia—4 wells, Taro—5 wells, Enza—3 wells, and Secchia—5 wells), and groundwater from springs (aggregated isotopic

values from springs belonging to 2 areas from the Secchia catchment, namely: Pietra di Bismantova and Montecagno). Firstly, a check on the assumption of stationary behaviour of each series of stable isotopes was carried out by means of conventional statistical tests coupled with inspection of standardised residuals.

Secondly, slopes and intercepts from δ <sup>2</sup>H–δ <sup>18</sup>O alignments were obtained by using 5 different regression approaches. Moreover, we further considered 2 different regressions applied to δ-values from rain gauges and rivers that had been preliminary weighted to the corresponding quota of precipitation and discharge, respectively. Thirdly, slopes and intercepts were visually inspected by means of heat maps to identify discrepancies among the several regression methods. Fourthly, slopes and intercepts were compared through bivariate (correlation matrices) and multivariate analyses (hierarchic clustering, i.e., dendrograms) to identify linear correlations and similarities. Fifthly, and with reference to the only surface water, we made a comparison between differences in slopes and intercepts with some selected catchments to verify linear or nonlinear correlations among these variables.

For convenience (see Figure 1 for location of the sampling points and Table 1 for further details), we report below rain gauges signified by letters (from a to d: "a"—Parma, "b"—Lodesana, "c"—Langhirano, "d"—Berceto); surface water locations as numbers (from 1 to 9: "1"—Trebbia, "2"—Nure, "3"—Taro, "4"—Enza, "5"—Secchia, "6"—Panaro, "7"—Reno, "8"—Lamone, "9"—Savio); groundwater from wells as capital letters ("A"—Trebbia, "B"—Taro, "C"—Enza, "D"—Secchia); and groundwater from springs as Greek letters ("α"—Pietra di Bismantova, "β"—Montecagno).

#### *3.1. Isotopic Datasets*

The dataset from 4 rain gauges and 9 rivers consists of monthly isotopic data while 17 water wells were characterised by grabbed four-monthly samples. All the data are derived from [7,16]. With reference to rain gauges, isotopic datasets lasted over the period from January 2003 to December 2006. Isotopic data from surface water and groundwater covered the period from January 2004 to December 2007. Further monthly, two-monthly, and three-monthly isotopic data from 2 groups of nearby low-yield springs were considered. These data were published in [17,18] and included the period January 2014 to December 2016.

The final dataset consists of 553 isotopic values, of which 91 are from precipitation, 264 from surface water, 145 from groundwater collected in wells, and 53 from groundwater collected in springs. Precipitation and river discharge that were used for further weighting procedures (see Section 3.3, "Linear Regression Types") come from [26].

As reported in the previous works of [7,11,16], the isotopic analyses were carried out by using isotope ratio mass spectrometry (IRMS) while instrument precision (1σ) was on the order of ±0.05‰ for δ <sup>18</sup>O and <sup>±</sup>0.7‰ for <sup>δ</sup> <sup>2</sup>H. With reference to groundwater from springs [17,18], the corresponding isotopic analyses were carried out by mixed technique involving IRMS for δ <sup>18</sup>O and cavity ring-down spectroscopy (CRDS) for δ <sup>2</sup>H. Instrument precision (1σ) was assessed as ±0.1‰ for δ <sup>18</sup>O and <sup>±</sup>1.0‰ for <sup>δ</sup> <sup>2</sup>H.

#### *3.2. Verifying the Stationary Behaviour of Isotopic Data Series*

As anticipated in the introduction, three of the five linear regression approaches are considered in this study, namely, ordinary least squares (OLS), reduced major axis (RMA) and major axis (MA), which are based on the assumption that a series of stable isotopes are stationary [27,28]. This means that statistical properties as mean and variance remain constant along each δ <sup>18</sup>O–δ <sup>2</sup>H alignment, i.e., modelling errors are normally distributed and homoscedastic. Moreover, the closest pairs of δ <sup>18</sup>O–δ <sup>2</sup>H must not be affected by autocorrelation phenomena as the modelling errors (i.e., residuals) from the regressions must be independent [27].

The presence of outliers or heteroscedasticity (i.e., modelling errors have not the same variance over the alignment) or autocorrelation lead the assumption of stationarity to be violated, thus slopes and intercepts from the abovementioned regression may not

be meaningful. In this work, we exploited conventional statistical tests for verifying multivariate normality (i.e., presence of outliers inducing non-normality; Doornik–Hansen test [28,29]) heteroscedasticity (Breusch–Pagan test [30]), and autocorrelation (Durbin– Watson [31]).

All of the abovementioned tests are based on a comparison of the corresponding statistics' *p*-value results with a threshold value (level of significance α set at 0.01) to decide whether the null hypotheses have to be rejected (*p* < 0.01) or not (*p* > 0.01). The following null hypotheses were selected: multivariate normality (Doornik–Hansen), homoscedasticity (Breusch–Pagan), no autocorrelation (at a lag of 1) in the residuals (Durbin–Watson). With reference to the Durbin–Watson test, it must be specified that values are included between 0 and 4. Values close to 0 indicate almost total positive autocorrelation while results in the proximity of 4 indicate a total negative autocorrelation. Values between 1.5 and 2.5 show there is no autocorrelation in the data.

As suggested by [27], in the case of rejection of the null hypothesis of multivariate normality or homoscedasticity, we further verified the standardised residuals for identifying outliers (i.e., those values for which standardised residuals fall outside the interval from −4 to 4). Moreover and always following [27], if residuals were found to be serially correlated, we made autocorrelation functions of the standardised residuals to further confirm the lag length of autocorrelation.

#### *3.3. Linear Regression Types*

In this work, 5 types of linear regressions were tested, namely, ordinary least squares (OLS), reduced major axis (RMA), major axis (MA), robust (R), and Prais–Winsten (PW). To avoid excessive mathematical details, we provide a cursory examination of the methods and the reader is referred to more specific literature [14,15,27,32]. For convenience, we report all the corresponding equations by replacing δ <sup>18</sup>O with x and δ <sup>2</sup>H with y, while n is the number of samples and r is the Pearson correlation coefficient.

The OLS regression assumes that the x values are fixed (i.e., it is commonly used when *x* values have very few associated errors) and finds the line that minimises the squared errors in the y values. The slope of the linear regression (slopeOLS) is calculated as follows:

$$slope\_{OLS} = r \times \frac{sd\_y}{sd\_x} \tag{1}$$

where *sd* represents standard deviations calculated for *x* variables (*sdx*) and *y* variables (*sdy*). Unlike OLS, RMA and MA try to minimise both the *x* and the *y* errors [33]. In the case of RMA, the corresponding slopeRMA can be obtained with:

$$slope\_{RMS} = sign[r] \times \frac{sd\_y}{sd\_x} \tag{2}$$

where *sign[r]* is the algebraic sign of the Pearson coefficient. The slopeMA is calculated by:

$$slope\_{MA} = -A + \frac{\sqrt{r^2 + A^2}}{r} \tag{3}$$

where *A* can be obtained as:

$$A = 0.5 \times \left(\frac{sd\_x}{sd\_y} - \frac{sd\_y}{sd\_x}\right) \tag{4}$$

As anticipated in introduction, the PW regression [34] has never been used for developing δ <sup>18</sup>O–δ <sup>2</sup>H alignments, as series of isotopic data have always been considered to the present as time-invariant (i.e., stationary). Recently, [19,35] highlighted that multiannual series of such isotopic data may be affected by nonstationary processes (such as, for example, trends in the means or presence of far-off values). In this case (nonstationary multiannual series of δ <sup>18</sup>O–δ <sup>2</sup>H pairs), the use of common regression methods such as OLS, RMA, and MA could induce residuals to be larger and characterised by stronger serial correlations. PW is commonly used for data with serially correlated residuals of the estimates [36]. As a matter of fact, this approach takes into account AR1 (i.e., autoregression of the first order) serial correlation of the errors in a linear regression model. The procedure recursively estimates the coefficients and the error autocorrelation of the model until sufficient convergence is reached. All the estimates are obtained by the abovementioned OLS.

As in the case of PW approach, the R method has also never been tested for δ <sup>18</sup>O–δ <sup>2</sup>H regressions. This method is an advanced Model I (in which x is always the independent variable) regression which is less sensitive to outliers than OLS estimates. Having less restrictive assumptions, R is recognised to provide much better regression coefficient estimates than OLS when outliers are present in the data. In particular, this approach has proven to be successful in the case of "almost" normally distributed errors but with some far-off values. This happens as outliers usually violate the assumption of normally distributed residual in OLS method. The algorithm is "least trimmed squares" reported by [37], in which the method consists of finding that subset of x–y pairs whose deletion from the entire dataset would lead to the regression having the smallest residual sum of squares. As in the case of PW approach, estimates of each subset are calculated owing to the OLS method. It must be added that, depending on the size and number of outliers, R regression conducts its own residual analysis and downweight or even these x–y pairs; this fact deserves an accurate inspection of the outliers made by the operator prior to any removal in order to decide whether these x–y pairs have to be considered or not.

For all 5 different regression approaches, the corresponding intercept is obtained with the following:

$$intercept = \frac{\sum\_{i=1}^{n} y\_i}{n} - slope \frac{\sum\_{i=1}^{n} x\_i}{N} \tag{5}$$

In all the abovementioned linear regressions, each observation has an equal influence of the orientation of the fitted line. As a matter of fact, it is well recognised that some isotopic data may be more important than others as related to a higher amount of water (for example, a flood in the case of a river or a high discharge event of a spring or high rainfall amount during a storm event). In this case, greater influence in the regression should be given to these isotopic data. In order to also take this effect into account, OLS were applied to isotopic datasets that had been previously weighted on the corresponding monthly amount of precipitation (rainwater) and discharge (freshwater and river water) by means of two different methods, i.e., the classical one that simply involves multiplying each *y<sup>i</sup>* by the water amount (see [28] for further details; hereafter called W) and as reported in [38] (see [7] for the formulation; hereafter called B).

#### *3.4. Comparison among Regressions*

Initially, slopes and intercepts from all the regressions were compared by means of heat maps. The heat maps are matrices of fixed cell size showing the magnitude of difference among values with a selected binary colour ramp (in our case from red to green, respectively, indicating the lowest value and highest value within the isotopic dataset), in which the colour intensities provide visual cues to the reader about discrepancies between the data. The goal was to verify any presence of clusters to be further investigated by means of bivariate and hierarchical cluster analyses.

As a bivariate analysis, we carried out scatterplot matrices to determine if linear correlation between multiple variables (slopes and intercepts obtained by the different regression approaches) were present or not. Tests were carried out highlighting the level of significance, which was set as *p* < 0.01.

Furthermore, a hierarchical cluster analysis was performed to identify similarities among the series of slopes and intercepts from the whole datasets. Clustering was done according to the unweighted pair–group average (or centroid) method, in which each group consisted of slopes (or intercepts) from a determined regression approach. The method was based on a step-by-step procedure in which series of slopes (or intercepts) were grouped into branched clusters (dendrogram) based on their similarities to one another. As a result, the two most similar series of slopes (or intercepts) were selected and linked based on the smallest average distance among the values of all slopes (or intercepts). Progressively more dissimilar series were linked at greater distances; in the end, they all were joined to one single cluster. The cophenetic coefficient was used as a measure of similarity between each pair of clusters; more than 2 time series being analysed, the dendrogram was supported by a cophenetic distance matrix. Further details on this method can be found in [28].

#### *3.5. A Focus on the Differences among Regression Approach from River Water: Comparison with Catchment Characteristics*

As suggested by [15], we investigated whether some selected catchment characteristics (also called descriptors) could have affected differences among values of slopes and intercepts as obtained by the different approaches reported in Section 3.3. In order to make all slopes and intercepts comparable, we followed the approach proposed by [15] that consisted of prior computed differences in the slopes (as slopeOLS–slopeRMA/MA/R/PW/W/B) and intercepts (as interceptOLS–interceptRMA/MA/R/PW/W/B). Then, and following again the procedure reported in [15], we applied the Spearman ranking correlation matrix in which the abovementioned differences in slopes and intercepts were compared with 9 catchment characteristics. It should be added that this approach is a nonparametric measure of rank correlation that provides a statistical dependence between the rankings of two variables at a time. Unlike the Pearson coefficient, the Spearman ranking assesses how well the relationship between two variables can be described using a monotonic function, even if their relationship is not linear [27]. In particular, several authors have highlighted that many hydrological processing occurring at both slope and catchment scales are nonlinear (see for instance [38,39]) and such behaviour was in turn seen in some descriptors calculated by means of time series of stable isotopes of water [7,40–42].

In order to take into account the linearity among the variables, we also considered the linear correlation by providing the Pearson correlation coefficients. Here, we recall that Pearson and Spearman matrices reflect the magnitude of similarity among the parameters by means of r (the Pearson correlation coefficient described in Section 3.3) and r<sup>s</sup> coefficients, respectively. Both correlation coefficients (r and rs) describe the strength and direction between the two variables and return a closer value to 1 (or −1) when the two different datasets have a strong positive (or negative) relationship. The significance probability (*p*-value) for both the Spearman and Pearson correlation coefficients calculated in this study was set at 0.01, meaning that *p*-values lower than 0.01 represented statistically significant relationships. Readers are referred to [28] for further details on statistical formulations.

The 9 catchment characteristics (or descriptors, see Table 2 and Supplementary Materials Table S1 for further details) were those already considered by [7], namely: catchment area (A); elevation (H); precipitation (P); flow length (F); specific mean annual runoff (q); specific river runoff exceeded for 95% of the observation period (q95; this is a well-known low flow index that is used worldwide for the regionalisation procedure and can be estimated even from a relatively short time series of daily runoffs [43]); and the young water fraction (Fyw proposed by [42]; this is considered the percentage proportion of catchment outflow younger than approximately 2–3 months and was estimated from the amplitudes of seasonal cycles of stable water isotopes in precipitation and stream flow that had been already calculated in [7]).

It should be noted that 4 descriptors (P, q, q95, Fyw) were obtained by processing daily precipitation (42 rain gauges homogeneously distributed over the study area for P), discharges (for both q and q95), and water isotopes time series (for Fyw) lasting over the same time period. Moreover, flow length (F) was derived from a 5 × 5 m gridded digital terrain model created by the digitalisation and linear interpolation of contour lines represented in the regional topography map at a scale of 1:5000.


**Table 2.** The 9 catchment characteristics included in the analysis. For further details on the catchment characteristics from a single catchment, see Table S1 in Supplementary Materials.

#### **4. Results**

#### *4.1. Stationary Behaviour of Isotopic Data Series*

Table 3 summarises the results from the three statistical tests (Doornik–Hansen for multivariate normality, Breusch–Pagan for homoscedasticity, and Durbin–Watson for autocorrelation) used for assessing the compliance with the stationary assumption. Isotopic series from rivers were those mainly affected by problems of non-normal behaviour (rivers "5,6,9") and autocorrelation at a lag of 1 (rivers "1,4,6,7"). The latter were positive (we recall here that values closer to 0 identify positive autocorrelation phenomena) and more intense in the case of rivers "4,6". Moreover, the Breusch–Pagan test suggested residual homoscedasticity for river "8". By considering the plots of standardised residuals (see Figure S1 in Supplementary Materials), the presence of outliers was further confirmed for rivers "5,6,9" (river "5", 3 outliers; river "6", 4 outliers; river "9", 1 outlier) as well as the increase of variance of standardised residuals along estimates for river "1" (i.e., heteroscedasticity). Autocorrelation functions carried on standardised residuals (see Figure S2 in Supplementary Materials) allowed for demonstrating the presence of serial correlations, although with different lag lengths (river "1", 2 lags; river "4", 2 lags; river "6", 3 lags; river "7", 2 lags).

#### *4.2. Slopes and Intercepts*

The δ <sup>18</sup>O–δ <sup>2</sup>H relationships are summarised in in Supplementary Materials containing slopes (a; see Table S2), intercepts (b; see Table S3), standard deviation of the estimates (c; see Table S4), standard deviations of the estimates and coefficient of determinations (d: see Table S5) coefficient of determinations R<sup>2</sup> . By viewing all the results reported in the form of heat maps (see Figure S3 in Supplementary Materials), the substantial invariance of slopes (from 6.9 to 13.1) and intercepts (from 7.2 to 8.4) from rain gauges located at lower altitudes (a, b, c), with high performance of the regression (R<sup>2</sup> always close to 0.99) was noticed. These results are in agreement with the GMWL (we recall that this line is characterised by slope and intercept equal to 8.0 and 10.0, respectively; see Figure S4 in Supplementary Materials) with no evidence of outliers. When the two weighting approaches were taken into account, no changes among the unweighted values of slope and intercept were found with the exception of intercepts in the rain gauge "a" (intercepts remarkably lower in the case of weighting procedures in the order of +4.0).


**Table 3.** Results from the three statistical tests aimed at verifying the compliance with the stationary assumption, namely: Doornik–Hansen (multivariate normality), Breusch–Pagan (homoscedasticity), and Durbin–Watson (autocorrelation). \* Null hypotheses rejected as *p* < 0.01.

With reference to the rain gauge "d", i.e., that located near the main watershed divide, the R approach provided remarkably higher values of both slope (8.1) and intercept (11.5) than those obtained with the other regression approach (we recall that all values of intercepts from "d" were negatives). It must be highlighted that the standard deviations of the estimates are slightly higher than those obtained with the other regressions (see Table S2 in Supplementary Materials).

By considering the surface water, the RMA and MA approaches almost provided slightly higher values of slopes and intercepts (up to 13.0 and 55.2, respectively, in the case of river "5"). In the case of rivers "1,2,3,4,7,8", values of slopes are in the range of those obtained by weighting procedures. On the contrary, intercepts showed a larger variability among the regression methods investigated. It should be highlighted that in the case of rivers "5,6,9" the values of slopes remarkably varied as well, in particular if MA and R were used. As in the case of the abovementioned rain gauge "d", the δ <sup>18</sup>O–δ <sup>2</sup>H alignments from "5,6,9" were characterised by the lowest values of R<sup>2</sup> and the larger values of standard deviations of the estimates (see Table S4 in Supplementary Materials).

It must be noted that the discrepancies reported for these points (i.e., "5,6,9") affected water with the presence of several outliers and/or serial correlations of residuals (see Figures S2 and S3 in Supplementary Materials and Section 4.1), which violated the stationary assumption.

Akin to the cases of rain gauges and rivers, RMA and MA approaches carried out on groundwater from wells and springs were characterised by larger values of both slopes and intercepts than in the case of OLS. With the exception of groundwater from "C" and "D", the R and PW approaches induced larger variations in both slopes and intercepts, which

were particularly marked in the case of groundwater from β (i.e., water whose δ <sup>18</sup>O–δ <sup>2</sup>H alignment was characterised by low values of R<sup>2</sup> and large standard deviations of the estimates, see Table S2 in Supplementary Materials).

In Table 4, the matrix reporting correlation coefficients between pairs of slopes indicated that the largest degree of association was found (*p* < 0.01) between OLS–PW and W–B. High values of correlation (with slightly lower degree of association but still *p* < 0.01) also characterised the following relations: OLS–W, OLS–RMA, OLS–B, RMA–PW, RMA–MA, and PW–W. A significant degree of association (*p* < 0.01) was also found for PW–B. It should be highlighted that in several cases regarding R and MA, the degree of associations was very low and, in some cases, even negative (i.e., an increase in the value of slope obtained with a regression corresponds to a decrease in the series obtained with RMA).

**Table 4.** Correlation matrix reporting associations among the slopes from different regression approaches considered in this study (namely: OLS, RMA, MA, R, W, B). Progressively darker green colour is associated with a higher correlation coefficient. \* Significant as *p* < 0.01.


By considering the intercepts (Table 5), the degree of associations already highlighted for slopes was further confirmed with the exception of OLS–W and OLS–B (here not significant as *p* > 0.01). It should be noted that almost all correlations were slightly lower than the corresponding ones from the slopes.

**Table 5.** Correlation matrix reporting associations among the intercepts from different regression approaches considered in this study (namely: OLS, RMA, MA, R, W, B). Progressively darker green colour is associated with a higher correlation coefficient. \* Significant as *p* < 0.01.


The hierarchic cluster analysis (see Figure S5 in Supplementary Materials) among the slope series from the several regression approaches (reported as different branches composing the dendrogram) demonstrated that MA was associated with none of the other regression methods, while two main group of pairs were clearly separated: the first is represented by OLS–PW while the second by W–B. The aforementioned first and second group were associated with each other while longer branches further linked them to R and RMA.

With reference to intercepts, the dendrogram confirmed the nonassociation of MA with the other regression approaches. Moreover, the two closest series were still those of OLS and PW, which were in turn associated to W. Contrary to the case of slopes, B series was strictly associated to RMA.

#### *4.3. Comparison between the Differences in the Slopes and Intercepts with Catchment Characteristics and Statistics*

By taking into account only δ <sup>18</sup>O–δ <sup>2</sup>H regressions from surface water, the Pearson rank correlation matrix (we recall that Pearson assesses linear relationships) comparing differences in slopes with catchment characteristics (see Table 6) did not provide significant (*p* < 0.01) correlations. If the Spearman rank correlation matrix (i.e., assessment of nonlinear relationships) was considered, we found positive and significant (*p* < 0.01) correlations with Fyw (SlopeOLS–RMA and SlopeOLS–MA) while negative ones with Hmin (SlopeOLS–W, SlopeOLS–B).

**Table 6.** Matrix of the Pearson (in grey) and Spearman (in green) rank correlation (values as r and rs, respectively; r value evidenced in grey while rs in green) between the differences in the slopes and the selected catchment characteristics considered for the 9 rivers. Relationship between differences in the slopes and coefficient of determination R<sup>2</sup> from δ <sup>18</sup>O–δ <sup>2</sup>H linear regressions are also reported. R <sup>2</sup> values from regressions were calculated starting from signed values of differences in slopes. \* Significant as *p* < 0.01.


In the case of the Pearson rank correlation matrix applied to differences in intercepts, we did not find significant correlations (see Table 7). On the contrary, and as in the case of slopes, significant (*p* < 0.01) nonlinear relationships were found for Fyw (SlopeOLS–RMA and SlopeOLS–MA) and Hmin (negative correlation for SlopeOLS–B).

**Table 7.** Matrix of the Pearson (in grey) and Spearman (in green) rank correlation (values as r and rs, respectively; r value evidenced in grey while rs in green) between the differences in the intercepts and the selected catchment characteristics considered for the 9 rivers. Relationship between differences in the slopes and coefficient of determination R<sup>2</sup> from δ <sup>18</sup>O–δ <sup>2</sup>H linear regressions are also reported. R <sup>2</sup> values from regressions were calculated starting from signed values of differences in intercepts. \* Significant as *p* < 0.01.


In all cases, the largest statistical performance (again significant as *p* < 0.01) was found for coefficient of determination R<sup>2</sup> from δ <sup>18</sup>O–δ <sup>2</sup>H linear regressions (SlopeOLS–RMA and SlopeOLS–MA; InterceptOLS–RMA and InterceptOLS–MA).

#### **5. Discussion**

We did not find significant discrepancies in the slopes and intercepts computed by the different regression methods in the case of precipitation data. On the contrary, marked variations were detected in the case of river water and groundwater (both from springs and wells in lowland aquifers) obtained using specific methods. Among others, such discrepancies were somehow reduced in the case of OLS, RMA, and PW. Because of different values of river and spring discharges and the corresponding changes in isotopic content of water during the year, weighting procedures (W, B) were characterised by diverse values of slopes and intercepts rather than the aforementioned OLS, RMA, and PW. Moreover, and as highlighted by both heat maps (Figure S3 in Supplementary Materials) and correlation matrices (Tables 4 and 5) and dendrograms (Figure S5 in Supplementary Materials), slopes and intercepts from the MA and R approaches were not comparable to others (nonsignificant statistical associations). Regardless of water type, the aforementioned discrepancies were promoted when δ <sup>18</sup>O–δ <sup>2</sup>H regressions were characterised by weak statistical performances (low values of R<sup>2</sup> and larger values of standard deviations of the estimates). With reference to rivers, the weak statistical performances were linked to the presence of outliers and/or serial correlation of the residuals violating the stationary assumption of OLS, MA, and RMA approaches.

The investigation carried out on the data solely from rivers highlighted that the magnitude of the differences in the slopes and intercepts was related in all cases (with the exception of R and PW) to the coefficient of determination R<sup>2</sup> characterising δ <sup>18</sup>O–δ <sup>2</sup>H linear regressions. The largest values of Pearson coefficients (see Tables 6 and 7) led us to consider R<sup>2</sup> as the main causal factor for such differences in slopes and intercepts.

In particular, the larger the correlations between δ <sup>18</sup>O and δ <sup>2</sup>H, the smaller the differences among slopes and intercepts detected by RMA, MA, W, and B within the specific sampling point (river, well, or spring). This is in agreement with the results reported by [15] and corroborated the hypothesis that statistical performance of the regression was the main driver of these slope and intercept variations. In any case, despite finding highly statistical significance with R<sup>2</sup> in our investigated dataset, no relations between differences of slopes (and intercepts) and the ranges in δ <sup>18</sup>O (and δ <sup>2</sup>H) along with the number of samples composing the dataset were noticed, thus indicating that extreme values of δ <sup>18</sup>O (and δ <sup>2</sup>H) were not significant causal factors.

With reference to RMA and MA, the Spearman rank correlation matrices involving differences in slopes and intercepts and catchment descriptors allowed us to find a significant nonlinear association with Fyw (we recall here that Fyw is the percentage of water younger than 2–3 months). In both cases (RMA and MA) the association (reported also as plots in Figure 2) indicated that the magnitude of the differences in the slopes and in the intercepts decreased along with the quota of young water.

This means that rivers showing low values of Fyw are likely to be more affected by differences in slopes and intercepts computed by different regression approaches. By examining the plots reported in Figure 2, it can be evidenced that nonlinearity is driven by two catchments (namely, the Secchia River "5" and Panaro River "6"). As already anticipated in Section 2, these two rivers ("5,6") were the only ones characterised by nival– pluvial discharges due to the melting of the snow cover in the upper part of the catchments during the spring months. Moreover, [7] stated that there were evidences of sublimation in several water samples collected in rivers "5,6" from January 2017 to April 2017. This was further confirmed by the remarkable number of snowfall events that occurred between December 2013 and April 2014 over the highest part of the catchments "5,6". Such events have allowed the consequent snowpack development alternating with partial snowmelt for a snow water equivalent higher than 600 mm [21].

**Figure 2.** Differences in slopes (a: OLS–RMA: c: OLS–MA) and intercepts (b: OLS–RMA; d: OLS– MA) for OLS–RMA and OLS–MA. Values of differences in slopes and intercepts (y-axes) are reported in modulus form. Codes for rivers (from 1 to 9) are also reported (for further details on river codes, see Table 1). **Figure 2.** Differences in slopes (a: OLS–RMA: c: OLS–MA) and intercepts (b: OLS–RMA; d: OLS–MA) for OLS–RMA and OLS–MA. Values of differences in slopes and intercepts (y-axes) are reported in modulus form. Codes for rivers (from 1 to 9) are also reported (for further details on river codes, see Table 1).

This means that rivers showing low values of Fyw are likely to be more affected by differences in slopes and intercepts computed by different regression approaches. By examining the plots reported in Figure 2, it can be evidenced that nonlinearity is driven by two catchments (namely, the Secchia River "5" and Panaro River "6"). As already anticipated in Section 2, these two rivers ("5,6") were the only ones characterised by nival–pluvial discharges due to the melting of the snow cover in the upper part of the catchments during the spring months. Moreover, [7] stated that there were evidences of sublimation in several water samples collected in rivers "5,6" from January 2017 to April 2017. This was further confirmed by the remarkable number of snowfall events that occurred between December 2013 and April 2014 over the highest part of the catchments "5,6". Such events have allowed the consequent snowpack development alternating with partial snowmelt for a snow water equivalent higher than 600 mm [21]. There, we recall that sublimation occurring during sunny days can modify the former isotopic composition of the superficial snow layers, allowing the release of a vapour phase from the solid skeleton to the atmosphere. In this case, the final snow cover does not preserve the isotopic composition of the original snowfall from which it was derived [44,45], a fact that also led differences in slopes and intercepts from δ <sup>18</sup>O–δ <sup>2</sup>H regressions to be enhanced. In detail, sublimation acting on a snow cover can lead to an enrichment of heaviest isotopes (such as <sup>18</sup>O and <sup>2</sup>H) within the solid skeleton and can induce a similar δ <sup>18</sup>O–δ <sup>2</sup>H pattern of that charactering the residual liquid subjected to evaporation (slope decrease of the δ <sup>18</sup>O–δ <sup>2</sup>H alignments for snowpack samples if compared to the water meteoric line (see for field and experimental studies: [46,47])). In particular, the slope decrease can be much more intense if the only late-season snowpack samples are considered (a value of 3.7 was found by [48]).

There, we recall that sublimation occurring during sunny days can modify the former isotopic composition of the superficial snow layers, allowing the release of a vapour phase from the solid skeleton to the atmosphere. In this case, the final snow cover does not preserve the isotopic composition of the original snowfall from which it was derived [44,45], a fact that also led differences in slopes and intercepts from δ18O–δ2H regressions to be enhanced. In detail, sublimation acting on a snow cover can lead to an enrichment of heaviest isotopes (such as 18O and 2H) within the solid skeleton and can induce a similar δ18O– δ2H pattern of that charactering the residual liquid subjected to evaporation (slope decrease of the δ18O–δ2H alignments for snowpack samples if compared to the water mete-In case the two rivers "5,6" are removed from the analysis, it is still possible to confirm such alignments, although linear, between y Fyw and differences in slope and intercept pairs. In this sense, such relations, still identifying an inverse association between differences in the slopes (and in the intercepts) and quotas of young water, may also be related to other hydrological processes taking place at the catchment scale. As already pointed out by [7], by checking both slopes (river water showed slightly lower values than those characterising rainwater; see Figure S6 in Supplementary Materials) and intercepts (that were negative compared to those from rainwater), all the river water considered underwent evaporation/evapotranspiration processes prior to their infiltration towards the aquifer.

oric line (see for field and experimental studies: [46,47])). In particular, the slope decrease can be much more intense if the only late-season snowpack samples are considered (a value of 3.7 was found by [48]). In case the two rivers "5,6" are removed from the analysis, it is still possible to confirm such alignments, although linear, between y Fyw and differences in slope and intercept pairs. In this sense, such relations, still identifying an inverse association between Albeit to a lesser extent, these variations also affected groundwater from wells (which were also fed by streambed dispersion and therefore by water isotopes already modifiedin the river water; see [11] and low-yield springs (these potentially characterised by preinfiltrative modification as slopes from δ <sup>18</sup>O–δ <sup>2</sup>H alignments were slightly lower than those obtained from precipitation water; see Figure S7 in Supplementary Materials)). On the contrary, the nonvariability of slopes and intercepts observed in the different δ <sup>18</sup>O–δ <sup>2</sup>H

differences in the slopes (and in the intercepts) and quotas of young water, may also be

alignments from precipitation was somehow expected, as these waters were unlikely to be affected by evaporative/sublimation processes once they had entered the rain gauges.

With reference to the effective role of young water fraction Fyw in influencing the differences in intercepts and slopes, we believe that further efforts have to be made for such catchments from the hilly part of the northern Italian Apennines and dominated by higher quotas of young water (i.e., Fyw > 25%; not analysed in this study for lack of isotopic data). The latter are characterised by intermittent discharge and wide outcrops of low permeable soils and bedrocks (prevailing clayey and marls materials; G<sup>C</sup> and G<sup>M</sup> in Figure 1). Further investigations should be isotopic-based in order to verify also the role of pre-infiltrative evaporation in isotopic deviation and, above all, in the change of slopes and intercepts.

Following [7], these processes were likely to be promoted in the clay-rich bedrock, where water molecules composing the soil moisture were slowed in percolation and thus kinetic fractionation processes were enhanced.

However, we can provide some preliminary recommendations for use of the different regression approaches for the four water types (precipitation, surface water, groundwater from wells, and low-yield springs) from the northern Italian Apennines:


#### **6. Conclusions**

We presented the comparison of five different regression approaches applied to δ <sup>18</sup>O– δ <sup>2</sup>H data from four different water types collected in the northern Italian Apennines. We found that all the tested approaches converged towards similar values of slopes and intercepts for only stable water isotopes from precipitation. Conversely, differences in slopes and intercepts from surface water and groundwater (collected from wells and low-yield springs) were often significant and related to the robustness of the regressions (i.e., standard deviations of the estimates) and their sensitiveness to outliers and autocorrelation. Moreover, and with reference to surface water, we found evidence of a relationship between young water fraction and the magnitudes in differences of slopes and intercepts, suggesting the control of kinetic fractionation processes (mainly related to sublimation acting on snow cover and, secondary, to active pre-infiltrative evaporation and evapotranspiration processes) on such discrepancies. These results allowed us to provide some recommendations for hydrological and hydrogeological studies involving δ <sup>18</sup>O–δ <sup>2</sup>H from the abovementioned water types collected in the northern Italian Apennines. Firstly, as no discrepancies were noticed between slopes and intercepts from all the methods applied to precipitation, the OLS approach is preferred. Secondly, and with reference to the other water types (surface water and groundwater from wells and springs), we warmly suggest carrying out conventional statistical tests coupled with inspection of standardised residuals for a preliminary check on the presence of outliers and autocorrelation phenomena. In the case of managing outliers, the MA and R approaches should be avoided as they are more sensitive to the statistical performance of the regressions and often provide unrealistic values of both slopes and intercepts. Thirdly, for surface water and groundwater, the OLS and PW approaches still showed the highest degree of robustness and produced the closest values of slopes and intercepts, thus resulting as the methods preferable for δ <sup>18</sup>O–δ <sup>2</sup>H regressions. PW would be more reliable in the presence of serial correlations of the residuals (which, in our case, often affected surface water). In the case of managing outliers, the possibility of removing them will have to be considered (as an example in the case of δ <sup>18</sup>O–δ <sup>2</sup>H pairs from marked low-flow periods).

Lastly, despite the presence of marked differences in the amounts of rainfall and their isotopic contents during the year, the convenience of using weighing approaches before applying OLS was not found.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/hydrology9020041/s1, Table S1: Catchment characteristics from the 9 rivers considered in this study, Table S2: Slopes from δ <sup>18</sup>O–δ <sup>2</sup>H regressions, Table S3: Intercepts from δ <sup>18</sup>O–δ <sup>2</sup>H regressions, Table S4: Standard deviations of the estimates from δ <sup>18</sup>O–δ <sup>2</sup>H regressions, Table S5: Coefficient of determinations from δ <sup>18</sup>O–δ <sup>2</sup>H regressions, Figure S1: Plots of standardised residuals for surface water affected by nonstationary processes, Figure S2: Autocorrelation functions of standardised residuals for surface water affected by serial correlations, Figure S3: Heat maps reporting slopes and intercepts values, Figure S4: δ <sup>18</sup>O–δ <sup>2</sup>H pairs from rain gauges, Figure S5: Dendrogram for slopes and intercepts series, Figure S6: δ <sup>18</sup>O–δ <sup>2</sup>H pairs from surface water, Figure S7: δ <sup>18</sup>O–δ <sup>2</sup>H pairs from groundwater.

**Author Contributions:** Conceptualization, F.C.; methodology, F.C. and A.T.; software, F.C.; validation, F.C.; formal analysis, F.C.; investigation, F.C.; data curation, F.C.; writing—original draft preparation, F.C.; writing—review and editing, F.C. and A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are freely available thorough Supplementary Materials.

**Acknowledgments:** Authors express their gratitude to two anonymous reviewers as their thoughtful and detailed comments allowed the early version of the manuscript to be greatly improved.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

