**Potential for Managed Aquifer Recharge to Enhance Fish Habitat in a Regulated River**

**Robert W. Van Kirk 1,\*, Bryce A. Contor 1, Christina N. Morrisett 2, Sarah E. Null <sup>2</sup> and Ashly S. Loibman <sup>3</sup>**


Received: 20 November 2019; Accepted: 18 February 2020; Published: 1 March 2020

**Abstract:** Managed aquifer recharge (MAR) is typically used to enhance the agricultural water supply but may also be promising to maintain summer streamflows and temperatures for cold-water fish. An existing aquifer model, water temperature data, and analysis of water administration were used to assess potential benefits of MAR to cold-water fisheries in Idaho's Snake River. This highly-regulated river supports irrigated agriculture worth US \$10 billion and recreational trout fisheries worth \$100 million. The assessment focused on the Henry's Fork Snake River, which receives groundwater from recharge incidental to irrigation and from MAR operations 8 km from the river, addressing (1) the quantity and timing of MAR-produced streamflow response, (2) the mechanism through which MAR increases streamflow, (3) whether groundwater inputs decrease the local stream temperature, and (4) the legal and administrative hurdles to using MAR for cold-water fisheries conservation in Idaho. The model estimated a long-term 4%–7% increase in summertime streamflow from annual MAR similar to that conducted in 2019. Water temperature observations confirmed that recharge increased streamflow via aquifer discharge rather than reduction in river losses to the aquifer. In addition, groundwater seeps created summer thermal refugia. Measured summer stream temperature at seeps was within the optimal temperature range for brown trout, averaging 14.4 ◦C, whereas ambient stream temperature exceeded 19 ◦C, the stress threshold for brown trout. Implementing MAR for fisheries conservation is challenged by administrative water rules and regulations. Well-developed and trusted water rights and water-transaction systems in Idaho and other western states enable MAR. However, in Idaho, conservation groups are unable to engage directly in water transactions, hampering MAR for fisheries protection.

**Keywords:** climate adaptation; stream temperature; streamflow; Henry's Fork; fisheries; Snake River; Idaho; water rights

#### **1. Introduction**

In the Western USA, important aquatic ecosystems and recreational fisheries often occur in river basins with large irrigated agricultural diversions, resulting in conflicts between water for irrigation and environmental streamflow needs [1–3]. Climate change exacerbates these conflicts, as precipitation regimes shift from snowfall to rainfall and evaporative demand increases, leading to flashier streamflow in winter and spring and reduced baseflow through summer and fall [4,5]. Climate warming and reduced baseflow work in tandem to warm stream temperature and are expected to reduce habitat for cold-water ecosystems [6]. Increasing streamflow, particularly during summertime baseflow conditions, cools stream temperatures by increasing the assimilative heat capacity of rivers [7]. Management options to increase streamflow include re-operating reservoirs [8,9], reducing diversions

through environmental water purchases [10], or conducting managed aquifer recharge (MAR) [11,12]. MAR is a promising strategy to enhance cold-water habitats while maintaining water resource benefits for people because excess water is intentionally recharged to raise aquifer levels, which increases baseflow, or may subsequently be pumped for irrigation. Groundwater additions to streams are particularly beneficial for cold-water fisheries because they create thermal refugia [13,14]. MAR is often recommended as a strategy to manage water for people and ecosystems with flashier runoff anticipated with climate change [12]. However, the potential of MAR to benefit cold-water ecosystems while maintaining irrigated agriculture requires (1) understanding the physical hydrology between the recharge site and the stream, (2) estimating temperature differences at groundwater seeps in the river and ambient temperatures, and (3) understanding administrative water rules to apply MAR to benefit cold-water habitat.

In streams that interact with local and regional aquifers, winter recharge enhances groundwater storage important for streamflow through the summer [14]. However, systems with shallow, unconfined aquifers are sensitive to climate variability [15] and may experience changes in the timing and magnitude of natural recharge [16], diminishing aquifer storage and groundwater-supported streamflow [16,17]. MAR can capture early rainfall or snowmelt and supplement late-summer return flows [18–20], by raising an aquifer's hydraulic head and creating groundwater seeps and shallow groundwater contributions that return to the stream in gaining reaches. Models have demonstrated that the lag time between MAR and return flow can delay the runoff peak [21–23], buffering a variable runoff regime [23] and potentially alleviating critical low-flow periods [24–26]. However, proportional contributions of MAR to streamflow depend, in part, on recharge site proximity [18,22,23].

Some studies have found that groundwater seeps and return flows mitigate the thermal effects of climate change on riverine habitats [27–29]. For example, measured water temperature at groundwater seeps have been 2–3 ◦C cooler than ambient river temperatures in the Pacific Northwest [21], and up to 4 ◦C cooler in Nevada's Walker River [12]. While shallow groundwater temperature is sensitive to long-term changes in air temperature, groundwater temperatures are less sensitive than surface water to changes in air temperature and are generally absent from heating by solar radiation [30–33]. Although studies note that MAR may increase summer baseflows [34], provide cool groundwater return flows to maintain cold-water fisheries during low-flow periods [26], and maintain aquatic ecosystems [18,35], field observations have yet to test these hypotheses. Furthermore, it is important to understand the extent and times that MAR can influence streamflow and stream temperature to maintain cold-water species in regulated rivers with climate change.

In the western USA, MAR must fit into the administrative rules of the Doctrine of Prior Appropriation, which allocates water for beneficial uses based on the seniority that water was first used. In most western states, the senior uses are mining and agriculture. Additionally, states must have well developed market and transfer mechanisms that provide administrative water for MAR within pre-existing allocation systems. However, western states that prioritize MAR, like Arizona, Colorado, California, and Idaho, each have different administration policies regarding which entities can implement MAR. Overall, implementation of MAR includes large-scale projects conducted by centralized public authorities, cities, and private companies in Arizona [36,37], smaller-scale projects implemented by landowners, local agencies, and counties in California [36,38], and a variety of MAR and recovery projects implemented by individual water right holders, local groundwater management districts, and cities in Colorado [37,39]. In Idaho, the state-run MAR program is primarily designed to increase discharge from the aquifer to the river for fulfillment of senior surface-water rights. Water transaction mechanisms in California and Colorado allow effective transfer of water to environmental uses, and conservation organizations can be direct participants in such transactions [40]. However, in Idaho, conservation groups cannot directly participate in transactions like water rental, which inhibits MAR for fisheries protection.

This study aims to understand the potential for MAR to benefit cold-water ecosystems while maintaining irrigated agriculture in Idaho's Henry's Fork Snake River. To do this, existing data and models were integrated and analyzed for a reconnaissance-level assessment. The following research questions were addressed: (1) What quantity and timing of streamflow response can MAR produce? (2) Does MAR increase streamflow by reducing channel loss to the aquifer or by increasing groundwater inflow to the stream? (3) Can groundwater inputs create local areas of decreased water temperature in the stream? (4) What legal and administrative hurdles exist for MAR to be used for cold-water fisheries conservation in Idaho? This research is instrumental to evaluate whether MAR is a water management strategy that has the potential to benefit managed fisheries. MAR is increasingly implemented throughout the western United States to meet human and environmental water objectives. There is a clear need to understand whether implementing MAR to benefit managed fisheries deserves additional effort and resources, and to identify current knowledge gaps when streamflow, stream temperature, and administrative water rules are considered for environmental MAR applications.

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

#### *2.1. Study Area*

The 93,000 km<sup>2</sup> upper Snake River basin in Idaho and Wyoming, USA, (Figure 1) is an ideal setting for assessing the potential of MAR to benefit cold-water ecosystems. The basin's water resources support an agricultural economy worth US \$10 billion [41], as well as many ecologically important stream systems and recreational trout fisheries, which contribute US \$100 million to local communities [42,43]. Mean annual surface water supply is 15,000 Mm<sup>3</sup> in the basin and 75% is withdrawn for irrigation. Irrigators also withdraw 1600 Mm3 of groundwater from the Eastern Snake Plain Aquifer (ESPA; Figure 1) and 500 Mm3 from tributary aquifers [41]. In a given year, over 10,000 km2 of irrigated land produce hay, wheat, barley, potatoes, and dairy products for global companies such as Anheuser-Busch, General Mills, and Clif Bar. The ESPA is a highly transmissive, unconfined, regional aquifer hosted in sediments interbedded within fractured Quaternary basalts [44]. Water generally flows through the ESPA from northeast to southwest and discharges to the Snake River near American Falls Reservoir and in a 100 km reach immediately upstream of King Hill (Figure 1). Water levels in and discharge from the ESPA have been declining for 60 years due to a combination of decreased recharge incidental to surface water irrigation and increased groundwater pumping [41,44,45]. Declining aquifer levels have caused costly legal disputes, increased reliance on reservoir storage to meet irrigation demand, increased groundwater pumping costs, and decreased streamflow for fisheries.

As part of a comprehensive plan to increase storage in and discharge from the ESPA, Idaho has implemented publicly funded MAR, with an annual objective of 330 Mm3 [41]. The primary management goal of the state's MAR program is to increase discharge from the aquifer to the river to fill senior surface water rights, rather than to store the water for future recovery via pumping. Increasing discharge over the long term requires increasing storage in the aquifer to maintain larger hydraulic gradients along connected river reaches. Higher storage volume, in turn, has ancillary benefits such as decreased pumping costs for groundwater users [41]. In addition, irrigation entities, cities, and private companies are using MAR on smaller scales to meet mitigation requirements of a 2015 legal settlement between senior surface water users and junior groundwater users. This settlement requires a specified reduction in groundwater pumping or mitigation with an equal amount of MAR. Concurrently, research describing benefits to aquatic systems from incidental and managed recharge in irrigated landscapes has motivated conservation organizations to consider MAR as a tool for maintaining and enhancing cold-water ecosystems in a changing climate [3,18,46].

An assessment was conducted in the Henry's Fork Snake River watershed (Figures 1 and 2), where the state of Idaho has recently invested US \$1.5 million to expand and improve a MAR site known as Egin Lakes (Figure 1). The Henry's Fork and its tributaries have an annual surface-water supply of 3200 Mm<sup>3</sup> and irrigators withdraw ~1500 Mm<sup>3</sup> to apply on ~1000 km<sup>2</sup> of agricultural land; very little groundwater is used for irrigation in the watershed [47]. When natural streamflow is insufficient to meet irrigation demand—usually early July through early September—streamflow is augmented by draft of Island Park Reservoir, near the river's headwaters (Figure 1). Nearly all irrigation water is delivered through unlined canals constructed in the late 1800s. Thus, these canals have provided a large amount of incidental recharge to local and regional aquifers via seepage for over a century [48,49]. Historically, irrigation water was applied via flooding or furrow irrigation, but most application was converted to sprinklers in the 1980s and 1990s [49,50]. The lower one-third of the Henry's Fork, shown as the "modeled reach" in Figure 1, is hydraulically connected to local and regional aquifers. Previous research has shown that this reach gains water seasonally in response to locally increased water tables during irrigation season, but loses water to the regional ESPA during the winter [48,49,51,52]. The conversion to more efficient sprinkler application has reduced both total diversion and groundwater return flows to the river by around 250 Mm<sup>3</sup> per year since the late 1970s [49].

**Figure 1.** Upper Snake River basin, USA, showing the ESPA, modeled river reach of Henry's Fork (in green), and a nearby MAR site. The red polygon delineates the Henry's Fork watershed, and the yellow rectangle shows the area enlarged in Figure 2. Black arrows indicate primary groundwater flow paths on the ESPA. Data credit: Idaho Department of Water Resources.

The "field study reach" is the 12 river-km of the Henry's Fork immediately downstream of U.S. Geological Survey streamflow gage 13050500 at St. Anthony (Figure 2). This gage is the streamflow management point in the lower watershed, triggering additional releases from Island Park Reservoir when streamflow drops below a specified target at this gage [53]. However, four diversions in the field study reach downstream of the St. Anthony gage substantially reduce streamflow during the summer (Figure A1). The study reach supports an increasingly popular and economically valuable recreational sport fishery for wild brown trout (*Salmo trutta*) [54,55], which has an optimal summer temperature range of 12 ◦C to 19 ◦C. Habitat suitability for brown trout decreases as temperature increases above 19 ◦C to the lethal limit of 27 ◦C [56]. Over the summers of 2016, 2017, and 2018, daily mean water temperatures during July and August ranged from 16 ◦C to 20 ◦C at a water-quality monitoring station

at the top of the field study reach and from 17 ◦C to 22 ◦C at a water-quality monitoring station at the bottom of the field study (Figures 2 and A2). Maximum instantaneous water temperature recorded at the lower station over this time period was 27.3 ◦C, and daily maxima frequently exceeded 22 ◦C. Due to high water temperatures, brown trout move either to local areas of groundwater input or out of the reach altogether during the summer [54].

**Figure 2.** Map of field study reach, showing locations of U.S. Geological Survey (USGS) streamflow gage at St. Anthony, temperature loggers deployed in 2010, water-quality monitoring locations, and stretch where springs were documented in 2019. Data credit: ESRI.

Summertime streamflow in the study reach could be increased by increasing draft of Island Park Reservoir 100 km upstream, but larger reservoir releases have numerous negative effects on other popular and economically important fisheries in the upper half of the watershed [3]. These include transport of suspended material out of the reservoir and resulting high turbidity during the peak fishing season, increased water temperatures downstream of the reservoir when it is drafted faster than thermal stratification can occur, and decreased trout survival during winter, when low outflow is required to refill the reservoir [57–59]. These effects do not propagate downstream to the study reach in the lower watershed. Thus, this study seeks to assess whether MAR has the potential to improve fisheries in the lower watershed without degrading those in the upper watershed. In particular, withdrawal of water for MAR at carefully identified times could increase groundwater inputs to the lower river during the summer, thereby increasing local trout habitat and water supply available for diversion there. In turn, increased summertime water supply in the lower river could limit reservoir draft, thereby simultaneously benefiting fisheries in the upper watershed.

#### *2.2. Streamflow Response*

We used an existing regional groundwater model to estimate response of streamflow in the Henry's Fork to MAR at the Egin Lakes site, located 8 km from the Henry's Fork (site location is shown in Figure 3, which also depicts the results). Modeling was done with the Idaho Department of Water Resources' Enhanced Snake Plain Aquifer Model Version 2.1 (ESPAM2.1) [52], a regional finite-difference flow model implemented in MODFLOW and configured with a single aquifer layer, monthly temporal

resolution, and roughly 11,000 1.6-km grid cells. The model supports both steady-state and transient simulations. Although storage coefficients are typical of unconfined conditions, the transient rendition of the model uses time-constant aquifer transmissivity, making model results additive and scalable. The model was calibrated to 1980–2008 hydrologic conditions, using the first five of these as a burn-in period [52]. Calibration used over 43,000 aquifer water levels, 2000 river gain and loss estimates, and 2000 spring discharge measurements and was performed using PEST version 12.0, a nonlinear parameter estimation program [60]. The model was built specifically to estimate effects of aquifer pumping and recharge on river reaches and springs in hydraulic connection with the aquifer, so calibration optimized groundwater-surface water exchanges rather than hydraulic heads. The model does not simulate solute transport nor thermal changes in the aquifer or its discharge. In the model, the hydraulically connected section of the Henry's Fork is treated as a single, 75-km reach, referred to as the "modeled reach" (Figure 1), whereas our field study reach is only 12 km in the center of the ESPAM2.1 modeled reach (Figures 1 and 2). Monthly model calibration residuals for stream gains and losses in the modeled reach of the Henry's Fork were generally on the order of ± 25%, but monthly residuals as large in magnitude as −100% were observed early in the irrigation season. Over the period 1985–2008, the model underestimated cumulative river gain by around 10%. Thus, the model is suitable for our reconnaissance-level assessment when applied over long time periods.

**Figure 3.** Steady-state discharge response in the modeled river reach to recharge conducted in a given model cell, as a fraction of total recharge volume. For example, a response fraction of 0.45 (white cells) indicates that 45% of the volume of water recharged in that cell will eventually contribute to streamflow in the modeled reach. The yellow rectangle indicates the field study reach.

The model was used in two ways. First, a steady-state simulation was used to calculate the fraction of total recharge in a given model cell that affects streamflow in the modeled reach of the Henry's Fork. Recharge was simulated in the model cells containing the Egin Lakes MAR site, as well as in other model cells in the vicinity of the lower Henry's Fork to assess whether developing MAR sites in other locations could increase streamflow response to MAR. Second, a 30-year transient simulation was used to estimate streamflow response in the modeled reach to a scenario in which 7.3 Mm<sup>3</sup> of water was withdrawn from the river and recharged at the Egin Lakes site in each of March, April, and October and 25.6 Mm<sup>3</sup> of water was withdrawn from the river and recharged at Egin Lakes during November. This annual scenario was similar to operation of the Egin Lakes site during 2019 and was repeated in each of the 30 years of simulation. Model output was used to estimate net change to streamflow in the study reach by allocating total streamflow response for the modeled Henry's Fork reach proportionally to the field study reach and including the effect of diversion for MAR upstream of the study reach. The median 2000–2019 hydrograph for streamflow at the bottom of the study reach was used as a baseline condition, although the effect of modeled MAR was also assessed relative to streamflow in 2016, the driest year in the basin in over 40 years.

Recharge proximal to the lower Henry's Fork increases hydraulic gradients between the aquifer and the river, but if these gradients were initially negative (i.e., water flows from the river to the aquifer), a positive streamflow response from recharge would occur through decreased river losses rather than through increased river gains. Although the resulting increase in streamflow would be equivalent between the two mechanisms, the first mechanism would not provide the benefit of decreased water temperature during the summer. Thus, summertime water temperature was measured upstream and downstream of a reach known to be hydraulically connected with the underlying aquifer. These measurements were conducted in 2010, a decade after conversion from flood to sprinkler irrigation but six years prior to initiation of MAR at the expanded Egin Lakes site. Canal seepage, which has been roughly constant since 2000, was the only source of local groundwater recharge in 2010. Water temperature loggers were deployed from 1 June 2010 to 31 August 2010 at two locations in the upper half of our field study reach (Figure 2) and secured underneath overhanging riparian vegetation at ~40 cm water depth. The upstream logger was located immediately downstream of a reach through which the river flows over basalt bedrock and has little interaction with shallow groundwater. The other logger was placed 5 km downstream, in a reach where the river is well connected with shallow groundwater. Mean daily upstream temperature was subtracted from downstream temperature to create a time series of temperature differences. After accounting for serial autocorrelation with lag-3 autoregressive terms, two statistical models were fit to the time series—one with hypothesized zero mean and another with a non-zero fitted mean. Statistical significance of the fitted mean was assessed with the likelihood ratio test at a 0.05 level of significance (Appendix B).

#### *2.3. Local E*ff*ects of Groundwater Inflow on Temperature*

Whereas the 2010 temperature observations were made to assess the nature of summertime streamflow response to recharge solely from canal seepage, a separate field study in 2019 documented the locations and temperature of specific groundwater springs to investigate the potential for MAR to provide cool groundwater return flow to the river. In July 2019, locations of groundwater springs contributing water to the river channel were documented by walking a 1-km length of the right bank of the river, in the lower half of the 5 km reach studied in 2010 (Figure 2). A steep bluff approximately 5–10 m in vertical relief forms the boundary of the active floodplain on the right side of the river. Springs emerged from the face of the bluff and along its base, often between 1 and 50 m from the channel bank, and most spring outputs flowed into a secondary river channel. Each spring site was classified as (1) a single discharge point, where water originating from the bluff face created a separate channel that actively flowed into the river, or as (2) a "wall seep", where water emerged from continuously saturated sediments along the bluff face and contributed unchannelized flow to the river. At each spring site, a FLIR T450sc thermal infrared camera (FLIR Systems Inc., Wilsonville, Oregon, USA) was used to document differences between the surface temperatures of incoming groundwater springs and the river. The FLIR T450sc camera senses radiant stream surface temperatures in the 7.5 to 13 μm range, with an accuracy of ± 1 ◦C or 1% of the range of the reading [61]. To complement the imagery, instantaneous temperatures were measured with a handheld thermometer at three lateral locations: spring emergence, 0.6 m and 6 m into the river channel from where the spring entered the river. At wall seeps, lateral temperatures were measured at the upstream and downstream extent of each seep. In total, three temperature measurements were recorded at each of 20 spring sites. Temperature

differences across the three lateral locations were analyzed with mixed-effects analysis of variance and Tukey's post-hoc test, treating spring site as a random effect. These tests were conducted at a 0.05 level of significance (Appendix B). The temperature analysis was not accompanied by assessments of whether physical habitat at these locations was otherwise suitable for or used by trout.

#### *2.4. Water Administration*

Potential streamflow and temperature benefits of MAR will not result in real changes in the river without sufficient availability of water for MAR at appropriate times. Thus, our assessment included analysis of physical and administrative availability of water for MAR in the upper Snake River basin within Idaho's prior appropriation system of water rights. This assessment relied on a formal review of Idaho's MAR program conducted for the Idaho Water Resource Board [62], to which two of the co-authors of this paper (RVK and CNM) contributed substantially. In addition, the state's water rights database, water-rights accounting manual [63], and water exchange procedures [64,65] were reviewed to identify opportunities for and limitations to conducting MAR for fisheries conservation purposes.

#### **3. Results**

#### *3.1. Streamflow Response*

The steady-state simulation using ESPAM2.1 predicted that 37% of the water volume delivered to the Egin Lakes MAR site will increase streamflow in the modeled reach of the Henry's Fork over the long term, and the balance will benefit other river reaches in the basin (Figure 3). If recharge were conducted closer to the river than the existing MAR site, the model predicted that >90% of recharge is realized as increased streamflow in the modeled reach. The modeled response fraction depended strongly on recharge location and decreased fairly rapidly with increasing distance between the river and recharge location (Figure 3). For example, 90% of water recharged in the red cells contributed to streamflow in the modeled reach, whereas less than 40% of the water recharged in the green cells contributed to streamflow in the modeled reach.

Transient simulation with ESPAM2.1 predicted that streamflow response to spring and fall recharge at Egin Lakes is relatively uniformly distributed over the year, with little month-to-month variability (Figure 4). Initial streamflow response to the spring-fall MAR scenario increased roughly linearly over time to reach 50% of its long-term value 6.5 years after first implementation of the annual MAR regime (Figure 4). Streamflow response increased more slowly after that, reaching ~90% of its long-term value 25 years after initial implementation of the annual MAR regime. Including the effects of water withdrawal from the river for delivery to the MAR site, the annual MAR scenario resulted in a 20%–25% decrease in streamflow during November and a 5%–10% decrease in streamflow during each of October, March, and April, relative to the current median hydrograph (Figure 5). Despite these decreases in streamflow due to withdrawal for MAR, median spring and fall streamflow still remained much higher than summertime lows, even after only five years of MAR. After 20 years of implementation of this annual MAR regime, median streamflow increased by 4%–7% during July and August of the median year (Figure 5), although increases during late June and early July were on the order of 10%–40% relative to streamflow in the dry year of 2016 (Figure A3).

Mean daily water temperature from 1 June 2010 to 31 August 2010 was 0.6 ◦C cooler at the downstream location influenced by groundwater inputs, and this difference was statistically significant (χ<sup>2</sup> = 5.3, df = 1, *p* = 0.02). This indicates that during summer, streamflow response to seasonal aquifer recharge results from inflow of groundwater to the river not from reduced loss of water from the river to the aquifer. Since this result was observed when canal seepage was the only source of aquifer recharge in the vicinity of the field study reach, additional recharge from MAR will further increase flow of groundwater into the river in the field study reach.

**Figure 4.** Recharge and discharge for annual MAR scenario of 7.3 Mm3 per month in each of March, April, and October, and 25.6 Mm3 per month in November, repeated every year for 30 years.

**Figure 5.** Net change in streamflow for annual scenario in which diversion for MAR from the study reach is 7.3 Mm<sup>3</sup> per month in each of March, April, and October, and 25.6 Mm3 per month in November. Top panel shows median water-year hydrograph prior to and 20 years after initiation of annual MAR regime. Bottom panel shows percent change in streamflow 5, 10, and 20 years after initiation of annual MAR regime, respectively.

#### *3.2. Local E*ff*ects of Groundwater Inflow on Temperature*

In late July 2019, thermal imagery identified areas of cool water in the main river and its side channels near the points of spring inflow (Figure 6). Mean instantaneous water temperature at the 20 spring sites differed significantly across the three lateral locations: spring, 0.6 m from the streambank, and 6 m from the bank (F = 29.7, df1 = 2, df2 = 38, p < 0.001). All pairwise differences among the locations were significant (Tukey's Honest Significant Difference, adjusted p < 0.001). Mean water temperatures at the lateral locations, respectively, were 14.4 ◦C, 16.0 ◦C, and 18.3 ◦C (Figure 7). Ambient water temperatures during the time of the field observations ranged between daily minima of 18 ◦C and daily maxima of 23 ◦C at the top of the field study reach and between 18 ◦C and 25 ◦C at the downstream boundary of the study reach (Figure A2).

**Figure 6.** A side-by-side comparison of a visual image (**left**) and thermal infrared image (**right**) of the area where outflow from a groundwater spring located at 43◦57'07.6" N 111◦43'26.7" W entered a side channel of the river (flowing right to left). The photo was taken from a point 1.5 m from the margin of the side channel, looking toward the spring confluence. The spring emerged from the ground ~30 m from the confluence point. Temperature is indicated in ◦C.

**Figure 7.** Water temperature at three locations measured at each of 20 distinct springs.

#### *3.3. Water Administration*

There are 84 decreed, permitted, or pending water rights for MAR in the upper Snake River basin, with a combined diversion rate of up to 725 m3/sec. However, every senior MAR right is very small, with combined diversion of only 0.59 m3/sec. The remaining 724.41 m3/sec are distributed among large water rights with priority dates of 1980 or later, in a prior-appropriation system where the majority of the irrigation rights have priority dates preceding 1910 and most reservoir storage priority dates precede 1940. Diversion for MAR allowed by the junior rights is available on an annual basis only in the winter and then only downstream of American Falls Reservoir (Figure 1). Considering only water rights specific to MAR, water is available for MAR at the Egin Lakes site in about half of all water years, usually between mid-May and early July. During irrigation season (1 April to 31 October), only water delivered to a designated, off-canal MAR site can be accounted as MAR, but during the winter, canal seepage also accounts as MAR, since canals would not customarily be delivering irrigation water then. Temporary transfers of senior water rights from irrigation or other uses to MAR occur through the state of Idaho's water supply bank, an administrative exchange bank rather than a physical storage bank such as those in Arizona and California.

A locally administered rental pool allows storage water held in Palisades Reservoir (Figure 1) to be rented for delivery to MAR sites anywhere in the basin, through administrative exchange or physical delivery, and these types of exchanges were used to provide the majority of water for MAR in the modeled scenario. Storage water rented in an administrative year (1 November to 31 October) must be delivered by December 1 and cannot be carried over any further into the subsequent year. Since the 2015 ESPA groundwater-surface water settlement was completed, new administrative rules have been enacted specifically to facilitate efficient but equitable water transactions for MAR. However, entities that do not hold water rights, including most conservation groups, cannot participate directly in water supply bank or rental pool transactions. Furthermore, the Idaho Water Resource Board is the only entity that can hold surface water rights for instream flow, regardless of whether those rights are permanent or temporary. Thus, even if conservation groups could participate in water transactions, there is no precedent for them to hold MAR rights specifically for environmental purposes.

#### *3.4. Limitations*

The importance of understanding the temporal response to water management actions is critical to address specific objectives. The one-month temporal resolution of ESPAM2.1 limits its ability to predict streamflow response during shorter time intervals that may be critical to trout survival. Furthermore, the model cannot distinguish between management actions that contribute groundwater and those that reduce streamflow losses to groundwater. Although our limited temperature observations suggest that summertime streamflow response to recharge near the field study reach is realized as groundwater inflow, the existing ESPAM2.1 model cannot predict whether MAR and other recharge strategies will change water temperature. The 1.6 km spatial resolution of the model also limits predictive use, especially in assessing response to MAR at hypothetical sites closer to the river, where response changes rapidly with distance away from the river. However, the greatest spatial-resolution limitation in ESPAM2.1 is the delineation of the river reach itself. The model cannot partition water across the 75 km lower Henry's Fork model reach to different locations, requiring simple proportional allocation to downscale model results, as was done here. Models with finer temporal and spatial resolution can be constructed [49], but calibrating higher-resolution models requires hydrologic observations made at the same scale, which are currently not available. A larger challenge to modeling groundwater flow in the Henry's Fork watershed is that the water budget cannot be closed using surface-water observations alone, whereas that for the larger ESPA can. This requires a priori assumptions about groundwater flux to calibrate parameters specifying boundary conditions. Groundwater flux is the most important model output in this case, so model output would essentially be pre-determined by assumptions required to calibrate boundary parameters.

An alternative approach to constructing finer-scale models is to estimate local groundwater fluxes using fine-scale piezometer data and temperature mixing models, and use those fine-scale models to downscale the coarse aquifer-river interactions predicted by the regional model. Models can be verified by conducting measurements of streamflow gain and loss across short stream reaches. Continuous, high-resolution temperature data are inexpensive to collect and would contribute not only to temperature mixing models but also to identification of thermal refugia across the whole field study

reach. Habitat surveys, observations of fish movement and habitat use, and water-quality analysis at and near areas of cooler water temperatures could then be used to determine whether reducing water temperature alone is sufficient to address factors limiting trout use of the study reach during summer.

#### **4. Discussion**

Overall, MAR can improve summer streamflow and stream temperatures for fish in localized areas of the basin. Thirty-seven percent of modeled Egin Lakes MAR returned to the study reach. Streamflow increased most (as a percentage change of baseflow) in the initial years following recharge, with 25–30 years needed to achieve steady-state response to an annually repeated MAR regime. Groundwater seeps confirmed that recharge was contributing to the river rather than merely reducing losses from the river to the aquifer. On average, stream temperature cooled 0.6 ◦C after traveling 5 km downstream during summer. July 2019 groundwater seep temperatures averaged 14.4 ◦C and were 16.0 ◦C about a half-meter from where groundwater seeps entered the river. These temperatures are in the suitable range for brown trout [56]. Average mid-summer ambient temperature was 20.1 ◦C for 2016 through 2018, which exceeds the optimal thermal brown trout threshold of 19 ◦C [61]. These findings suggest that MAR provides thermal refugia for managed fish species during summer. A review of administrative MAR rules for fisheries conservation provided more equivocal results. Some water is available for MAR in spring and fall of the wettest half of years, but substantial water for recharge is not consistently available, and canal capacity to the Egin Lakes MAR site limits recharge when water is available. Conservation organizations cannot participate directly in water transactions, but must partner with irrigators or water rights holders on MAR projects.

#### *4.1. Physical and Administrative MAR Implications for Idaho and Henry's Fork*

To maintain cool summer water temperatures at the reach scale for cold-water species, MAR volumes will likely need to offset declines in recharge that have occurred from improved irrigation efficiency, as sprinklers and center pivots have replace flood irrigation over the past few decades. Improved irrigation efficiency typically does not increase streamflow as more land is put into production or junior water rights come into priority [66,67]. Around 250 Mm<sup>3</sup> is needed annually to offset lost incidental recharge, which exceeds the physical capacity of the existing MAR site. This volume could be attained if recharge occurred year-round, so that winter canal seepage also contributes. Modeling showed that spring and fall recharge contributes to summer streamflow. Since only 37% of the total volume recharged at the Egin Lakes MAR site returns to Henry's Fork, the benefits to summer streamflow must be weighed against the negative effects of withdrawing larger amounts of water from the river at other times of year. These include negatively impacting aquatic habitat availability, species life history expressions (e.g., spawn timing), or fluvial geomorphic processes (e.g., floodplain maintenance). Developing MAR sites closer to the river could increase streamflow in the study reach per unit of water withdrawn.

Fine-scale field observations conducted in 2019 showed that at seep locations, groundwater inputs can cool ambient stream temperature by over 2 ◦C during summer, a difference that is biologically significant for trout. Even if cool groundwater inputs are not widespread across river reaches or the contributions are not large compared to river flow, springs may create local thermal refugia for fish, allowing greater survival throughout the summer than would otherwise occur given the same physical habitat availability and streamflow [68]. Further identification of groundwater inflows, their hydrogeologic properties, and their use by fish could inform water and fisheries management actions to enhance groundwater springs for fish populations. Understanding the connectivity of thermal refugia will also help managers understand where fish can become trapped or bottlenecks occur for movement [13,68].

However, it is also important to understand the quality of these groundwater return flows on these groundwater-dependent ecosystems. When conducted on working agricultural lands, MAR risks mobilizing nutrients and increasing chemical constituent loading to streams and riparian soils [20,69–71]. MAR can also facilitate groundwater contamination via crop-mediated aquifer contamination of fertilizers

and pesticides [72]. Whereas aquifer recharge may introduce water quality concerns elsewhere, return flow from the ESPA is of high quality that is suitable for instream habitat uses [73]. This research is supported by monitoring conducted at Egin Lakes in 2019 by the Idaho Department of Environmental Quality that showed no increase in nitrogen or fecal coliform as a result of MAR (Aaron Dalling, Fremont-Madison Irrigation District, 2019, presentation to Henry's Fork Watershed Council, October 22). Thus, MAR operations in the Henry's Fork avoid such groundwater contamination and pollutant leaching by conducting recharge via canal seepage and infiltration at Egin Lakes.

While improved modeling and detailed field work can provide technical understanding of MAR benefits to summer habitat for trout, administrative and logistical hurdles must be overcome for implementation [74,75]. The most basic of these is the junior priority dates of water rights for MAR. In the larger context of prior appropriation and development of Idaho's water resources, MAR is a relatively new administratively-recognized beneficial use of water in a system in which water rights for agriculture and mining date back to the mid-19th century. The Idaho Water Resources Board implemented aggressive groundwater and conjunctive management policies and procedures in the 1980s and 1990s, including obtaining large-volume MAR rights with 1980 and 1998 priority dates. Idaho also has flexible and well-established water transaction mechanisms to facilitate transfer of senior water rights to MAR. In the Henry's Fork watershed, the state's MAR water rights are in priority in only half of all water years and then usually during irrigation season, when the canal system is already near capacity delivering irrigation water. Although new MAR infrastructure has been built throughout the basin since 2009, the majority of conveyance to MAR facilities occurs through the existing irrigation canal system. Additional canal capacity could alleviate this limitation during wet years but would go unused in the other half of years.

Because of summer canal capacity limitations, costs of new infrastructure, and the junior priority of MAR rights, storage water rental and other exchange mechanisms offer the greatest potential to increase MAR volumes. For example, in 2018 and 2019, reservoir storage water rented by groundwater users to meet mitigation requirements of their legal agreement with surface water users was not needed for direct delivery to the surface water users because of good water supply in those years. Instead, the rented water was assigned to the state for MAR, allowing recharge of an additional 84 Mm3 in 2019 over what would have been available using the state's junior MAR rights alone (Wesley Hipke, Idaho Department of Water Resources, 2019, data distributed to stakeholders via email, December 6). The need to manage and administer groundwater and surface water conjunctively to meet the legal requirements of the ESPA settlement creates a pseudo-market to fund such exchanges. In 2017, 2018, and 2019, some of the water recharged at the Egin Lakes site was made available for MAR through water exchanges, including those described above, and delivered in spring and fall, outside of peak irrigation season.

Ideally, conservation groups could facilitate incentive-based irrigation reduction, rent the saved storage water for MAR, and keep that water in reservoirs throughout the summer, thus reducing negative effects of reservoir drawdown. The rented storage water could be diverted for MAR in the off-season, when natural streamflow is sufficient without reservoir releases. The separation of physical and administrative water that makes this possible is routine in Idaho under current administrative procedures. Canal capacity to deliver water for MAR is greater in the winter, when canals are not also delivering irrigation water, although the basin's sub-freezing temperatures present challenges in managing winter canal delivery and frozen soil can impede infiltration [76]. Canal seepage during the winter adds recharge that would not otherwise occur and is not simply an administrative replacement for historical seepage incidental to irrigation that occurs during the summer. Rental pool water is available in larger quantities and at lower prices during wet years, potentially allowing for more MAR during wet years, which would then contribute to streamflow in subsequent dry years. However, changes to water rental rules are needed to allow storage water rented in one administrative year to be diverted for MAR several months into the next administrative year [64]. In addition, conservation organizations cannot implement MAR projects, but instead must partner with irrigation entities or individual water rights holders. To fully capitalize on the high value anglers place on upper Snake

River fisheries, new administrative and transfer mechanisms are required to allow conservation organizations to participate directly in water transactions.

#### *4.2. Physical and Administrative MAR Implications for other States in the Western USA*

Most states in the western USA possess some of the physical and administrative features required for MAR to benefit cold-water ecosystems, but few have all of the requisite ingredients. Arizona's progressive and flexible administrative systems have been highly successful in facilitating recharge of Colorado River water, but the sole purpose of this MAR is to store the water for later recovery, not to enhance streamflow [36]. On the other hand, restrictive administrative rules in Colorado limit use of MAR in headwater alluvial aquifers [37], where snowmelt could be captured and recharged with the intent of enhancing streamflow later in the summer, and where Colorado's progressive water markets would allow conservation organizations to obtain this water for environmental uses [40]. California is now conducting "flood-MAR" in depleted aquifers using stormwater [38], but this water is junior to other rights and is available for MAR only when existing water rights and required environmental uses are fulfilled. Some regions of California are considering creating an environmental water account using MAR, although this idea is currently untested.

#### *4.3. MAR as Climate Adaptation Strategy for Fisheries Conservation*

Cold-water fish habitat is anticipated to decline substantially with climate change. In fact, brown trout are expected to lose 48% of their habitat in the interior western US by the 2080s [68]. These changes are driven by warmer stream temperature and increasing winter floods, and could have major repercussions for a local economy reliant on a US \$100 million recreational fishery. Additionally, an increase in extreme climate events—i.e., a higher frequency of wet and dry years, with fewer 'normal' years is also expected [77]—which will alter instream conditions for biota [78]. Together, these changes suggest that climate adaptation strategies that provide mechanisms to reduce winter flooding, increase summer baseflow, cool summer stream temperature, and enhance thermal refugia are warranted. MAR for fisheries conservation is one such strategy. MAR is a promising climate-adaptive water management strategy because winter flows can be recharged to underlying aquifers to maintain baseflow and cold-water fish habitat throughout the year. However, junior water right holders will have considerable uncertainty from the increased inter-annual variability inherent with climate change [79]. Although this does not inherently reduce the utility of MAR for fisheries conservation in a warming climate, it does suggest that relaxing current administrative rules for greater flexibility to carry over reservoir rental water between years would improve the utility of MAR for fisheries conservation in a changing climate. Since many important recreational trout fisheries in the Western U.S. are located downstream of reservoirs, renting reservoir storage not used for irrigation in a given season and using it for MAR during the subsequent off-season is an innovative conservation mechanism that could have wide applicability.

#### **5. Conclusions**

MAR during spring and fall is expected to increase streamflow by around 5% during mid-summer, but only after 20 years of consistent MAR. Developing MAR sites closer to the river than the existing site would provide greater streamflow benefit per unit of water withdrawn and on shorter time frames. However, even relatively small increases in streamflow could have disproportionately greater benefits to trout, as streamflow response in our study area occurs in the form of increased groundwater inputs rather than decreased losses to the aquifer. July water temperatures were locally 2 ◦C cooler where groundwater flowed into the river. Because MAR is a new and junior water use in a priority system with irrigation rights dating back to the 19th century, MAR water rights are generally in priority only during late spring and early summer of years of above-average supply. However, administrative exchange that allows reservoir storage to be used for MAR can make water available during spring and fall, when MAR infrastructure capacity is greatest. Changes to current administrative rules could increase the effectiveness of such exchange mechanisms in providing water for MAR.

**Author Contributions:** Conceptualization: R.W.V.K.; methodology: R.W.V.K., B.A.C., C.N.M., and A.S.L.; software: B.A.C.; validation: R.W.V.K., B.A.C., and C.N.M.; formal analysis: R.W.V.K.; investigation: A.S.L., R.W.V.K., and C.N.M.; resources: R.W.V.K. and S.E.N.; data curation: R.W.V.K., B.A.C., and C.N.M.; writing—original draft preparation: R.W.V.K., B.A.C., C.N.M., and S.E.N.; writing—review and editing: B.A.C., C.N.M., S.E.N., R.W.V.K., and A.S.L.; visualization: C.N.M., B.A.C., and R.W.V.K.; supervision: R.W.V.K., C.N.M., and S.E.N.; project administration: R.W.V.K. and C.N.M.; funding acquisition: R.W.V.K. and C.N.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** Collection of the 2010 water temperature data was funded by U.S. Department of Agriculture grant 2008-51130-19555 to R.V.K. C.N.M. received support from National Science Foundation grant no. 1633756. A.S.L. was supported during summer 2019 by an internship from Colgate University. Additional funding was provided by the Federal Highway Administration and Fremont County Idaho.

**Acknowledgments:** Matt Hively assisted with field data collection. Input from the academic editor and five anonymous reviewers greatly improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Appendix A. Streamflow and Temperature Graphs**

**Figure A1.** Streamflow at the St. Anthony gage (top of field study reach) and downstream of all diversions (roughly at the upstream extent of the springs identified in Figure 2).

**Figure A2.** July–August water temperature at the St. Anthony water quality station (top of field study reach) and at the Parker-Salem water quality station (bottom of field study reach). Shaded polygon is optimal thermal range for brown trout. Horizontal line is lethal temperature limit for brown trout. The dates of the 2019 temperature observations are shown for reference.

**Figure A3.** Net change in streamflow in a dry year for annual scenario in which diversion for MAR from the study reach is 7.3 Mm<sup>3</sup> per month in each of March, April, and October, and 25.6 Mm3 per month in November. Top panel shows the 2016 water-year hydrograph and hypothetical effect of MAR implemented 20 years prior. Bottom panel shows percent change in observed 2016 streamflow 5, 10, and 20 years after initiation of annual MAR regime, respectively.

#### **Appendix B. Statistical Methods**

#### *Appendix B.1. Time Series Analysis*

Statistical hypothesis tests require independent observations for correct distribution of test statistics [80]. In time series such as daily water temperatures, observations are not independent of one another because of correlation between a given observation and the observations that precede it in the time series. This is referred to as serial or temporal autocorrelation. Its effect must be removed to obtain independence of observations before conducting hypothesis tests on time series. Autoregressive (AR) models are used to accomplish this [81]. The simplest AR model is the first-order model:

$$y\_t = \mu + \phi\_1 (y\_{t-1} - \mu) + \varepsilon,\tag{A1}$$

where *yt* is the observation at time *t*, μ is the mean of the time series, φ<sup>1</sup> is the first-order autoregressive coefficient, and ε is random, independent, normally distributed error. The autocorrelation term φ1(*yt*−<sup>1</sup> − μ) removes the dependence of *yt* on *yt*−1, allowing hypothesis tests to be conducted on the mean μ. In our case, serial autocorrelation was high enough that observations were correlated with those one, two, and three time steps prior. The resulting third-order model has the same form as Equation (A1), but with three autocorrelation terms. Because observation *yt* in our time series was the difference in temperature between the two locations, our null hypothesis was μ = 0. If rejected, we infer μ - 0. We used a significance level of 0.05, which is the default standard in statistical hypothesis testing. It represents the probability of having made an error in rejecting the null hypothesis, referred to as Type 1 error [80].

#### *Appendix B.2. Tukey's Post-Hoc Test*

Analysis of variance tests the single null hypothesis that all group means are equal. If the null hypothesis is rejected at a given significance level, the alternative hypothesis is simply that at least one group mean differs from at least one other. Additional tests must be done to assess which group mean(s) differ from which others. The probability of committing Type 1 error is compounded each time an additional test is performed. In our case, testing all possible differences between three group means requires three tests. If each is performed at a 0.05 significance level, the probability of committing at least one Type 1 error across the three tests is roughly 0.14. Tukey's post-hoc test is a method for conducting the tests for differences across all pairs of groups while maintaining the desired level of significance across the whole family of tests [80].

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Suitability Mapping for Managed Aquifer Recharge: Development of Web-Tools**

**Jana Sallwey 1,\*, Robert Schlick 1, José Pablo Bonilla Valverde 2, Ralf Junghanns 1, Felipe Vásquez López <sup>3</sup> and Catalin Stefan <sup>1</sup>**


Received: 12 August 2019; Accepted: 21 October 2019; Published: 28 October 2019

**Abstract:** Suitability maps for managed aquifer recharge (MAR) are increasingly used and hold the potential to be integrated into sustainable groundwater management plans. However, the quality of the maps strongly depends on the input data quality as well as the expertise of the decision-maker. The maps are commonly derived through GIS-based multi-criteria decision analysis (GIS-MCDA). To date, there is no common understanding of how suitability mapping should be conducted, as there is considerable variability concerning used GIS data and MCDA methodology. This study presents two web-tools that were conceptualized based on a review of GIS-MCDA studies in the context of MAR suitability mapping. The data retrieved from the review was compiled into a web-based query tool making the MAR- and MCDA-relevant information easily accessible. Based on the most commonly used MCDA practices in the assessed studies, we conceptualized and implemented a second web tool that comprises a simplified web GIS as well as supporting tools for weight assignment and standardization of the criteria. Both web tools will enable decision-makers to engage in MCDA for MAR mapping in a more structured and informed way. As the tools are open-source and web-based, they can facilitate the collaboration between multiple stakeholders and the easy sharing of results.

**Keywords:** managed aquifer recharge; web GIS; web tools; multi-criteria decision analysis; suitability mapping

#### **1. Introduction**

The application of managed aquifer recharge (MAR) is continuing to grow worldwide as a measure for sustainable groundwater management [1,2]. Before MAR schemes can be developed, comprehensive planning is required to ensure their long-term sustainability. While guidelines on the planning of MAR schemes exist [3–6], they mostly focus on their design and operation and put less focus on site selection. The selection of sites suitable for MAR is a critical step in the planning phase of a MAR project, as the location influences the recharge technique as well as the operational and maintenance parameters, such as the infiltration quantity and the recovery efficiency [7–10]. Site selection for MAR application is mostly conducted through field investigations. Suitability maps that show the potential of a foreseen area for the application of a certain MAR type can be generated as a preliminary step to field investigations. These maps are increasingly being used [11] and may fill a void in missing strategic MAR site planning. Their advantages for water management plans lie within the spatial display through maps [12], the quickness and simplicity of the analysis [13], the possibility to include projections of climate scenarios, population growth or land-use changes [14] as well as the assessment of different MAR techniques and their location [15].

While these maps are increasingly being used, there are no common guidelines on how the suitability mapping process should be conducted. The maps are generally generated by combining geoinformation of the surface and the subsurface with socio-economic criteria. This can be achieved by integrating multi-criteria decision analysis (MCDA) for solving spatial problems with GIS software [16]. A set of geospatial data must be chosen based on the study's objectives. The different GIS criteria are then weighted based on their importance for the study and combined into a suitability map. MCDA comprises a variety of methods for criteria weighting and combining [17]. A study showed that the GIS maps used and the methodologies applied for MCDA in the context of MAR, show a great variety [11]. The choice of GIS data and how important each dataset is seen for each study is dependent on the data availability and the local characteristics but also on the expert opinion and the problem statement. Finding common ground is near impossible, as these aspects are highly case-study dependent. However, the methodologies used for suitability mapping of MAR could potentially be synthesized. Rahman et al. [8] developed a GIS-based tool for MCDA site selection analysis and structured the methodology of site suitability mapping, making a first effort to standardize the GIS-MCDA methodology for MAR site selection.

This study continues the work of Rahman et al. [8] to structure and simplify the decision-making process for MAR suitability mapping. Our work is based on the knowledge generated from a previously published review on GIS-MCDA application for MAR suitability mapping [11]. Findings from the review were taken into consideration to design and implement the web tools presented in this paper. Two related tools were designed to help with the standardization of the MAR mapping process. All data collected from the review were implemented in a web-based query tool that makes the information easily accessible. From the review, the most frequently used methodologies for map generation were determined and included in a web GIS tool. This tool takes a systematic approach, engaging decision-makers in the MCDA process in a structured way. Links between the web GIS tool and the query tool support the decision-making process as they readily depict GIS criteria, criteria weighting, and MCDA methodologies.

While the previous work was dedicated to structuring the suitability mapping process, the present paper focuses on the development of user-friendly tools and their web-based implementation. Thereby, this work aids the decision-makers in undertaking a standardized mapping procedure and, thus, can help to increase the reliability of the method application for the generated maps. As the tools are web-based, they enable the collaboration of multiple stakeholders, thus, potentially improving the decision-making process as well as facilitating easy sharing of results. The web-based nature, as well as the open-access availability, thrive to enhance the usability of the tools, a particularity distinguishing them from existing desktop solutions.

#### **2. Implications from Reviewing GIS-MCDA Studies**

This section is based on an already published review on MAR suitability mapping [11] and focusses on analyzing the relevant parameters and methods for the comprehensive approach to GIS-based suitability mapping by Rahman et al. [8]. Their approach divides suitability mapping into four main steps. It follows the scheme of (a) problem definition, (b) screening of suitable areas (constraint mapping), (c) suitability mapping including the classification of thematic layers or criteria, standardization, weighting of the criteria, and layer overlaying by decision rule, and (d) sensitivity analysis.

Problem definition is the basis of choosing relevant GIS maps and weighting their importance for solving the problem statement. As this part of the approach is highly problem-specific, no general statement or implications could be formulated from the review.

Constraint mapping identifies the parts of the study area that are not suitable for the application of MAR or need to be excluded from the analysis as they are, for instance, natural reserves or private land. This is achieved through threshold values for the GIS criteria and by applying Boolean logic to clip respective areas from the final map. Half of the analyzed studies used constraint mapping as a tool to exclude unsuitable areas. As this methodology is widely used, options to constrain single GIS datasets or complete areas from the resulting map were envisaged for the web GIS tool.

Suitability mapping is the core of the MCDA process as it ranks the study area based on its suitability for the application of MAR. This step comprises the standardization of GIS maps, the assignment of weights to every map, and the combination of the weights and the standardized maps by decision rule. The most commonly used weight assignment methods are the rating method, the ranking method, the multi-influence factor (MIF), and the pairwise comparison. The rating and ranking methods are very simple methods comprising manual weight assignment on a predetermined scale [17]. MIF is a graphical weight assignment method where linkages between GIS datasets are drawn, and the weights are calculated based on the number and the importance of linkages between the criteria [18,19]. Pairwise comparison is the most used method for GIS-MCDA in the context of MAR. The weights are calculated through a matrix-based comparison of pairs of criteria [20]. The methods range from simple (rating method) to more complex (pairwise comparison). The advantage of the simple methods lies in the easiness of use, whereas the complex methods, such as pairwise comparison, offer a coefficient indicating the consistency of the decision-maker's choices. To account for the advantages of both the simple and more complex methods, the rating and ranking method, MIF and pairwise comparison were chosen to be incorporated into the web tool.

The decision rule states how the standardized datasets and their weights are combined to obtain the suitability map [17]. This integration can be based on threshold values (Boolean logic) or more elaborate integration rules, such as weighted linear combination (WLC). WLC comprises the summation of the weighted and standardized criteria and is the most commonly used decision rule. It has been developed further to its derivative analytical hierarchical process (AHP). AHP is the more structured approach that categorizes the GIS maps into hierarchical levels before aggregating und summing up the weighted criteria. It is used to solve more complex decision problems. The two most used methodologies, WLC and AHP, were chosen to be incorporated into the web tool.

To verify the obtained map, a sensitivity analysis should be conducted. It is used to display the effect of different standardization and weights on the final suitability map and indicates the robustness of the obtain suitability map [16]. While this is an important factor for strengthening the reliability of suitability maps, only 21% of the reviewed studies conducted this step.

All reviewed studies used desktop GIS software for their analysis. Some studies created their own tools for the analysis, e.g., Rahman et al. [8] developed the GIS-based Gabardine desktop decision support system. Other studies used tools available through standard software, e.g., an AHP tool has been implemented as an extension for ArcGIS [21], which has been applied by Anane et al. [22].

#### **3. Materials and Methods**

The developed tools are embedded into the INOWAS online platform (https://inowas.com). The INOWAS platform is an open-source collection of empirical, analytical, and numerical web-based models focusing on the planning, management, and optimization of MAR applications. The INOWAS platform and the tools are accessible with any state-of-the-art web browser. The platform works account-based, enabling the user to store and share their work. One key feature of the platform and its tools is the intuitive graphical user interface, guiding the users through the application of the tools.

The technical infrastructure of the platform is based on three components: the CLIENT, which is the user's internet device and browser, the SERVER, which is a standard Linux Server, and the WORKER, which is a Linux cluster with connection to the server. These components communicate with each other via a TCP/IP connection using standard protocols such as HTTP/HTTPS. The REST interface developed by the INOWAS group specifies the individual API calls and their function between the components.

For MAR suitability mapping, two tools were developed for the INOWAS platform. A query tool to filter information from the database that resulted from reviewing related studies and a web GIS system to obtain the suitability maps by following an integrated workflow.

The web-based query tool was designed to grant easy access to the information gathered from the GIS-MCDA review, namely data on different MAR related aspects as well as the MCDA methodology used in the studies. It is based on a pivot table approach. Decision-makers can sort, average, or sum up the database content by creating tables and graphs. Filters can be used to make specific queries, e.g., search studies for a certain MAR technique. The review tool is developed with ReactJs and is based on an open-source 3rd party project (https://github.com/nicolaskruchten/pivottable).

The INOWAS platform and web GIS system interfaces were created with ReactJS and Semantic UI (https://semantic-ui.com/) using some open-source 3rd party projects. Geodata is displayed in leaflet maps (https://leafletjs.com) and uses open street map layers as base layers (https://www. openstreetmap.org). Charts are displayed with ReCharts (http://recharts.org/en-US/), sliders with rc-slider (http://react-component.github.io/slider/) and network diagrams with visJS (http://visjs.org/). Raster calculations are done in JavaScript, using the mathJS library (http://mathjs.org). Users can save and share projects via a connection to the INOWAS backend server and their API entry points.

The full code of the INOWAS platform and the tools for MAR suitability mapping with version history is accessible through GitHub: https://github.com/inowas/inowas-dss-cra.

#### **4. Results**

All tools developed are available through https://inowas.com/tools/ using the "Start using now!" button. The tools are open-access and free of charge, but user registration is required. The tools can be accessed through the personal dashboard, which shows all available tools and stored projects. The projects can be shared with other users or can be made publicly available so that various users have access to the project and can edit it.

#### *4.1. Database Query Tool*

The database query tool is listed as tool T04 in the toolbox of the dashboard. Its basis is a database with information accumulated from the reviewed GIS-MCDA studies. The tool enables the user to research MAR specific information, such as the MAR type used, the water source, the objective of MAR application, or the location of the study. The main information stored is focused on MCDA related data, for example, the number and type of criteria used in the study, the weights assigned to the criteria, and the criteria standardization. Furthermore, information on weight assignment methods, decision rules, and the use of constraint mapping or sensitivity analysis by the authors of the study has been accumulated.

The different attributes from the database can be chosen to be displayed by dragging them into the column or the row fields of the tool, also enabling the display of combinations of attributes (Figure 1). Each attribute is equipped with a filter function where specific information queries can be chosen through class selection or conditional and numerical operators. The results can be visualized in different forms of tables and heat maps as well as be exported. They can be further modified through conditional and mathematical operators, including "Count", "Count Unique Values", "List Unique Values", and "Sum". The tool is supported by documentation explaining the underlying database and displaying three examples that help to get the user acquainted with the functionalities of the tool.

The database query tool can aid decision-makers at different steps along the MAR mapping procedure. During the problem statement step, it allows the decision-makers to investigate GIS criteria that are used often for certain MAR techniques or certain recharge water resources. During the weight assignment step, the importance that other decision-makers have given to a criterion can be assessed.


**Figure 1.** Interface of the database query tool showing the query for the most used criteria for the suitability mapping of in-channel modifications.

Figure 1 shows an example of the tool where the query was set to find the most used GIS criteria for suitability mapping of the in-channel method. Here, the filter was set to display the number of studies that used each criterion for this MAR method and to visualize them as a table heatmap. The results show that slope was the most used criterion, followed by land use and geomorphology. By adapting the attribute selection, it is also possible to show the 14 papers using the slope criterion. This enables the user to further study them regarding the use of the slope criterion, e.g., how it was standardized or classified in the studies.

#### *4.2. Web GIS for Suitability Mapping*

The web GIS is listed as tool T05 in the toolbox of the dashboard. The user is guided through the MCDA workflow with the help of a systematic approach indicated by the menu in the left column of the tool structure (Figure 2). It follows the workflow introduced by Rahman et al. [8] but excludes the steps of the problem statement and sensitivity analysis. The process starts with (1) the choice of GIS criteria, and continues with (2) the weights assignment, (3) the data upload, the constraint mapping and the reclassification, (4) the additional global constraint mapping, (5) the suitability mapping, and (6) the results visualization. The user is guided through the workflow by small green or orange circles that indicate whether the steps have been completed successfully. In case preceding information is required for the subsequent steps, the link in the navigation menu might be disabled until all requirements have been fulfilled.

A simplified example for MAR suitability mapping was prepared to depict the capabilities of the tool. It comprised suitability mapping for surface infiltration methods in southern Africa with four GIS criteria: geology, soil, land cover, and slope. A comprehensive example for MAR mapping in southern Africa was prepared in [23], including geoinformation on water sources and water demand. In this manuscript, the number of used criteria was reduced to account for better readability of the figures.

**Figure 2.** Interface of web GIS tool, with workflow display on the left and the weight assignment tool multi-influencing factor method. Here, the user can draw arrows between entities representing major or minor influences of a criterion upon another entity.

The problem definition step suggested by Rahman et al. [8] was not directly included in the tool approach. As this is a conceptual, case-specific step framing the project objectives and choosing criteria based on those objectives, it could not be incorporated into a tool. However, the database tool can deliver indications on the criteria choice based on the evaluation of previously conducted studies with similar objectives. In the case of surface infiltration methods, the four aforementioned GIS criteria are the most used criteria according to the database query tool.

Starting the tool, the user must choose whether to use the AHP or the WLC method as the decision rule, the latter being the default method of the tool. In both cases, all GIS criteria to be used for the suitability mapping need to be listed. Here, a link to tool T04 was included, enabling the user to analyze other MAR mapping studies regarding their choice of criteria. For each GIS criterion, the user needs to specify whether the criterion data is discrete or continuous and may set units to support their visualization. If AHP is activated, the main criteria classes need to be defined, and all GIS data need to be sorted by those main criteria classes. For this, a hierarchical tree is used and graphically presented. Weight assignments are done for the main criteria and ensuing for each branch with the respective sub-criteria. AHP is useful for complex problems with many GIS criteria. Dividing the procedure into separate branches simplifies the weight assignment step.

The weight assignment step allows performing any number of weight assignments through the integrated methods: ranking, rating, MIF, and pairwise comparison. This enables the user to try out different methods and compare the results, choosing the method and the subsequent weights most suitable for their project. All weighting methods are implemented in a visually appealing and user-friendly way.

For the MIF method, the user must draw connections between the different criteria and the project itself. Those connections represent the dependencies between the entities (Figure 2). The user can adapt the direction and strength of the connections. This information is then used to calculate the associated weights. In the example shown in Figure 2, the connections drawn between the four used criteria and the project itself correspond to weights that emphasize the importance of slope and geology.

For the ranking method, intuitive approaches, such as arrow buttons or the drag and drop method, help to order the criteria in a list, starting with the most important criterion and ending with the least important one. All ranks are summed up and then converted to weights in relation to the other criteria input values. The rating method is the most basic method of all, enabling the user to choose their own weight for each criterion. The fourth available weighting method is the pairwise comparison. On a predefined scale, the criteria are compared to each other by moving a slider towards one side, indicating the importance of a criterion compared to a second criterion (Figure 3). Based on these pairwise preferences, the criteria weights are calculated. The sliders indicate a preference of geology and slope over the other criteria, which is then depicted in the higher resulting weights. From the weighting choices, a consistency coefficient is calculated. The coefficient indicates the consistency of the user's preferences, and if it surpasses a threshold value, the users are asked to re-check their pairwise comparison choices. Again, a link to T04 is provided to the user, enabling the investigation of previously conducted studies regarding their criteria weighting choices.

**Figure 3.** Interface of web GIS tool, with weight assignment tool pairwise comparison where user can set criteria preferences via sliders and indication of the robustness of the decision is given via consistency ratio.

During the third step, "Criteria data", data upload, constraining, and reclassification are performed for each GIS dataset. At the beginning, the final grid size of the project must be set. Then for each criterion, a GIS raster file must be uploaded. Uploaded raster files may be resampled through nearest-neighbor interpolation if the raster size differs from the final grid size. Constraints can be set in the second step, indicating the criterion classes that will not be used for the suitability mapping process. With discrete data, the automatically created classes can be disabled. For continuous data, ranges not to be considered in the final calculation can be defined by Boolean logic. Afterward, reclassification and normalization of the criterion data are performed. For continuous criteria, classes can be defined by indicating the minimum and maximum value and their respective reclassified value or by using a reclassification function. The reclassified values should be between zero and one. For discrete datasets, every existing criterion class is assigned a normalized value, with the possibility of several criterion classes forming one normalized class. All classes can be named and given a color for geo-visualization. Finally, the results of the reclassification are displayed. It is possible to switch between the original, reclassified, constraint, and resulting data.

In the fourth, optional step, global constraints can be set by drawing polygons in a project area, which will then be disregarded from the final calculation of the suitability map.

For the fifth step, "Suitability mapping", the user must select one of the obtained weight assignment calculations. Then, the calculation is performed, combining all information on constraint mapping, criteria standardization, weighting, and decision rules. For the resulting map, the suitability classes can be redefined or left as default. Finally, the suitability map is displayed and can be downloaded as a text file for further processing in other GIS software (Figure 4). The map for southern Africa shows the geographic distribution of areas that, based on the criteria evaluated, are more suitability for the implementation of MAR schemes.

**Figure 4.** Interface of web-GIS tool, showing the final suitability map with redefined suitability classes.

The tool is supported by documentation explaining the tool functionalities as well as the underlying concepts and methodologies (https://inowas.com/tools/t05-gis-mcda/). One example case study is used to help to get the user acquainted with the functionalities of the tool. It is incorporated into a tutorial that provides the user with a step-by-step guide on how to generate a simple suitability map, also providing an example dataset (https://inowas.com/#tutorials, Tutorial 4).

To validate the correctness of the web-tool and the methodology incorporated, a case study was both prepared with ArcGIS and our web tool for comparison [23]. Minor differences occurred with a normalized root mean square error below 0.1. These differences were not attributed to the MCDA methodology but to the resampling during the data upload step. The nearest-neighbor algorithm was found to shift the raster by one half-pixel [24]. This caused a slight divergence between the map obtained from this tool and the map obtained through ArcGIS.

#### **5. Discussion**

This study indicates a trend in MCDA methodology applied for MAR suitability mapping, namely constraint mapping, suitability mapping by using pairwise comparison, and WLC or AHP, and less often a subsequent sensitivity analysis. Based on these findings, we designed an open-source web tool to guide the user through the MCDA process. We included several weight assignment methods so that next to pairwise comparison, the decision-maker can use rating and ranking method as well as MIF. While we kept those methods for their simplicity or advantage in visual decision-making, the combination of pairwise comparison with AHP must be highlighted as the methodology with the highest increase in usage [11] and the most benefits.

AHP offers a clear, systematic procedure that represents all aspects of the problem statement enforcing robust decision-making. In combination with pairwise comparison, the decision-maker must only give priorities considering two criteria at a time. Through a designated index, the methodology offers an indication of the consistency of the decision-maker's choices. The shortcomings of this method include the inability to include threshold values, no direct measure to assess the robustness of the criteria standardization values, and possible rank reverse issues [25]. Nevertheless, it can be asserted that, currently, AHP, in combination with pairwise comparison, is the state-of-the-art methodology for MAR suitability mapping. This can be further underlined by its frequent use in other environmental MCDA studies [25–27].

The established web tools cover all aspects of the suitability mapping process apart from the sensitivity analysis. A structured GIS-MCDA process should include a sensitivity analysis as it can help to generate more robust and reliable suitability maps. A future prospect of the web tool will be the integration of the sensitivity analysis process, which is relatively more complex and computation-intensive and, thus, was not included in the workflow yet. Furthermore, it is planned to improve the tool by increasing the possible maximum resolution of projects and by providing more flexibility for handling input raster data as well as supporting vector files.

The web-based implementation of the tools offers advantages over standard desktop GIS solutions. The user does not require advanced GIS-specific knowledge and does not need to install any GIS software unless pre-processing of the GIS data is necessary. Nevertheless, the platform provides an interface for data exchange with standard GIS software to allow for pre- or post-processing with conventional desktop-based software. The entire system works in a standard web browser, with no specific system requirements. The input data, as well as the resulting maps, are available online and can be easily shared among stakeholders allowing for collaboration on the mapping procedure as well as flexible sharing of the results. The web tools can be run in two operating modes (private and public), which offers a high degree of flexibility in data sharing as well as adequate privacy. The entire workflow is very transparent, with the possibility to revise and reverse steps. It provides a pre-defined workflow for inexperienced users while offering a comparison of several MCDA methods for more advanced users.

While MAR suitability mapping can be seen as a viable source for strategic water management, it is not a sufficient technology to point down to actual locations for MAR implementation. Maps can deliver an indication of areas of interest, but those need to be further assessed by numerical modeling [10,28] or on-site measurements of the local hydrology and hydrogeology [13,29,30]. For the numerical modeling and MAR scheme design and optimization, the INOWAS platform can be used as well, as it further offers numerical groundwater flow modeling tools as well as algorithms for MAR scheme optimization.

#### **6. Conclusions**

The developed web tools can help planners of MAR sites by increasing knowledge on MAR suitability mapping as well as by engaging in the MCDA processes in a structured and cooperative way. The tools are further envisioned to aid in capacity building measures as well as the education of water practitioners by accumulating knowledge on GIS-MCDA in the context of MAR and translating them into easy-to-use web tools. The clearly outlined process of map generation enforces standard methodology and can help to generate maps that are comparable due to a common methodological approach. Since MAR mapping has been increasingly used in recent years, the quality of the maps produced should be critically evaluated, and analysis and categorization of the methodologies used is one first step to improve the reliability of the maps.

While the tools can outline the map generation process, they cannot standardize one of the main sources of uncertainty—the datasets and respective weights assigned. Putting individual choices into perspective with similar studies retrieved from the database tool is a step towards decreasing the subjectivity of their weighting and standardization process. However, the problem statement and its specifics define the importance of a GIS dataset for each individual case study. Thus, the weights assigned to the criteria cannot be defined by rules. Furthermore, data availability and quality are major constraints in the mapping process. Thus, the tools at hand standardize and simplify the MAR suitability mapping process but cannot substitute the decision-maker's expertise in choosing relevant datasets and their importance for the specific study.

**Author Contributions:** J.S. and J.P.B.V. designed the overall study; J.S., F.V.L. and J.P.B.V. conducted the review, J.S., R.S. and F.V.L. devised the tools; R.S. and R.J. implemented the tools; J.S. and R.S. wrote the paper; C.S. contributed with discussion and advice.

**Funding:** This study was supported by the German Federal Ministry of Education and Research (BMBF), grant No. 01LN1311A (Junior Research Group "INOWAS").

**Acknowledgments:** We acknowledge the financial support of the graduate academy of TU Dresden as well as the Open Access Funding by the Publication Fund of the TU Dresden.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Potential Benefits of Managed Aquifer Recharge MAR on the Island of Gotland, Sweden**

**Peter Dahlqvist 1,\*, Karin Sjöstrand 2,3, Andreas Lindhe 3, Lars Rosén 3, Jakob Nisell 1, Eva Hellstrand <sup>1</sup> and Björn Holgersson <sup>1</sup>**


Received: 6 August 2019; Accepted: 16 October 2019; Published: 17 October 2019

**Abstract:** The Island of Gotland (3000 km2), east of mainland Sweden, suffers from insufficient water availability each summer. Thin soils and lack of coherent reservoirs in the sedimentary bedrock lead to limited reservoir capacity. The feasibility of Managed Aquifer Recharge (MAR) is explored by identifying suitable areas and estimating their possible contribution to an increased water availability. MAR is compared to alternative water management measures, e.g., increased groundwater abstraction, in terms of costs and water availability potential. Results from GIS analyses of infiltration areas and groundwater storage, respectively proximity to surface water sources and surface water storage were classified into three categories of MAR suitability. An area of ca 7700 ha (2.5% of Gotland) was found to have good local conditions for MAR and an area of ca 22,700 ha (7.5% of Gotland) was found to have moderate local conditions for MAR. These results reveal the MAR potential on Gotland. The water supply potential of MAR in existing well fields was estimated to be about 35% of the forecasted drinking water supply and 7% of the total water demand gap in year 2045. It is similar in costs and water supply potential to increased surface water extraction.

**Keywords:** MAR; groundwater; mapping; Sweden; decision-support

#### **1. Introduction**

The Island of Gotland (3000 km2), situated in the Baltic Sea 100 km from the mainland of Sweden (Figure 1a), suffers from insufficient water availability to supply the ever-increasing demand from society, especially during the tourist season (June–August, [1]). The annual precipitation on the island (ca 550 mm/year) is sufficient to cover a forecasted increase in water demand. However, intensive drainage of arable land, thin soil layers, and relatively impermeable rock lead to precipitation run-off and limited reservoir capacity in both surface water and groundwater reservoirs [2]. The already constrained water supply will be further aggravated in the future because the total water demand on the island is estimated to increase by 40% by 2045 [1]. The current water resources on the island will not meet this projected increase in demand. A high availability of water during the winter and a high demand for water in the summer makes MAR a suitable way to increase the water resources. Due to these factors, it is important to investigate the potential for MAR on Gotland.

The bedrock on Gotland consists mainly of Silurian limestone and marlstone, which represents the upper part of a 250–800-m-thick sequence with Palaeozoic sedimentary rocks overlying the crystalline basement [3]. The Quaternary overburden is generally thin (less than 2 m) and is largely composed of till and postglacial sand deposits. The relief is low, and the highest point is 82 m a.s.l. The main land uses are agricultural and forestry. The main aquifers on Gotland are situated within the bedrock where cracks, fractures and dissolution cavities store and transport the groundwater. Nevertheless, the soil layers have an important role to play because areas with soil, especially sand and gravel, act as infiltration and storage systems for the bedrock aquifer. The island may be considered one large aquifer, but with groundwater divides (created by relief) producing 7 (sub)aquifers, according to the European water framework directive. Saline groundwater is a problem because relict saltwater occurs under the entire island at a depth of 20–100 m b.s.l.

The total water demand on the island is estimated to increase by more than 40% by 2045, with increases of 30% in tourism, 20% in domestic demand, 20% in animal keeping, 15% in industry, and 100% in irrigation [1]. To enhance water resource security, and close the water supply and demand gap, several alternative water management measures are being examined. Managed Aquifer Recharge (MAR) is one of them. Today, the public water supply on Gotland relies on 14 well fields, two surface water catchments and a desalinisation plant. MAR is currently not used in any public water supply on Gotland, but may play an important role in the future if suitable areas can be found. Conversely, on the Swedish mainland, MAR has been in use for over 100 years, and accounts for approximately 20% of the public water supply [1,4]. MAR can be explained as the intended recharge into and storage of water in an aquifer [5]. It may be used to increase water security for uses including drinking water supply, irrigation, preventing saltwater intrusions, as well as providing environmental benefits [6]. MAR is widely distributed and applied on various scales around the globe, as well as in Europe [7]. The water source can be of varied origin, e.g., river water, seawater or sewer water. In some cases, there is a need for pre-treatment before groundwater recharge to minimize the risk of pollution or aquifer clogging [8]. The recharge can be made by spreading methods in areas with high infiltration capacity; by deep infiltration direct into the aquifers via wells; or as induced infiltration due to withdrawal [9]. Since the different MAR types are suitable for different conditions within hydrogeological settings (e.g., confined or unconfined aquifers), treatment opportunities, and land use, the selection of suitable recharge sites is crucial [10,11]. In this study, we focus on areas suitable for recharge through infiltration basins and natural conditions for storage. Conditions for well or induced infiltration are expected to be of minor importance because of the geomorphology, geology and hydrology of the island. Hence, the suitability of those MAR types is not investigated in this study.

To prioritize between alternative measures to improve water resources security (e.g., increased groundwater abstraction and desalination), useful decision support is needed. GIS-MCDA (Geographical Information System Multi Criteria Decision Analysis) [12] is a regularly applied method in MAR suitability assessment [13]. There are several possible criteria for mapping MAR suitability, the three most common being aquifer storage capacity, geomorphology and soil [13]. There are also concerns on limitations and discussion on the uncertainties of these GIS approach made visible by e.g., [14]. Although a GIS analysis will show where MAR might be successful, field work and numerical modelling will be important tools for increasing the success of MAR [15–19]. The economic assessment of MAR is an important question that has been studied previously [20–22]. The aim of this paper is to explore the feasibility of MAR on the island of Gotland by: identifying potential areas for MAR in proximity to a fresh water source and which are available for recharge; estimate the possible increase in groundwater recharge and groundwater extraction at existing wellfields; and compare MAR to other alternative measures in terms of costs and water availability potential.

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

#### *2.1. Materials*

Mapping of potential MAR locations on Gotland is based mainly on the existing data listed as follows: Intensive hydrogeological investigations from two campaigns of airborne transient electromagnetic surveys (2013–2015, SkyTEM resistivity measurements along flight lines covering 30% of the island, with 200 m spacing and geophysical soundings every 30 m) [23,24]; a site-specific overview in 2016 of groundwater catchments of the 14 existing well fields that showed that 30% of those well fields had a favorable geology and hydrology for spreading and induced infiltration-based MAR types [2]; 3D geological and hydrogeological models (Geoscene 3D by I-GIS) of the entire island (2015–2019, available online in 2020) based on resistivity models from the SkyTEM survey; existing geological information such as bedrock and soil maps (regional scale); seismic profiles and information from water wells; and national scale mapping (2018) of groundwater recharge and storage capacity [25] provide comprehensive additional data for assessing the potential of MAR on Gotland.

**Figure 1.** GIS maps and results from GIS analysis to identify areas with potential for MAR. (**a**) Map of Sweden and location of Island of Gotland. (**b**) Groundwater storage capacity (mm/year) [25]; (**c**) Ratio of groundwater recharge/storage capacity; only values below 1.0 are shown (GS); (**d**) Closed depressions (>1 ha) in the bedrock with no contact with the Baltic Sea (CD); (**e**) Areas from the lithological 3D model, blue: areas >4 m sand and/or gravel (IA) suitable for infiltration, red areas >4 m till and/or clay (ST) suitable for construction of surface storage dams; (**f**) Surface waters, i.e., raw water source (S).

#### *2.2. Methods*

A GIS-based (Boolean logic) approach was used to find suitable locations for MAR systems (in this study focusing on areas suitable for recharge through infiltration on the surface) on Gotland. No parameter weighting was included. Several mapping projects assessing suitability for MAR have

been made globally [13]. GIS was used for analysis of 5 surface and subsurface datasets, which are presented below. Three of these (1–3) concern the MAR location, and two (4 and 5) further explore the MAR location to determine if there are local sources of water supply for the MAR site.


In addition to the above analyses, two further GIS analyses were completed with a combination (overlap with no priority weight) of data on favorable areas for infiltration (Figure 2a) and surface water source and storage (Figure 2b).

**Figure 2.** Maps of the combined GIS analysis. (**a**) Areas suitable as infiltration areas and for groundwater storage; (**b**) Areas close to a natural source (intermittent stream, perennial river, lake) and for construction of surface storage areas (dams). Data has been reclassified into classes; lower numbers more suitable.

To make estimations of possible increases in groundwater recharge and groundwater extraction through MAR, in relation to current well fields, favorable groundwater catchments were identified [5]. This was done using a GIS-aided and analytical approach. The values presented here for potential infiltration rate and increased withdrawal volumes are estimations with an inherent uncertainty. The estimates are based on access to surface waters of adequate size, presence and sufficient thickness of permeable sediments, the possibility of creating dams with local material, and current withdrawal capacity.

MAR was compared to other alternative measures (i.e., increased groundwater abstraction, enhanced water reuse for irrigation, increased surface water extraction, metered leak detection and desalination) in terms of annual water availability potential and economic viability. The measures were selected for inclusion in the analysis based on the outcome of a multidisciplinary stakeholder workshop, in which the participants were asked to identify measures with potential to improve the water resource security on the island. The comparative method was based on marginal abatement cost curves [27,28], including cost–benefit and cost-effectiveness analyses. The costs of the measures included investment costs, operating costs and cost savings. The measure costs were described by present values (PVs) [29], analyzed with a 3.5% discount rate over the 27-year time horizon from year 2019 to 2045 (corresponding to the current water plan period applied at Gotland) and based on 2018 prices. The NPVs were then expressed as equivalent annual costs (EACs) in SEK (Swedish krona) per year [30]. A theoretical maximum level of implementation was assumed for each measure category, except for desalination, which instead was based on estimates of one new desalination plant. The estimated cost and water availability input variables were based on a combination of literature data, expert judgements and GIS-based analyses. The unit cost of each measure was calculated as the ratio of the measure's EAC and annual water availability potential.

#### **3. Results**

#### *3.1. Mapping of Suitable MAR and Source Areas*

Results are described in three parts; areas where the geology is favorable for infiltration and/or groundwater storage (IA (Infiltration Areas) + GS (Groundwater Storage) + CD (Closed Depressions)); areas with proximity to a surface water source and/or suitable for surface water storage (S (Source) + ST Surface water storage)); and areas where these overlap (IA + GS + CD + S + ST).

#### 3.1.1. Infiltration Areas and Areas for Groundwater Storage (IA + GS + CD)

Three data sets regarding the possibility for an area to be suitable for artificial groundwater recharge were combined into a raster set (Figure 2a). The resulting raster was divided into three classes: areas with three (class 1), two (class 2) or one (class 3) of the included favorable attributes. A raster with class 1 has a good potential for both infiltration and groundwater storage. Areas with class 2 or 3 are less suitable because one or two of the included raster sets (IA, GS, CD) is absent. Class 1 is present in 2068 ha (0.7% of Gotland), class 2 in 14,437 ha (4.8% of Gotland), and class 3 in 43,453 ha (14.4% of Gotland).

#### 3.1.2. Source and Suitable Areas for Surface Water Storage (S + ST)

The source for artificial groundwater recharge in this investigation is natural surface waters from streams, rivers and lakes. To regulate and decide where there are favorable infiltration conditions there is also a need for a seasonal surface storage, in this investigation in the form of man-made dams. In this study, we do not discuss present land use and slope—two factors that might influence outcomes—but we regard them as having a minor influence (see Discussion) on Gotland. The combination of these sources with possibilities of constructing storage capacity is shown in Figure 2b. The resulting dataset is divided into three classes: (1) areas with both a nearby raw water source and good conditions for storage in surface dams (10,271 ha, 3.4% of Gotland); (2) areas with good conditions for dams, but more than 0.2 km (smaller intermittent streams) or 0.5 km (perennial rivers and lakes) from a source (10,345 ha, 3.5% of Gotland); and (3) areas with a source but the distance to a suitable storage area (dam) is more than 0.2 km (smaller intermittent streams) or 0.5 km (perennial rivers and lakes, 112,468 ha, 37% of Gotland).

3.1.3. Areas with Combination of Infiltration, Groundwater Storage, Source and Surface Water Storage (IA + GS + CD + S + ST)

To narrow down the selection of promising areas for MAR a raster set with different combinations of IA + GS + CD + S + ST is shown in Table 1 and Figure 3b. These are the best-adapted areas for MAR construction (recharge through infiltration) based on this study. There are a few stream catchment areas that appear more promising to work in. Please note that because the raster set with surface water source/storage is made with a buffer there will be areas with overlapping datasets (Figure 3a).

Class 1 (Table 1) is not present on Gotland in this analysis. The reason is in the construction of the analysis these conditions cannot exist in the same raster (ha scale). Class 2 is present in ca 7719 ha (2.5% of Gotland). Areas in this class have good conditions for successful MAR. Class 3 is present in 22,710 ha (7.6% of Gotland). In these areas, there are probably moderate conditions for successful MAR. Class 4 is present in 2765 ha (0.9% of Gotland). In these areas, there is no surface water source within the chosen distance.


**Table 1.** Synthesis from final step in GIS analyses. GS = Ratio >1.0 groundwater recharge/storage capacity; CD = Closed depression; IA = sand/gravel >4m; S = Raw water Source; ST = Storage. Location of mapped classes see Figure 3b.

**Figure 3.** Maps from combined GIS analysis. (**a**) Combination of overlapping datasets from Figure 2a,b. Please note that since the raster set with surface water source/storage is made with a buffer there will be areas with overlapping datasets. (**b**) Map with the best areas for infiltration/groundwater storage and source/surface water storage. Classes are defined in the text and summarized in Table 1.

#### *3.2. Estimation of Increased Groundwater Recharge and Groundwater Extraction at MAR Favorable Groundwater Catchments in Use Today*

The presented values for potential infiltration rate and increased withdrawal volumes are estimations with considerable inherent uncertainty. Furthermore, there are large differences between the well fields. Site-specific conditions at existing abstraction areas differ due to, e.g., hydrogeology, number of wells, local water supply demand, quality, abstraction volume, etc. The abstraction volume in different well fields varies between 25 and 3300 m3/day. There are also differences in estimation of infiltration capacity between 50–660 m3/day (mostly based on local infiltration capacity) and the estimated increase in abstraction volume between 25–330 m3/day (also includes the possibility of adding new wells). The percentage of predicted increased abstraction volume versus mean abstraction

volume varies between 20–300% for estimated infiltration volume, and 10–150% of increased abstraction volume, in comparison to numbers for each well field, respectively.

#### *3.3. Comparative Study of Alternative Measures*

Figure 4 shows a marginal abatement cost curve for the analyzed measures, in which each measure is represented by a bar showing its unit cost (bar height) and annual water availability potential (bar width). As displayed, increased groundwater extraction and desalination had the largest potentials to improve water availability on Gotland. The water availability potential of desalination can be much larger, but the calculations here were based on assumptions of one new desalination plant. Increased groundwater extraction was associated with the lowest costs per cubic meter water provided, whereas desalination was associated with the highest unit costs. One reason for the high unit cost of desalination was the long pipelines needed to reach demand centers. In this comparative analysis, MAR was limited to groundwater recharge in the municipality's existing well fields. Hence, the water availability potential of MAR on Gotland may be significantly higher when not constraining the analysis to those areas. The unit cost of MAR in existing well fields on Gotland was in the same range as that for increased surface water extraction.

#### **4. Discussion**

The presented data and analysis represent an early stage of mapping of MAR areas (focused on spreading methods in areas with high infiltration capacity) and estimates of potential and feasibility of this type of MAR on Gotland. Our project also includes mapping of good local conditions for a local source, e.g., water supply, which we believe further increases the utility of these study results for water management on the island. Mapping of suitable MAR areas with GIS is a widely used method [13]. There are uncertainties in both data and accuracy in the analysis, but within these limitations there is now a detailed picture on possible MAR and source areas which can be used by the municipality, farmers, and other stakeholders. Concerns regarding the limitations of these GIS analyses have been raised by, e.g., [14], who suggested that the use of sensitivity analyses of the factors used for MAR

feasibility studies. In 2016, SGU made a first attempt to apply an overall assessment of the possibilities of MAR at existing abstraction areas [2]. Results from the present study can now improve the prior prioritization between MAR and alternative water measures for the island. The results should also be addressed on more than the water quantity. Increased groundwater recharge by MAR may influence the quality of the groundwater, e.g., by dilution effects, but also change the salinity levels in the bedrock aquifers. The municipality can analyze the results and compare with areas that are not used today but where the presented method shows good potential for MAR. There will be need for validation with numerical modelling and field tests to increase the strength and substance of the results, as shown in, e.g., [15–19].

This paper does not use any restrictions on unsuitable areas due to land use, a criterion that has been used in several analyses [13]. This can preferably be made by water management authorities, who are more suited to deciding between conflicting interests. Mainly because of the low relief of Gotland, the often-used parameter "slope" [13] is also not used here. A potentially much more important criterion is "closed depressions". This criterion uses the bedrock surface and the slope to calculate where there are bedrock depressions where water has more time to infiltrate and be stored as groundwater.

Even though there is a clear picture of the best areas from the perspective of infiltration/groundwater storage and source/dam, other areas can also take advantage of the presented data. Local conditions—for example, where the distance to the public water network is long—can make MAR solutions profitable in those more remote areas. The outcome of favorable areas (from 2000 ha (20 km2) up to at most 43,000 ha (430 km2)) in the analysis of infiltration areas, source areas and the combination of these should be compared to the area of the island (3000 km2). The designated areas constitute only a small part of the island, and therefore care must be taken so that these are not destroyed by over-exploitation.

The degree of detail in the results is determined by the available data sets. SGU is working on an update on a few of the data sets, which will improve the certainty of the result. There are also several ongoing and future investigations associated with some of the data sets. For example, a few of the designated areas with a closed depression may be particularly suitable areas for groundwater dams. This is also the case with depressions that are not completely closed, not used in our investigation, and hence an interesting subject for future analysis for the viability of this MAR technique on Gotland. A groundwater dam is a man-made structure that obstructs the natural flow of groundwater and thereby can store larger quantities of water in the aquifer [31]. The results from this study are not validated by field studies. The data sets are; however, delivered to the local authorities for water resources and water management, and for analyses. Once the data sets are updated, the results may be integrated into the hydrogeological 3D model and calibration of parameters from field studies may improve future work with MAR on Gotland. The presented data sets may be tested in the newly developed online tools for suitability mapping, e.g., https: //dss.inowas.com/tools [32,33]. The resulting suitability maps will be shared at the international MAR portal (https://apps.geodan.nl/igrac/ggis-viewer/viewer/globalmar/public/default).

The economic aspect of implementing MAR systems to improve potable and agricultural water supply has previously been investigated in different parts of the world (e.g., [21,22,34]). The associated capital costs are highly system specific, influenced by, e.g., hydrogeological, socioeconomic and legal factors [35]. As the economic analysis of MAR in this paper was based on rather small-scale complimentary infiltration of surface water at existing municipal well fields, no additional costs for new wells, treatment plants or pre-treatment were considered. The economic analysis was based on cost estimates associated with infiltration basins, raw water intake and new piping, resulting in a unit cost of approximately 5 SEK/m3. This can be compared to cost estimates for MAR in Spain ranging between €0.08–0.58 per m<sup>3</sup> [20] (approximately 0.8–6 SEK/m3).

The total water demand on Gotland is forecasted to increase by more than 40% by the year 2045 [1]. This will require water currently not available on the island. To make well-founded decisions on how to meet this forecasted demand and concurrently increase the preparedness for water scarcity

situations, thorough decision support is needed. The presented data and analyses can be used to inform decision-making on measures to increase the amount of water that can be recharged on the island. Considering the entire island as a single groundwater aquifer would permit a holistic approach to groundwater management, and the water situation more robust. This is particularly important for the forecasted increase in public water demand but also for individuals, industry and farmers relying on private wells. By comparing MAR to alternative measures, in terms of costs and water availability potential, this paper also provides support in assessments of the measures' economic viability, usually an important decision criterion for municipalities, corporations and individuals alike.

#### **5. Conclusions**

This paper contributes results from analyses of possible MAR areas and their potential to increase water availability on the island of Gotland, Sweden. The method can be used to evaluate the MAR potential in similar areas with the same data sets. The results can be used on different scales, by authorities, the public water producer, farmers, industry and by people with private wells, for improved water resource security and further validated with field tests and more detailed models such as the afore mentioned hydrogeological model. Comprehensive field tests are probably the best for better understanding the problems in general. However, they are time consuming, and only a limited number of sites are available to accommodate the tests. In contrast, the GIS analysis allows us to explore and assess multiple sites with relative ease, yet the validity needs to be carefully checked. The main results are listed below.


**Author Contributions:** Conceptualization, P.D., K.S., A.L., E.H. and B.H.; Data curation, J.N.; Formal analysis, P.D., K.S., J.N. and E.H.; Funding acquisition, K.S.; Investigation, P.D., K.S., A.L. and B.H.; Methodology, P.D., K.S., A.L. and J.N.; Project administration, P.D.; Supervision, L.R.; Visualization, P.D. and K.S.; Writing—original draft, P.D. and K.S.; Writing—review & editing, P.D. and K.S.

**Funding:** The work was funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 754412; Region Västra Götaland; Region Gotland; The Swedish Agency for Economic and Regional Growth; and the Swedish Research Council Formas contract no 942-2015-130.

**Acknowledgments:** The authors would like to thank Mikael Tiouls and Lars Westerlund at Region Gotland for contribution with local expertise in the comparative study.

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).
