Received 16 March 2015; accepted 14 July 2015; published 17 July 2015

ABSTRACT

This research is concerned with new developments and practical applications of a physically- based numerical model that incorporates new approaches for a finite elements solution to the steady/transient problems of the joint ground/surface water flows. Python scripts are implemented in Geographic Information System (GIS) to store, represent and take decisions on the simulated conditions related to the water resources management at the scale of the watershed. The proposed surface-subsurface model considers surface and groundwater interactions to be 2-D horizontally distributed and depth-averaged through a diffusive wave approach for surface flood routing. Infiltration rates, overland flows and evapotranspiration processes are considered by a diffuse discharge from surface water, non-saturated subsoil and groundwater table. Recent developments also allow for the management of surface water flow control through the capacity of diversion on river beds, spillways and outflow operations of floodgates in weirs and dams of reservoirs. Practical application regards the actual hydrology of the Mero River watershed, with two important water bodies mainly concerned with the water resources management at the Cecebre Reservoir and the present flooding of a deep coal mining excavation. The MELEF model (Modèle d’ÉLÉments Fluides, in French) was adapted and calibrated during a period of five years (2008/ 2012) with the help of hydrological parameters, registered flow rates, water levels and registered precipitation, water uses and water management operations in surface and groundwater bodies. The results predict the likely evolution of the Cecebre Reservoir, the flow rates in rivers, the flooding of the Meirama open pit and the local water balances for different hydrological components.

Keywords:

Integrated Surface/Subsurface Flows, Numerical Modelling, Finite Elements, Watershed Hydrology, Geographic Information System, Water Resources Management

1. Introduction

At the present moment there is an increasing need for integrated surface and ground water numerical modelling. The philosophy and role of hydrological models in water resources have been widely described (Abbot, M.B. and Refsgaard, J.C. [1] , Burnash, R.J.C. et al. [2] ).

Several process-based numerical models for the simulation of coupled surface and subsurface flow, at the scale of a particular region, have emerged in recent years. Examples include Querner, E.P. [3] , who mainly has developed the surface and groundwater flow model MOGROW; Restrepo, J.I. et al. [4] , who developed a new module for the MODFLOW groundwater model in order to be applied to wetland regions; Panday, S. and Huyakorn, P.S. [5] , with a fully coupled physically-based spatially distributed model for surface/subsurface flow; Erduran, K.S. et al. [6] , whose integrated model focuses on the coupling of shallow surface and saturated groundwater flows; Gunduz, O. and Aral, M.M. [7] , who focused on the interactions between rivers and aquifers; Langevin, C. et al. [8] , whose integrated surface water-groundwater flow model is applied to a coastal wetland and estuary; Kollet, S.J. and Maxwell, R.M. [9] , whose coupling approach is based on the continuity of the pressure head and fluxes at the ground surface; Jones, J.P. et al. [10] , whose coupling of Richards equation and various approximations of the Saint-Venant equation relies on the hypothesis of a first order diffusion of water between the surface and the subsurface and an interface layer with a finite thickness; Shan, H. et al. [11] , who integrated flow and transport in surface water, groundwater, and overland regimes; Spanoudaki, K. et al. [12] , whose 3-D integrated surface water-groundwater model behavior was verified; Weill, S. et al. [13] , who developed a model based on the generalization of Darcy-Richards equation that allowed implicit coupling between surface and subsurface; Sparks, T. et al. [14] , based on recent extensions, developments and verifications of DIVAST model, in order to more accurately describe the interactions between surface and groundwater flows; and García-Rábade, H. et al. [15] , who developed a new coupled finite element model for the joint resolution of the shallow water and the groundwater flow equations in order to solve problems related to surface and subsurface water flows.

Recent numerical methodologies were also developed on combined watershed and ground-water applications to the whole of the water resources of a particular river basin (Ross, M. et al. [16] ; Sophocleous, M. and Perkins, S.P. [17] ). MIKE SHE and MIKE BASIN are two examples of numerical and physically based modeling systems developed by DHI [18] for describing the major flow processes of the entire land phase of the hydrological cycle which integrates the Saint-Venant surface equation to a vertical 1-D Richards equation for unsaturated flow and a 3-D finite elements solver for saturated flow (Graham, D.N. and Butts, M.B. [19] ); Camporese, M. et al. [20] , whose catchment’s hydrology model couples a finite element solver for the Richards equation describing variably saturated porous media, and a finite difference solver for the diffusive wave equation describing surface flow.

Standard features of detailed process-based numerical models for water resources management at the scales of the watershed could include: a large variety of numerical conditions and sub-models for the different hydrological components, water uses and managements. These characteristics will not be described in any case here; this research will focus instead on some new features and developments, summarized below, of an integrated surface-subsurface flow numerical model for suitable applications to practical water resources management in regulated watersheds.

In this field, a new methodology of finite element modeling has been developed (MELEF model) which considers novel modeling features of the joint surface and ground water regional flows of continental freshwater and coastal salt water intrusion from the sea by means of an immiscible fresh/salt water interface reliable for applications on different kinds of watersheds, water resources management and environmental problems (Padilla, F. and Cruz-Sanjulián, J. [21] ; Padilla, F. et al. [22] ).

The model considers surface and groundwater interactions to be depth-averaged through a novel interpretation of a diffusive wave approach for surface flood routing. Infiltration rates as well as overland flows generation processes are assessed by the concerned sub-models. Evaporation and evapotranspiration processes are considered by a diffuse discharge sub-model from surface water, non-saturated subsoil and groundwater table (Hernández, J.H. et al. [23] ).

Therefore, this research deals with the numerical developments being required to consider new practical applications of this finite elements solution, in particular those concerned with water uses and management of surface water, through the water flow control of the capacity of diversion works on river beds, spillways and outflow operations of floodgates in weirs and dams of reservoirs. Furthermore, Python scripts were developed and implemented in a Geographic Information System (ArcGIS 9.3) in order to manage the new capabilities by polygons (shapefiles). Finally, it creates a new shapefile of nodes with all the information, by spatial join geoprocessing, to generate the model simulation input files.

The practical application regards the actual hydrology of the Mero River watershed where there are two important water bodies. These are concerned with the water resources management at the Cecebre Reservoir and the present flooding of the Meirama open pit, a deep coal mining excavation, being filled with freshwater coming from the upper Meirama sub-basin, all of it in the context of the water resources uses, management and freshwater hydrological aspects of groundwater and surface waters at the Mero River catchment (~250 Km²), Coruña, Spain. The model behavior is checked, during a simulation period of almost five years (2008/2012), featured by a complex geology and the particular geomorphology and hydrology of the opencast Meirama mine, being presently restored as a deep and regulated lake, as well as the Cecebre Reservoir, for the historically registered precipitation, water uses and water management operations in surface and groundwater bodies.

2. Model Description

2.1. Surface-Groundwater Interactions

The MELEF model for continental and coastal catchment applications couples subsurface and surface regional hydrology, by a diffusive wave approach, with a joint finite elements solver for the saturated porous media flows of fresh and saltwater through an immiscible sharp interface (Huyakorn, P.S. et al. [24] ; Padilla, F. et al. [22] ). Moreover, the application of the model is extended to fractured media and hydraulic anisotropy.

In this sense, this numerical approach (Figure 1) was mainly developed in order to analyze the surface and groundwater behavior at a watershed scale, and it can also consider other water components such as the evapotranspiration process, that is a diffusive discharge from surface water and soils within the unsaturated zone by a

root water uptake sub-model, as well as the overland flow by a rainfall-runoff sub-model based on an exponential method for assessing the infiltration rates (Hernández, J.H. et al. [23] ).

2.2. New Numerical Conditions

2.2.1. Diversion Capacity

Diversion works may be considered for bypassing the flow from a river bed toward another place, by transferring it mainly through an impervious channel, pipeline or tunnel. The place where the water is transferred could be located in the same catchment area or in a different one. Ordinarily, it is not always possible to bypass all flows; this would depend on the required capacity of the diversion works for a particular weir.

In order to take into account the diversion capacity of a weir built on a river bed, the following equation has been implemented, defining the velocity of diversion V, in terms of the flow rate being diverted per unit area of flooded surface.

(1)

where z_s is the thickness of the stream or overland flow (in meters), and V₀ is a parameter defining a particular velocity of the surface water diversion (in terms of flow rate diverted per unit area of flooded surface).

Therefore, it is being considered that the diversion works that control the regulation capacity of a particular weir from a river bed could be reproduced by properly defining the parameter V₀, that is, the particular velocity of the surface water diversion that would make sense for a surface water depth of 10 meters in the stream, and that in terms of the flow rate being diverted per unit area of flooded surface. Then, Equation (1) has a double purpose: firstly, the velocity V₀ defines the flow rate that is diverted per unit area of flooded surface, mimicking an efficiency curve for the diversion rate that is a function of the surface water depth; and secondly, at the same time, it provides a smooth increase of overflows in order to avoid numerical fluctuations.

2.2.2. Pipe Spillways

Nowadays, the major part of water transmission is carried out through pipes, which has a considerable increase in the number of new types. Pipe spillways are often used with dams in order to regulate part of the water flow at the outlet of a surface reservoir. They have relatively low capacities, when compared to floodgates, but they carefully regulate the water stored behind the dam and reduce the peak rates of outflow.

The inflow capacity of a pipe spillway at full flow can be hydraulically estimated by the flow rating curve Q as follows:

(2)

where: D (m) is the diameter of the pipe; g (m/s²) is the gravity acceleration; H (m) is the total head; and ±ΔQ (m³/s) could optionally be the unexpected flow rate losses or leakages of water through the dam, which will also probably be dependent on the hydraulic head above the spillway mouth (Figure 2).

In this context, the above flow rating curve of the pipe spillway needs to be adjusted to the practical flow rating curve of the pipe itself, according to the curves provided by the operator, for the practical diameter of the

pipe (D). The latter would take into account the different openings (d) of the valve installed in the pipe, and a way of defining the discharge coefficient C_d (where D = d∙C_d), that is, the hydraulic losses of the pipe system due, for instance, to the installed valve gate itself (sliding gate, butterfly, etc.), as well as to the pipe inlet, the turbulences and other frictions in the pipe itself (inlet, pipe bends, pipe lengths, etc.). Therefore, there is a need to calibrate the above flow rating curve to the real one, for the different openings of the pipe, through the appropriate choice of the corresponding and calibrated openings (D) of the actual pipe.

Optionally, another parameter that could be calibrated is the unexpected flow losses or leakages of water through the dam (±ΔQ), whose dependence on the hydraulic head would also be subject to calibration. In order to supply right away an answer to this issue, the following water loss function is proposed:

(3)

where: α is a coefficient of adjustment; and ±ΔQ^* (m³/s) could optionally be the maximum unexpected flow rate losses or leakages of water through the actual dam (+) or, instead, through the modeled dam (−), for an approximately maximum total head H (m).

2.2.3. Floodgates

Dam floodgates at the outlet of surface reservoirs are devices that usually have relatively high capacities to regulate the water outflow rates. The types of floodgates and how they operate are nowadays quite different.

To this concern, the crest floodgate, or overflow spillway, is the most appropriate to be managed from the point of view of a numerical model. The main reason is because they operate as weirs and one of the controlling variables will be the gate crest itself. This way it could be possible to assess the water depth at the gate anytime, provided that the surface water levels being reached nearby in the reservoir are known. In order to properly operate with these variables, and to correctly define the outflow rating curves through the dam gates, a good numerical approximation of the gate shape would also be required, as well as an accurate assessment of the surface water levels and flows around the operating gates. Therefore, a fine surveying of the floodgates and their transient evolution is of main concern. Then, if we consider an overflow spillway with a type of gate similar to the vertical lift crest floodgate at the outlet dam of a reservoir, we would not require any special treatment, as far as the concerned shape is provided, the heights are continuously defined, and the free water levels are known.

However, a vertical lift crest floodgate is not usual as a reservoir outlet, mainly because it requires a high superstructure. It is much more common to install a sluice gate, which will work submerged with its lower edge more or less close to the sill of the dam. This type of floodgate is also called a vertical lift gate, its installation is simpler and its maintenance is easier. The latter is a submerged device that reduces the water depth, by decreasing the free water level, and it may be opened just by rising, leaving an opening (or orifice) through which the water flows, underneath an usually higher water level.

If we try to find an equivalence between the opening h of a vertical lift gate and the water depth h' at a suppressed rectangular sharp-crested weir, we could easily find the required height at the vertical lift crest gate, and that in order to achieve an equivalent water outflow rate (Mays, L.W. [25] ). The appropriate hydraulic expression could be the following:

(4)

where, H is the water head above the sill of the vertical lift gate, which can be calculated from the known free water level; h is the opening of the vertical lift gate; and h' is the searched water depth over the vertical lift crest gate. Therefore, we can continuously assess the height of the corresponding vertical lift crest gate, provided that we know the free water level at the floodgate anytime: height = free water level − h' (Figure 2).

This means that it would also be possible to properly assess the water outflow rate at a reservoir outlet which has a dam provided with vertical lift floodgates, whether the openings of the gates are known and we are able to assess the free water levels at the gates anytime. This would also imply a continuously moving height at the gates of the dam. Therefore, the floodgates being actually supported in the MELEF model are of those two types: lift crest gates (or overflow spillways) and vertical lift gates (or sluice gates). In this sense, it can be concluded that an appropriate numerical modeling of the gates position and their shape is also of main concern, but nevertheless it could be a task of great practical difficulty.

2.3. Boundary Conditions and Resolution

MELEF, the current tool used in water resource modeling, is a two dimensional finite elements model for regional surface and groundwater flows through drainage basins, developed for a temporal implicit (Eulerian) centered (Crank-Nicholson) and spatially centered (Galerkin) numerical approach.

The inputs of freshwater are exclusively due to rainfall, however pumping and injection wells of prescribed flow rates have been also implemented, as well as irrigation and artificial recharges, which could be also prescribed inside the considered domain. Therefore, the boundary conditions of the modeled region, at the scale of the watershed, are mainly those of the outflow or discharge conditions. To this respect, prescribed water heads, outflow face and open boundary conditions are commonly used at the discharge zones of the drainage basin (Padilla, F. and Cruz-Sanjulián, J. [21] ).

Triangular elements of three nodes are used for the analytic integration of the corresponding numerical formulation for steady and transient conditions. The preconditioned iterative algorithm GMRES (Saad, Y. and Schultz, M.H. [26] ) provides the solution to the system of equations by means of a reduction in the computer memory requirements and allowing for a simple processing of the numerical mesh. There is no difference in time-scales between the surface and ground water flows because they are solved together in a whole lumped and depth-averaged way. This is a clear advantage because the lumping properties allow for good numerical results in practical applications with affordable meshes as well as time scales. Hydraulic heads and gradients are the solutions. Then the model is depth-averaged solved, but the results are interpreted and approached in a multilayered way, giving solutions for water flows and mean horizontal velocities above the soil surface, runoff and overland flows, as well as under the soil surface, saturated subsurface and groundwater flows.

As part of the present surface/subsurface flow numerical approach, the model incorporates capabilities to assess the drainage layout of surface runoff and the freshwater ground levels, overland and subsurface flows, with punctual and diffuse surface and ground water recharges and discharges, thickness and velocities of surface and ground waters, as well as several types of rivers diversions, water balances and flooding of water bodies of high depth level.

2.4. Simulated Conditions and GIS

First of all, a Geographic Information System (GIS) is implemented by a Digital Elevation Model (DEM, Maidment, D. and Djokic, D. [27] ) which is required as a pre-processor to define the watershed study area. Afterwards both ArcHydro Tools and Agree-DEM surface reconditioning system method (Hellweger, F.L. [28] ) are carried out. Then, flow accumulations and flow directions are considered for stream definition and delimitation of surface-groundwater interaction buffer zones. Thereafter the watershed study area is implemented to create a triangular finite elements mesh with major density nodes at buffer zones. Buffer zones, or polygons, are generated by a python script that uses the streams layer, for an initial width of these streams stored in the attribute table and a seed value that controls the incremental distance between buffers by an exponential relation, according to which the script creates separated buffers and their vertex are redistributed. The seed value is used to define the faster or slower evolution from the zones of higher density nodes to the zones of lower density nodes. Thereafter a polygon layer of the watershed study area, with buffer zones, is implemented to create a triangular finite elements mesh with zones of different density nodes controlled by the vertex distance of polygons.

Then, MELEF code simulated conditions and the discrete model mesh are managed with Model Builder and Python scripts in ArcGIS 9.3.1. The Mesh definition and the nodal properties are controlled by shapefiles (points, polygons and polylines).

Finally the shapefiles with the triangular elements and nodes of the mesh, pumping wells, cross sections discharges, surface water diversions zones, boundary conditions and rain gauged zones are generated. A database, with measured or acquired data, is related with the table of the shapefile attributes through a common field to generate spatially and temporally transient conditions.

3. Model Application

Since 1980 Meirama coal mine had been exploited by Lignitos de Meirama, S.A. (LIMEISA) as part of an important activity that finished in March 2008. As part of the environmental plan of closure the opencast pit hole was to be flooded to finally create a large lake of ~1.86 × 10⁶ m² up to a maximum depth level of some 200 m and a capacity of ~146 × 10⁶ m³. In this context, the forthcoming evolution is analyzed as a function of the actual strategy of flooding. The presently forming pit lake of Meirama is to lie on the drainage basin of the Barcés River, a tributary in the Mero River catchment (247 km²) leading the flow towards the Cecebre Reservoir, north- western Spain, 20 km from the city of La Coruña (Figure 3).

In order to improve the calibration of the stream flow model, continuous water depth measurements were carried out with data logger transducers on 8 of the main streams (Schlumberger Mini and CTD DIVER) at the Mero watershed (labeled with S1 to S9 in the Figure 3) from June 2011 until late September 2012. Water depth measurements were corrected by barometric pressure. Water discharge measurements by velocity area and dilution methods, for different regimes, were made on these main streams to obtain rating curves and to transform water depth measurements into water discharges.

With respect to the geological aspects, there are three main regions in the Mero River catchment, one in the NW featured by a fractured and altered granite massif being partially kaolinized, a second one in the S-NE characterized by a schist substratum partially covered by alluvial Quaternary sediments, and a third one marked by a tectonically complex sedimentary basin with Tertiary materials and mining refills which is located between the two former zones of granite and schist and where most of the mining works of the lignite exploitation took place (Figure 4). For these materials, the main hydrologic parameters (hydraulic conductivity, drainage porosity, infiltration rate, soil thickness, capillary fringe, …) have been adjusted during the simulation period (going from January 2008 to September 2012) mainly based, on one hand, on the field measured parameters and on the tectonics and geological features of the region, as well as, on the other hand, on the measurements of the surface water flows, the free surface levels at the Cecebre Reservoir and at the pit lake of Meirama, as well as on the known groundwater phreatic levels and the water balances inside an area delimited around the pit (Figure 5).

The climate in this area is typically Atlantic with an annual rainfall of about 1300 mm and a mean evapotranspiration of 600 - 800 mm, nevertheless the rainfall in the catchment presents a positive linear gradient conditioned by altitude. To this respect, the catchment might be split up into three different zones by rain gauging station altitudes (Figure 5), where historical daily precipitation rates, temperatures and insolation were considered during almost five years 2008-2012. Thereafter, daily historical precipitation was used and potential evapotranspiration rates were estimated with climate variables over the drainage basin.

During the last exploitation period (until March 2008) some twenty pumping wells were placed at Meirama pit mine contours so as to drawdown the water table. The water pumped from the bottom of the pit was then led to a treatment plant, being afterward drained by means of two perimeter channels towards the Barcés River, tributary of the Mero River (Figure 5). In order to prevent or to allow the surface water getting into the pit during the flooding period (after March 2008), the inlet streams were gradually diverted through adjustments of flow rates in weirs (Figure 5). Other water uses within the Mero River catchment area of main interest are, for instance, the registries of the regulation operations provided by the water supply company of La Coruña, EMALCSA (mainly the flow rating curves and the openings of spillways and floodgates) for the surface water outflows at the dam outlet of the Cecebre Reservoir (Figure 5), as well as others surface water diversions for water supply for industries and neighboring municipalities.

With respect to this, the present work applies the MELEF model to the fresh water resources of the whole Mero River catchment, which includes, among others, the regulation operations of a free surface reservoir and the exploitation, closure and hydraulic restoration of the opencast mine, with the aid of the historically registered hydrological parameters and data. Thus, the simulation strategy considers a period of almost 5 years, going from 1-January-2008 to 30-September-2012. The simulated conditions consider measured daily rates of precipitation, and the main management of groundwater and surface water at the Cecebre Reservoir and the Meirama coal exploitation and flooding of the mine. Once the flooding was initiated on 18-March-2008, most of the existing wells stopped pumping (Figure 5) and the flooding of the mine was started with the groundwater flow. After 3-October-2008 the mine is being flooded with water diverted from some of its neighboring streams (Figure 5). After 15-June-2011, the surface flow of the main stream (Pereira creek), which up to that date had been diverted to the open pit, was thereafter split up, by means of a weir constituted by a spillway (with a fixed opening) and an overflow system (with fixed position), into both the flooding of the Meirama lake and an attempt to reach the minimum environmental flow which was to be supplied downstream through a tunnel to the Barcés River. Thereon, the mean monthly value of the ecological discharge was estimated as about 0.2 m³/s.

In this context, the calibration of the hydrological parameters has been a difficult task, mainly because of the lack of homogeneous groundwater levels in the most fractured rocks materials in the upper basin of the Meirama Lake, where there are quite vertical hydraulic gradients, as well as isolated and perched aquifers. These are not the most appropriate hydrological conditions for the present horizontal, vertically averaged, porous media modeling approach.

Otherwise, it is quite a challenge to reproduce the outflow regulation operations of Cecebre Reservoir dam at its outlet, through the particular and changing openings of its five floodgates and its three submerged pipe spillways (one near the middle, and two deep). In respect of this, the flow rating curves of the pipe spillways were calibrated, on one hand, to the actual flow rating curves of the middle pipe, and on the other hand, to the two deeper pipes being considered together. Therefore, the diameters of the pipes were calibrated for the different openings of the valves installed at the pipes and the corresponding discharge coefficients, and this according to the flow rating curves provided by the operator, (Figure 5 and Figure 6).

With respect to the five sluice gates (vertical lift floodgates), they were placed together in the numerical model, as if they were just one bigger floodgate (with lumped position, shape and openings), during the almost five years of operation, and that in order to try to reproduce their hydraulic behavior from the point of view of the surface discharges as well as the surface and groundwater levels in the reservoir and its surrounding area (Fig- ure 6).

The operating policy of this quite small Cecebre Reservoir requires that there should be a safeguard to ensure the filling up to about 70% of its total capacity, i.e. 14.5 × 10⁶ m³, from November to March, both of them inclusive to prevent eventual floods. From March on, the management is focused on filling the reservoir as much as possible, about 20.5 × 10⁶ m³, in order to meet the demands that take place during the summer.

In order to properly achieve this whole issue, the initial groundwater and surface water levels of the simulation period in the whole watershed need to be a steady-state solution, and this is of main concern, because those simulated water levels required to be quite similar to the actual ones. This is the most important reason why the simulation period starts on 1-January-2008, which was a short period of actual severe drought in the region, where the groundwater levels and surface water discharges were very low (with a quite dried Meirama open pit, still exploiting the coal, and the Cecebre Reservoir having very low free water levels), quite close to a steady- state numerical solution for those particular hydrological conditions.

Then, the aim is to reproduce the actual regional hydrology of the Mero River catchment, during these almost five years of water management (2008-2012), focusing on the actual evolution of the present flooding of the Meirama open pit, and on the estimated surface water discharges at the main streams, as well as on the evolution of the free water levels at the Cecebre Reservoir.

4. Results and Discussion

The results of the application of the model to the whole of the water resources at the Mero River catchment, during a total period of almost five years, that includes the recent hydrological history of the destination of surface and ground waters, related to the Cecebre Reservoir and the flooding of the Meirama open pit, provide a good variety of hydrological data based on the precipitation, the water usages, the calibrated parameters and the geology of the region.

In this respect, a trial and error method was employed to estimate the hydrological parameters, on a 6 hour time step system resolution, and it was emphasized at those materials closer to the surface of the soil, where the regional hydrogeology is considered more relevant.

In other respects, the results in the surface domain are compared with some of the measured discharges at different sub-catchments and superficial drainage sections at the Mero River (Figure 3 and Figure 7). Although the discharges were not measured continuously throughout the period of simulation, the analysis of the comparison of surface runoff results shows in general a quite good behavior. It can be seen that the predicted surface runoff in the sub-catchments, with natural surface flow regimes, is much better assessed than, for instance, in Section 6 (S6) downstream in the Barcés River. In any case, the main diversion upstream in the Barcés River is to the Cañas Treatment Plant for Drinking Water, which has a constant allocation of 100 liters/second. Any other diversions from the Barcés River are not still well known or registered. Although, it is also interesting to note the evolution of surface flow rates in Section 8 (S8) far upstream in the Barcés River, where a quite good behavior was obtained after calibrating the diversion capacity of the weir installed at the Pereira stream. Since 15-June-2011, this weir splits up the diverted flow into the flooding of the Meirama Lake and the environmental flow supplied downstream through a tunnel to the Barcés River. This former flow was calibrated for a diversion capacity of 18 meters/day, in terms of a flow rate being diverted per unit area of a flooded surface with 10 meters of water depth (Equation (1)).

The flooding of the open pit began on 18-March-2008 with groundwater flow only. Meanwhile, the surface water drained by the Meirama basin was diverted to the Barcés River. After 3-October-2008 the mine is being flooded with water diverted from some of their drainage streams, including the Porta Antiga and the Pereira streams (Figure 5), and this was done until 11-June-2011, when the diversion is shared with the Barcés River in order to better satisfy its environmental flow. The evolution of the flooding, illustrated in Figure 8 (up left), is depicted from the point of view of the comparison between the registered and the assessed water levels during the simulation period (it includes calibration and validation stages), as well as that of the evolution of the relative errors for the maximum water depth (about 150 meters). Although the general behavior of the flooding is

quite acceptable, some of the differences are understood as a consequence of the variability of the surface and ground water hydrology, focused on the upper basin of the Barcés River, where the flooding of the open pit mine of Meirama takes place.

In addition to the analysis of the hydrological evolution of the flooding in the open pit by means of the water elevations that have been reached, other water variables, like the cumulative volumes of discharges and hydrological components around the open pit mine, are also analyzed (Figure 5, Figure 8 (up right)). This cumulative volume analysis considers: the surface water flow, the groundwater flow, the overland flow, the rainfall, the evaporation from the free surface water and the bare soil, as well as the transpiration due to the local vegetation. The total balance of all the components of this particular hydrological cycle are compared with the volume of water associated with the simulated flooding evolution (surface water storage in the open pit). Through this analysis, it is possible to assess the water storage in the involved porous materials of the analyzed area as the difference between the surface water storage and the total balance. Then, the multiple sources of water getting into the open pit are mainly involved in filling up the pit lake and in saturating the surrounding materials, as can be seen by the behavior of the surface water cumulative volume. An averaged groundwater discharge of about 3.1 cubic hectometers per year, inside a perimeter of about 8 km (Figure 5), is in good agreement with the volume pumped within the mine during the last exploitation period, until 18-03-2008.

It is also of great interest in this study to analyze the regulation capacities of the outlet works installed in the dam of the Cecebre Reservoir, although it could bring to unavoidable cumulative errors in the reservoir water storage. On this regard, some calibrations have been required.

First of all, concerning the floodgates, the numerical model is obviously not sufficiently enhanced, neither from the point of view of the developed numerical approximation, nor from the point of view of the mesh refinement at the floodgates and the dam. Therefore, a calibration of parameters would be needed, in order to properly assess the discharging curves of the concerned sluice floodgates. The behavior of the floodgates is supported, firstly, onto the average openings of the five gates (transformed on the equivalent transient heights of the corresponding overflow gates), secondly, onto the coarse numerical approximation for the dam and the five gates shape (as if they were only one) and, thirdly, onto the continuously computed water heads, hopefully appropriated anytime around the dam gates. In this context, and after studying the general behavior of the modeled floodgate, with the help of the water levels reached in the reservoir and the discharged flow rates, the averaged

openings (Equation (4)) were corrected by the following function that proved to be satisfactory:; where, h and h^* are, respectively, the corrected and the estimated mean openings (Figure 6) of the vertical lift gate. This function carries out the calibration adjustment between the actual and the modeled openings of the dam floodgates. As can be seen, it works in a way that increases the openings above 0.1 m, and reduces the openings below this value, a feature that was necessary after performing the calibration analysis.

Afterwards, concerning the pipe spillways of the dam, it looked like the dam needed an additional correction for the calculated flow rates. After studying the general behavior of the discharged flow rates when only the pipe spillways were operative, it had been found that there was a general maximum leakage in the modeled dam, which was assessed as about −0.08 m³/s with respect to the modeled middle pipe spillway mouth (Equation (3), the adjustment coefficient for the hydraulic head dependency was established as 0.004).

After the calibration of the behavior of the reservoir outlet works was considered satisfactory, for the estimated discharge flow rates at the dam, the comparison was made between the registered water levels reached in the reservoir, as provided by the water supply company of La Coruña (EMALCSA), and the simulated water levels, at two points of the reservoir, one downstream and the other upstream (Figure 8 (down)).

Despite how tough the calibration process of the parameters was for the involved devices, which were operating during almost five years (January-2008 to September-2012), it looks like the comparisons are in general in quite acceptable agreement. Nevertheless, it should be noticed that when surface flows are involved in the reservoir, the present modeling approach gives higher surface level gradients than those which could be hydraulically expected. The main reason for this is because the simulated water levels are approximated by the water heads, which are vertically averaged and lumped with the groundwater heads. It can also be noted that during the last hydrological year (mainly 2012), there is a poorer agreement between the observed surface water levels and the simulated ones. This is probably because during this period there was quite a severe drought in the region, and the shortage of water resources would enhance the storage accumulated errors related to the simulation of the water management mainly concerned with the appropriate free surface water levels and the water components involved in the hydrology and the operating policy of this small Cecebre Reservoir. A part of the errors are specific to the numerical simulation itself, but we can not forget those related to the lack of information (data, usages and management) in the Mero River catchment.

Finally, a past scenario of the surface and groundwater layouts which have been simulated is depicted in Figure 9, where the high level of interaction and variability of the hydrology of surface and groundwater flows in the Mero River watershed is clearly shown during a specific hydrological event of the close to five years simulation period. Among these results we would like to underline, for instance, the drainage layout of surface runoff and the freshwater hydraulic heads around the Meirama open pit, and the Cecebre Reservoir.

The main conclusion of this analysis is based upon the fact that the conservation of the water balance is verified on both a local and global basis. This responsibility lies on the consistent surface and ground water model interactions and on the here implemented numerical conditions that allow an integrated analysis of the freshwater resources management at the scale of this regulated watershed.

5. Conclusions

Some new developments have been performed on a process-based integrated hydrological model of 2-D saturated subsurface and surface flows (MELEF). This numerical model for continental and coastal catchment applications couples surface and subsurface regional hydrology through a finite element resolution of the saturated

flows of fresh and saltwater. Hydrological properties, simulated conditions and numerical results are managed by GIS.

Recent numerical developments are presented in order to extend this finite element solution to new practical applications dealing with water uses and management of surface water in widely regulated watersheds. This comprises the management of surface water flow control through the capacity of diversion works on river beds, spillways and outflow operations of floodgates in weirs and reservoir dams.

One application to the freshwater management in the Mero River catchment (La Coruña, Spain), during a period of almost five hydrological years, was presented to illustrate a high variety of interactions between surface and subsurface water, where regional hydrology allows, among others, the water management and flooding of a deep surface water body (Meirama pit lake), as well as the regulation operations at the dam of the Cecebre Reservoir through its pipe spillways and floodgates.

Firstly, the comparison of observed and simulated runoffs, surface water levels and water balances seems quite appropriate with respect to the observed regulation capacities, the environmental flow rates and the flooding evolutions in the open pit lake. Secondly, the comparison between observed and simulated surface flow rates and free water levels around the Cecebre Reservoir seems quite acceptable.

It can be concluded that the numerical model MELEF, developed for the present joint surface and subsurface regional flows, can provide useful and quite precise results concerning all the freshwater resources management of relatively large and regulated watersheds, like that of the Mero River of about 250 km², with a complex, anisotropic and fractured geology, where regional hydrology allows the flooding of a deep surface water body and the practical water regulation in weirs of river beds and common reservoir outlet works as pipe spillways and floodgates.

Acknowledgements

The work reported in this paper has been funded by the Spanish Ministry of Education and Science CICyT (CGL2009-11258), the Xunta de Galicia (Maria Barbeito Program), the companies LIMEISA and EMALCSA, as well as the European Science Foundation.

Conflict of Interest

The authors declare that they have no conflict of interest.

Cite this paper

FranciscoPadilla,J. HoracioHernández,RicardoJuncosa,Pablo R.Vellando, (2015) A Numerical Solution for the Integrated Analysis of Water Resources Management: Application to the Mero River Watershed, La Coruña, Spain*. Journal of Water Resource and Protection,07,815-829. doi: 10.4236/jwarp.2015.710066

References

NOTES

^*Solution and application to integrated watershed management.

^#Corresponding author.