Comparative Analysis of Runoff and Evaporation Assessment Methods to Evaluate Wetland–Groundwater Interaction in Mediterranean Evaporitic-Karst Aquatic Ecosystem

: This work compares the applicability of several free-surface evaporation and runoff equations in simulating water level variations of small Mediterranean wetlands. The Amarga and Jarales wetland are two pilot sites with an evaporite-karst genesis located in southern Spain. The water level was continuously recorded in both wetlands, and exhaustive weather monitoring was performed. The combined datasets have permitted quantiﬁcation of the surﬁcial elements of their water budget (precipitation, runoff, and evaporation). Several campaigns of groundwater level measurements were also done to characterize the direction of groundwater ﬂows. The morphometrical analysis of the Jarales wetland was accurately performed based on a LiDAR dataset. A total of 225 limnimetric simulations of the Jarales (90) and Amarga (135) wetlands were performed, combining different evaporation and runoff equations. During the study period, the curve number method, coupled with the Penman equation, reached the Jarales wetland’s best calibrations. The Vardavas–Fountoulakis modiﬁcation of the Penman model ﬁt better with the Amarga wetland record. The obtained results permit speciﬁcation of the water budget of both wetlands during several years and conﬁrm that the groundwater–surface water relationship affects the wetland hydric dynamic to different degrees. Nonetheless, the limnimetric models were calibrated for a short period, including dry years, making it necessary to extend the control period longer and validate the models under different hydroclimatic conditions. Finally, the differences between wetland functioning are explained in a conceptual hydrological model that can be useful for wetland conservation and management of related aquatic ecosystems. The understanding of the origin and fate of water in wetlands permits assessment of how future scenarios would affect hydric functioning and suggests adequate conservation measurements.


Introduction
Wetlands are valuable aquatic ecosystems with a permanent or temporal presence of water at the surface or in the soil [1] that provides a broad list of environmental and human services (e.g., [2][3][4]). Wetlands are frequently related to groundwater, particularly in semi-arid areas, where rainfall presents a seasonal pattern and its annual amount is lower than the evaporation rate [5]. Groundwater-surface water interactions in wetlands are complex, as they are temporarily and spatially variable and sensitive to topographical, geological, and climatic conditions [6][7][8]. During the last decades, some methods have been developed to quantify the groundwater-wetland exchange. However, this remains challenging due to the temporal variation of the water storage and fluxes in wetlands and their sediments' heterogeneity and low permeability [9,10].
Such an approach is common in areas where the groundwater fluxes are difficult to determine with direct measurements. However, its accuracy dramatically depends on the uncertainties in estimating the rest of the budget components: direct precipitation, overland runoff, evaporation, surface water inflow and outflow, and the changes in water storage [10]. As the estimation of volumes (storage, inputs, and outputs) is dependent on the wetland area, the error can be more extensive in wetlands located in low-gradient basins, where the flooded surface can broadly vary (e.g., [9,20]). The uncertainties of the wetland water budget can be partially reduced using precise digital elevation models based on modern tools for elevation data acquisition, such as LiDAR (Laser Imaging Detection and Ranging), which results in more reliable wetland bathymetry [21].
In semi-arid areas, such as the Mediterranean region, wetlands are particularly vulnerable to climate change due to the expected increase in the evaporation rates and the precipitation pattern changes [22]. Thus, defining the hydrogeological functioning of such aquatic ecosystems is vital for their management and conservation. For that reason, it is necessary to develop models that accurately estimate the weight of each of the components in the budget to compute future hydrological changes in this sort of wetlands. In that sense, several comparisons of estimation methods for free-surface water evaporation are found in the literature [23][24][25][26], but equations for overland runoff are more rarely discussed. However, this is a crucial aspect to consider in small wetlands in semi-arid regions [10,27].
The Amarga and Jarales are two wetlands located in the province of Cordoba (southern Spain, Figure 1). The regional government protects them as Natural Reserves, they are part of the Ramsar List of Wetlands, and they are a Special Protection Area (SPA) under the European Union Directive on the Conservation of Wild Birds. The genesis of both wetlands is related to dissolution/karstification processes affecting the underlying evaporite rocks (gypsum, anhydrite, and halite) embedded in a Triassic clayey matrix that constitutes the best part of the so-called Chaotic Subbetic Complexes (CSC) [28]. Despite having a similar origin, both wetlands have different hydrological functioning. Previous research carried out in the area [12,[29][30][31][32] described their hydrological functioning and quantified their water budget components. Most of the research reached similar conclusions regarding the groundwater-wetland interaction, even though they used diverse methods for evaporation and runoff estimation. They performed limnimetric simulations based on utterly different soil parameters, partly because neither the observation period nor the frequency of control/modeling was the same.
This study aimed to reach more profound insights into the hydrologic functioning of evaporite karst wetlands in semi-arid regions, emphasizing the surface water-groundwater relationship. The applicability of several methods used in estimating the water budget component is also discussed, particularly those for runoff and evaporation computation. The final goal was to create a hydrological conceptual model that helps to understand groundwater's contribution to the wetlands and its temporal variations as a tool for environmental management of the dependent aquatic ecosystems.

Geological and Hydrological Settings
The Amarga and Jarales wetlands are located on a plateau of 80 km 2 between the rivers Genil and Anzur (Figure 1a). The outcropping rocks mainly belong to the CSC, formed by a clayey-evaporitic megabreccia in which olistoliths (blocks) of diverse nature and size are embedded, including limestones, dolostones, sandstones, and massive gypsum enclaves, among others [28]. Discordantly over the CSC rocks, post-orogenic deposits are constituted by sandstones, conglomerates, marls, calcarenites, and unconsolidated sediments from Miocene to Quaternary ages (Figure 1b  The landscape is characterized by gentle slopes and plenty of endorheic areas, some of which are occupied by wetlands (Figure 1a) with different hydroperiods. Despite the abundance of low permeability materials in the CSC, the presence of small outflows, wells, and boreholes evidence groundwater flows in the media [30,34]. Moreover, the plateau is drained through three brine springs emplaced by the Genil and Anzur Rivers (Figure 1a), between 247 and 300 m above sea level (asl) [35]. Dissolution and karstification processes promoted the subsidence and collapse of the terrain, giving rise to closed-depression landforms. Furthermore, based on the alignment of fractures and slopes, previous research suggests that tectonics could also play an essential role in the development of karst depression [29,36].
The Jarales wetland is located in the central part of the area (Figure 1c) at 406.6 m asl and has a flat-bottom basin. Its average flooding surface is 5 ha, while the maximum water level recorded is about 3 m. Its flooding pattern is irregular because it commonly gets dry during the summer of average years (from the standpoint of rain amount) but not when the hydrological year is considered wet. On the contrary, no flooding is observed during dry years, not even between autumn and springtime, when most of the annual rainfall occurs. The Amarga wetland is located 3 km NW of the Jarales wetland (Figure 1a), at an altitude of 366 m asl. It has a deep basin and a permanent hydroperiod. Its water stage largely varies from 2.5 m to nearly 6 m high, and its average flooding area extends 4 ha.
The watershed of the Jarales wetland extends 90.4 ha over post-orogenic deposits (Miocene-Quaternary) outcrops (mainly marls, sandstones, and conglomerates) but also on clayey-evaporitic rocks of the CSC (Figure 1b,c). The same materials are on the surface of the Amarga wetland watershed (256.9 ha), although it presents a higher proportion of clays of the CSC, as well as dolostone olistoliths. The morphological characteristics of both wetlands' watersheds (Table 1) reveal similar shape characteristics, although Amarga's is relatively less elongated and has a deeper bottom. The protected area around the wetlands is defined by a 500 m buffer around the maximum flooding surface's perimeter, although only the first 50 m have the highest level of protection given by the Natural Reserve ( Figure 1).  [37]: R c = 4πA/P 2 , where A = area of the basin (m 2 ), P = perimeter (m). b Schumm's elongation ratio [38]:

Methods
Four field campaigns of water table measurements in 39 wells and boreholes were carried out under different hydroclimatic situations (May and September 2015, January and March 2016). Those measurements permitted the drawing of schematic water-table contour maps that characterize the direction of groundwater flows in the study area.
Two automatic water-depth data loggers were installed next to respective staff gauges in the Amarga and Jarales wetlands to obtain daily records of water surface fluctuations. The register extended from October 2014 to September 2017 for the Amarga wetland, while it lasted from March 2014 to July 2015 in the Jarales wetland, when it got dry. Additionally, single staff lectures were done in the Jarales wetland until September 2017, when occasional flooding occurred.
The changes in the volume of water stored in a wetland (V) can be approximated as a function of the variations in the wetland stage (h) and the flooding area (A) [10]: Thus, the change in the wetland storage between two heights, h 1 and h 2 , can be calculated given the size of their corresponding flooding areas (A 1 and A 2 , respectively): The flooding surface related to different wetland stages was calculated, considering height intervals of 0.5 m. To that end, a Digital Terrain Model (DTM) of the area, with a resolution of 0.5 m × 0.5 m, was created, based on LiDAR data [40] obtained from the National Geographical Institute of Spain (IGN). This dataset has a mean square error in vertical measures below 10 cm and a nominal data density of 0.5 points/m 2 . This density could be low for the DTM resolution; however, minor deviations are observed when interpolating the DTM, and therefore this issue has low implications on the volume calculation. The point cloud cover was generated from a photogrammetric flight performed in summer 2015 when the Jarales wetland was dry, allowing for the whole wetland basin's morphological characterization. In Figure 3a,b, the considered h values are shown versus their corresponding stored volumes (V) and flooding areas (A). The equations defining those relationships (the filling curve and the hypsometric curve, respectively) were obtained using the software CurveExpert. During the recording period, the Amarga wetland's stage descended below the water table's height registered by the photogrammetric flight. Consequently, LiDAR data do not provide helpful elevation information for assessing the Amarga wetland storage changes during the control period. Nevertheless, García-Ferrer et al. [41] published bathymetric data of the Amarga wetland used by Aljibe Consultores [31] to obtain the filling curve and the hypsometric curve defining the wetland stage relationship with the storage and flooding surface (Figure 3c,d). Additionally, the delimitation of the watershed of both wetlands was created based on the DTM obtained from the LiDAR dataset. To that end, flow direction and flow accumulation covers were calculated, and all of the micro-watersheds flowing toward the wetlands were merged later.
As the drainage network in the study area is not well defined and there are no permanent surface water inflows to any of the wetlands, the only two surficial water inputs considered were direct precipitation and diffuse overland runoff. Direct rainfall inputs were obtained by multiplying the wetland surface at each time and the precipitation was measured at a weather station located in the central part of the plateau ( Figure 1). Surficial runoff was estimated using simple rainfall-runoff models that are commonly used in wetland water budget studies: the curve number (CN) method [16,42,43] and the soil water balance (SWB) [12,22,32,44,45]. The CN method was developed by the US Soil Conservation Service [46]. It considers that runoff starts when a rainfall event exceeds a particular threshold value called initial abstraction (I a ), independent of whether the soil is saturated with water. Then, runoff (Q) depends on precipitation (P) and I a : The I a value depends on land use, slope percentage, and soil texture. In general terms, a low value of I a (nearly 0 mm) corresponds to low permeability soils (clayey), scarce vegetation, and steep slopes. On the contrary, I a has high values (>100 mm) if the conditions are less propitious for runoff generation: permeable soils (sandy), dense vegetation, and smooth relief. The I a values were selected according to the criteria of the Spanish Ministry of Public Works [47], and they are relatively low (10, 15, and 19 mm), since the wetland watersheds have gentle slopes and low permeability soils, occupied mainly by sparse natural vegetation or olive orchards. The agricultural practices related to olive growing include plowing and eliminating other arboreal or herbaceous plants.
The SWB [48] assumes that the soil has a specific water storage capacity in the root zone, called "available water-holding capacity" (AWC). The water held within the soil initially varies depending on the input due to precipitation and the output provoked by the actual evapotranspiration (ET A ), which withdraws water from the soil until running out or until its maximum potentiality, called potential evapotranspiration (ET P ). When the volume of water stored in the soil and the precipitation exceeds ET P , there is an excess of water that produces either surficial runoff, infiltration, or both. Such excess is known as effective precipitation (P e ). To compare different methods, the values of daily ET P were estimated by applying different equations: Thornthwaite [49], Hargreaves [50], and Blaney-Criddle [51]. The selected AWC values were 50, 60, 75, 100, and 125 mm. As the permeability of the bedrock in the study area is generally low, the runoff coefficients (RC) applied to P e were high: 50%, 75%, and 100%.
The computation of free-surface water evaporation was made using the climatic data recorded in the weather station placed in the center of the study area ( Figure 1a) by applying the Penman equation [52]. This evaporation model is one of the best for estimating evaporation from shallow wetlands (<2 m of water column) according to a meta-study by McMahon et al. [26], in which many previous studies were reviewed ([23,25,53-56], among others). However, the Penman equation does not work so well when applied to water bodies more than 2 m depth, as it does not take into account the heat transfer toward the water mass, which buffers the evaporation effect [57][58][59][60][61]. Vardavas and Fountoulakis [62] proposed an easy-to-apply modified version of the Penman model that considers the net energy exchange between the atmosphere and the water body, based on the monthly temperature variations measured in the wetland water and the monthly mean wetland stage. Since the water column in the Amarga wetland is several meters high, the Vardavas and Fountoulakis modification was also considered in the present work. Additionally, the Penman-Monteith equation [63] was included in this study to compare the results obtained from different methods.
The simulation of the daily changes on wetland storage during the observation period was made using their respective hypsometric and filling curves (Figure 3), rainfall recorded at the weather station (Figure 1a), and overland runoff and evaporation obtained from the different methods mentioned. The underground components of the wetland water budgets were inferred from the differences between the recorded and the computed limnimetric evolutions.

Results
During the four groundwater-level measurement campaigns, minor variations of the groundwater table were observed, with an average drop of 0.6 m from May 2015 to March 2016. In most wells, a general decreasing trend was observed (Figure 4), although in some cases, a groundwater level rise was registered between September 2015 and March 2016. That was mainly seen in wells near the Jarales wetland, particularly those emplaced at higher altitudes than the wetland itself ( Figure 4).  Figure 1 shows the location of the wells (in red). The evolution of the Amarga and Jarales stages is also shown.
The water table contour outlined in January 2016 is presented as an example in Figure 1a since the results obtained from the rest of the campaigns and the data retrieved from previous research [30] are relatively similar. According to the outline, the water table's morphology is somehow adapted to the land surface, typical of low permeability media. Thus, a potentiometric dome exists in the plateau's central area, creating a hydrogeological divide from which groundwater flows in a radial and centrifugal way toward the edges of the system.
The fit of each limnimetric simulation to the wetland stage record series was evaluated using the root-mean-square error (RMSE) and the correlation coefficient (R 2 ), whose values are shown in Figure 5. The RMSE measures the magnitude of the predictive error and evaluates the fit accuracy based on the differences between the computed values and the observed values. The R 2 tests the relationship between the two data series, indicating the covariance of both variables independent of the scale.
Regarding the runoff models, the worst results were obtained using the Thornthwaite equation into the SWB, which provides high RMSE values and low R 2 ( Figure 5). Hargreaves and Blaney-Criddle methods give more relevant results, although they are slightly better in the second case. The simulations derived from the CN equation are generally well correlated with the limnimetric record, and their RMSE values are low ( Figure 5). As for the evaporation models, the Penman equation generates better results in the Jarales wetland than those obtained from the Penman-Monteith method, while in the Amarga wetland, the best fittings were obtained with the Vardavas-Fountoulakis correction to the Penman model.  The JW3 simulation has the lowest RMSE value (0.19 m) and one of the best R 2 (0.957). However, in that simulation, the water stored within the soil never exceeds the AWC (nor in the rest of the simulations with identical RMSE and R 2 values in Figure 5). Consequently, no runoff generates during the heaviest rain events, contrary to in situ observations. If the field capacity is lowered to 75 mm and a RC of 50% is applied (JW4), the model simulates a certain runoff amount, but it does not replicate the flooding starting in November 2015 ( Figure 3). The two simulations that include the CN runoff (JW1 and JW2) better replicate the water level rises recorded, particularly in November 2015 ( Figure 6). The RMSE and R 2 values of JW1 are slightly better than those of JW2, partly because the last model does not predict the wetland's desiccation toward the end of the hydrological year 2014/15. Figure 7 shows the four selected simulations of the evolution of the Amarga wetland stage. Three of them were obtained applying the CN method with an I a of 15 mm and the evaporation according to Penman (AW1), Penman-Monteith (AW2), and Vardavas-Fountoulakis (AW3). In this case, the best fitting corresponds to an I a lower than the one in the Jarales basin (19 mm). That is in agreement with the lithological characteristics of the catchment area (Figure 1b,c), which is less permeable in the Amarga wetland case due to a more significant outcropping proportion of CSC rocks that confer the soils a clayey texture. The other one (AW4) was computed based on the SWB, with the ETP after Blaney-Criddle (AWC: 50 mm, RC: 75%) and the evaporation of Vardavas and Fountoulakis. AW4 is one of the simulations with the lowest error ( Figure 5), as it only has slight differences with the actual values for most of the recording period, except for part of 2016 ( Figure 7). However, the computed water table rises do not resemble the observed ones, and some of the actual ascents were not replicated. The latter occurs because the runoff generating such water level changes is caused by rainfall occurring when the soil was not saturated. Contrary to the SWB, the CN method (used in AW1, AW2, and AW3) predicts the runoff generation in those circumstances, as it does not depend on the previous state of the AWC but on the total amount of rain. Based on the previous analysis, JW1 and AW3 simulations were closer to the actual records. Thus, each wetland's water budget (Tables 2 and 3) was calculated on a monthly and annual scale, using direct rainfall precipitation inputs and the evaporation and runoff simulation used in JW1 and AW3. The difference between the actual and simulated water storage volume was attributed to the net groundwater-wetland exchange (∆G), meaning the difference between groundwater inputs and outputs due to infiltration in the wetland bed. Thus, positive values of ∆G stand for a net underground input. During 2014/15, evaporation was the most significant component in the Jarales wetland budget (34.9 dam 3 ; Table 2). Direct precipitation inputs were 13.5 dam 3 , while the runoff contribution was 1 order of magnitude lower (1.     (Table 3), while precipitation inputs are 20.3 dam 3 and 15.0 dam 3 , respectively. Runoff contributions vary notably, from 13.1 dam 3 in 2015/16 to 2.3 dam 3 in 2016/17. This last value must be an underestimate, as the simulated level rise in November 2016 is much lower than the actual one ( Figure 7). That underestimation could relate to the fact that the intensity of rainfall was lower in this period than in the previous years; thus, the CN method did not predict much runoff. The total simulated variation of stored water (∆V) between 2015/16 and 2016/17 is −90.2 dam 3 (−40.8 and −49.4 dam 3 , respectively), which would reduce the stored volume at the end of the simulation period (September 2017) to 37.7 dam 3 (Table 2), in contrast to the actual storage (91.7 dam 3 ).

Choosing the Runoff and Evaporation Equations for the Limnimetric Simulation
The correlation degree observed between the actual limnimetric record and the simulations made with the CN method is generally good, and the magnitude of the errors is low ( Figure 5). That could be due to the low permeability of most of the rocks outcropping in the wetland watersheds. Thus, water does not infiltrate rapidly into the soil, resulting in runoff generation during heavy rain events even if it is not completely saturated. Such a process fits with the theoretical basis of the CN. The dry hydroclimatic conditions registered during the study period would have favored that the water held in the soil rarely exceeded the AWC so that the SWB was not able to predict the runoff generation in certain moments. The best simulations obtained applying the SWB considered low values of AWC (50-75 mm, Figure 5) that were coherent not with the characteristics of the outcropping rocks but with the limited thickness of the soil and its limited radical depth. Both factors combined with the semi-arid regime of the study area, and the climatic conditions (dry years) would have enhanced the weight of the ET P , which impedes the SWB to be suitable to the studied reality. Based on all that, the CN method is the most adequate for estimating runoff in the explained context.
Regarding the estimation of free-water surface evaporation, the simulations integrating the Penman model fit slightly better in the Jarales wetland ( Figure 5). Additionally, the JW1 simulation predicts the wetland's desiccation in the summer of 2015 using this equation, while JW2 does not replicate this issue using the Penman-Monteith evaporation ( Figure 6). The good applicability of Penman equation is in line with other research comparing methods for the computation of free-water surface evaporation [24,26]. In shallow lakes, like the Jarales wetland, the Penman equation can be used to estimate lake evaporation, as the advected energy and changes in seasonal stored energy can be ignored. However, the Penman-Monteith method was designed to calculate reference-crop evaporation, and it includes simplifications of the atmospheric pressure and the latent heat of vaporization [63] that affect open-water evaporation.
On the contrary, the RMSE values derived from the simulations using Penman-Monteith in the Amarga wetland are generally lower than those obtained with the Penman model, although the best adjustment was reached using Penman plus Vardavas and Fountoulakis correction ( Figure 5). The Penman equation's assumptions are not entirely applicable to deep water bodies [26,60], such as the Amarga wetland. On one hand, the runoff inputs to the Amarga wetland can create rises of tens of centimeters in the water column. Such inflows may have a different temperature than the average water temperature in the water body, so advected energy cannot be neglected. On the other hand, the water body can store and transport heat, so heating and cooling must also be considered as energy fluxes [64]. As the depth of the mixing varies in time and space, it is hard to estimate changes at a short time step, although a number of methods offer satisfactory approaches [57,59,61]. The Penman modification by Vardavas and Fountoulakis is one of them [62] and is based on monthly surface water temperatures, which are easy to acquire. The method reached satisfactory results in Australia, including case studies with Mediterranean and semi-arid climates. In the Amarga wetland, with a similar climatic setting, the limnimetric simulation is also improved when applying the modification. Thus, it is highly recommended to check its applicability in other wetlands in similar contexts.
Previous work carried out in the study area [12,29,32] performed limnimetric simulations and wetland water budgets using the SWB. They reached satisfactory results using AWC values between 103 and 260 mm, which are high considering the limited soil thickness. Said works considered average monthly data and/or wet years from the standpoint of precipitation amount. The use of those values for the period 2014/15-2016/17, integrated with dry and medium years, would not have resulted in daily-stepped runoff. In wet years, the rain widely exceeds the AWC, even if it has relatively high values. Additionally, the monthly-stepped simulations use accumulated values of precipitation and ET P , so they are less precise in predicting the time distribution of the soil saturation status. Therefore, the hydroclimatic conditions and the control and simulation steps are essential for selecting both the equations and the fitting parameters. Thus, generating long daily data series of wetland stage variation would be advisable for obtaining a model able to predict the limnimetric evolution of the wetlands under any climatic conditions.

Integrating the Results into a Hydrogeological Conceptual Model
During most of the observation period, the Jarales wetland's hydrological dynamic was mainly influenced by the surficial components of its water budget, particularly evaporation and direct precipitation to a lesser extent ( Table 2). The groundwater-wetland interaction only represents a small volume of the whole budget but explains the difference between the actual limnimetric record and the computed series ( Figure 6). These results agree with the functioning proposals based on average data generated by previous research [29,30]. According to the monthly water budget of the Jarales wetland during the first six months of control (May to October 2014), the loss of water storage is greater than the evaporation outputs ( Figure 6), which suggests that there is net wetland water infiltration (∆G < 0, Table 2). On the contrary, throughout most of the year 2014/15, net groundwater inputs were estimated. This input would explain the nearly steady state of the wetland water level from December 2014 to May 2015, which does not fit the downward trends simulated, particularly since March ( Figure 6).
The Jarales wetland is placed in an intermediate position with respect to the water table contour outline (Figure 1a), between the potentiometric dome in the central part of the study site and the discharge areas situated at lower altitudes. Thus, Jarales would be a flow-through wetland concerning groundwater [5,7,65], meaning it would receive groundwater inputs in some locations (SW), while infiltration losses would occur in others (NE). The water-table changes nearby the wetland can modify the hydraulic gradient around it, making it possible to explain variations in the net groundwater-wetland exchange direction [6,66]. The high water level in the Jarales wetland at the beginning of the control period ( Figure 6) resulted from essential inputs during the previous year. This would have caused a hydraulic gradient rise, favoring water infiltration toward areas with a lower table. That explains the negative values of ∆G from May to October 2014 and the difference between the actual and the simulated limnimetric evolutions. In 2014/15, the evaporation rate extensively exceeded precipitation and runoff inputs, causing a pronounced fall of the wetland stage. The low permeability of the terrain would have smoothed the water table's descent, and therefore a higher gradient would have been established between recharge areas (potentiometric dome, Figure 1) and the wetland. The result is a positive volume of net groundwater-wetland exchange (groundwater input) from November 2014 to May 2015 (Table 2). Eventually, the water table would drop downflow the wetland during summer 2015 (Figure 4), favoring water infiltration and leading to wetland desiccation. The subsequent water table rise upflow the wetland in the first half of the hydrological year 2015/16 may have helped the intermittent swamp of the area, although the general descending trend observed in the downflow wells would have favored infiltration and impeded long-lasting flooding. In any case, this whole interpretation is made based on a short control period, and it will be necessary to confirm it with a more extended observation period that includes both dry and wet years and more detailed groundwater table monitoring.
The permanent hydroperiod in the Amarga wetland is only possible if a considerable net groundwater input exists, as deduced from its ∆G values (Table 3). The estimated volume of that contribution ranges between 20.9 dam 3 /year (2015/16) and 33.0 dam 3 /year (2016/17), in agreement with previous research. Aljibe Consultores [31] calculated 28.94 dam 3 of average groundwater contribution, and Moya [67] estimated 33.6 dam 3 for the hydrological year 1983/84 (average from the rainfall amount standpoint). Nonetheless, an overestimation of the groundwater input in 2016/17 is possible due to the undervaluation of runoff mentioned before. The Amarga wetland's location toward the potentiometric dome border, together with its deep basin, favors the water table to be always above the wetland bed, as observed in the nearby wells ( Figure 4). As a result, groundwater input is permanent and exceeds the water infiltration outputs. Though the gross groundwater inflow to the Amarga wetland was not computed, it could be higher than the net value, as the presence of discharge pointing toward lower altitudes would enable water infiltration downflow. Additionally, locals declare to have observed water infiltration through karst swallow holes formed in the wetland side, and the hydraulic connection between nearby endorheic areas and one of the springs was proven by dye tracer tests [36].
Thus, the differences in the Jarales and Amarga wetlands' hydric functioning are due to hydrological aspects related to their watershed characteristics and hydrogeological ones linked to wetland-groundwater interactions. A hydrogeological conceptual model ( Figure 8) can explain the wetlands' functioning based on a previous model proposed by Andreo et al. [65] and proved after by Gil-Márquez et al. [68], which describes the general hydrogeological functioning of the CSC. The higher emplacement of the Jarales wetland and its flatter bed impede the permanent wetland-groundwater interaction. When the water table reaches the wetland bed, either net groundwater input or infiltration outputs may occur, depending on the wetland's water level. While groundwater inflows would be related to local flows, the water infiltration occurring in the wetland bed would be incorporated into the regional groundwater flow net. On the other hand, the Amarga wetland position toward the potentiometric dome border would permit permanent groundwater inputs from nearby areas but would also recharge the region at the highest part of the plateau. The infiltration that may occur in the Amarga wetland would incorporate into the regional flow system that would eventually be drained through the brine springs existing at the plateau's border [65,68]. Previous research has collected hydrobiogeochemical evidence that agrees with the proposed conceptual model (Figure 8). The Clconcentration and the enrichment in δ 18 O of both wetlands' waters are coherent with the hydrochemical and isotopic signature of groundwater from nearby wells after undergoing evaporation [35]. Thus, wetlands are likely receiving groundwater inflows. On the other hand, the analysis and joint interpretation of 3 H, 3 He, 4 He, and 14 C in groundwater samples indicate that several flow components with different ages converge at the brine springs located in the system's border [68]. Finally, the determination of δ 15 N and δ 18 O of dissolved NO 3− and δ 13 C DIC in groundwater samples revealed that nitrate removal processes occur at regional scale, linked to organic carbon oxidation taking place in the wetlands [69]. Therefore, water that is infiltrated in the wetlands is incorporated into flowpaths of different lengths until reaching the discharge areas. Not only does the result derived from the present work fit with the above statement, but it also permits an explanation of the hydric variability of the wetlands.
Several factors may condition the wetland-groundwater connection, including the basin geometry [75], the lithology [80], the wetland position with respect to the groundwater table [5,18,74], and its location within the regional groundwater flow networks [73,81]. In this last sense, Smith and Townley [82] simulate 39 flow regime geometries and found that flow-through wetlands are dominant within regional flow systems, particularly when shallow lakes intersect high-gradient regional flow systems, typical of low-medium permeability media. Neff et al. [83] proved that the presence of depressional wetlands within a broader landscape modifies the regional groundwater flow regime, causing the appearance of flow-through or discharge wetlands far from the regional discharge area. The presence of numerous depressional wetlands within the CSC and the relatively high groundwater gradient can favor the wetland-groundwater interaction, as is suggested by this research.
Previous research proved the wetland-groundwater interaction in the Amarga wetland [12,31], but the present work also reveals the groundwater connection with the Jarales wetland. In both cases, the protected areas are buffers that do not fit with the actual watersheds (Figure 1), so the management measures are not able to act on processes that can affect the wetland through surficial runoff, particularly the carryover of agricultural pollutants and sediments that can potentially lead to the wetlands' clogging. Moreover, as both wetlands receive regional groundwater inflows in different proportions, it is crucial to advance on the definition of the wetland catchment area, as human activities can affect their hydrological functioning beyond their surficial watershed. Future management, conservation, and restoration policies must take that into account.

Conclusions
This work analyzed the relation between groundwater and two Mediterranean wetlands with evaporite karst origin through the interpretation of water table contour outline and the simulation of wetland level fluctuations. The latter has allowed quantification of the wetland water budget components and estimation of the existence and weight of the wetland-groundwater interaction. A hydrological functioning model was proposed in which the relation of the wetlands (mainly Jarales) with the regional groundwater flow network much condition its hydric behaviors. Thus, it is necessary to control the hydrological factors affecting its water budget, including the groundwater surface, so that the hydrogeological context can be fully understood.
The use of different methods for the estimation of runoff and open-water evaporation permitted a contrast of their applicability to the simulation of water level variations of wetland in a semi-arid context. Results evidence that the Penman model is better for evaporation computation, although the correction proposed by Vardavas and Fountoulakis improves it in deep wetlands. The curve number method provided better results than the soil water balance method in dry years, in which the water stored in the soil rarely exceeds the field capacity, but runoff generates nonetheless when rainfall overcomes the soil's capacity (predominately clayey) to infiltrate water.
The difference between the calibration parameters included in the simulation presented in this work and previous studies evidence the need for performing a continuous limnimetric record in wetlands for long and uninterrupted periods. The obtained series can be the base for developing and calibrating hydric functioning models that apply to any hydroclimatic conditions and that help reach a more accurate knowledge of the wetlandgroundwater relationship. That kind of model will be essential for the management, protection, and restoration of wetlands, particularly in areas that are particularly sensitive to climate change, such as the Mediterranean area. Only through an understanding of the origin and fate of water in wetlands is it possible to assess how future scenarios would affect hydric functioning and to propose adequate conservation measurements.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ongoing research using it.