domain: spi math phys

text:
Luiz F T Fernandez

Murat Bronz
email: murat.bronz@enac.fr

Thierry Lefebvre
email: thierry.lefebvre@onera.fr

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Multi-Vehicle Simulation Framework for Heterogeneous Unconventional MAVs

The growing interest in using aerial vehicles for diverse missions creates the need for higher autonomy and dexterity. Unconventional configurations can potentially satisfy such requirements, at a cost of new challenges in design and control. Simulation environments are often employed to tackle such challenges while minimizing resources expenditure. In this paper, an existing simulation framework has been extended to allow for the modeling of new UAV configurations. Aerodynamic forces and moments generated by propellers and airframes are defined in separate functions to ease the understanding and allow future addition of new vehicle configurations. Incremental nonlinear dynamic inversion and vector field guidance methods are implemented as vehicle agnostic control and guidance solutions. Vehicle models are extracted and implemented into Paparazzi open-source autopilot system, which makes a seamless switch towards real flights and hardware in the loop tests, facilitating control law development. The simulator can also be used within vehicle and mission design environments. With that, the entire closed loop behavior of the system can be addressed from the design stage, allowing for the vehicle and fleet optimization while considering control and operation constraints. Finally,

INTRODUCTION

Simulators have become indispensable tools in aeronautical and robotics research. Their use is one of the factors that led to the increasing interest in deploying aerial robots for real life missions. In a fully computational environment, researchers can evaluate vehicle performance and different control strategies without expending resources to manufacture and test robots. By employing a simulator to train a learning based solution, Loquercio et al. [START_REF] Loquercio | Learning high-speed flight in the wild[END_REF] showed how quadcopters can safely and nimbly navigate through complex environments such as forests and disaster zones. Petracek et al. [START_REF] Petracek | New era in cultural heritage preservation: Cooperative aerial autonomy for fast digitalization of difficult-to-access interiors of historical monuments[END_REF] also used simulations to assure that two unmanned aerial vehicles (UAVs) could be used to document historical buildings without endangering them. Thanks to the successful application of UAVs in different tests or even in real life scenarios, there is an increasing interest in using them in different tasks. Such missions are increasingly complex, covering demands where time plays a major role, such as health [START_REF] Hiebert | The application of drones in healthcare and health-related services in north america: A scoping review[END_REF] related tasks, and the ones involving more interaction with humans, asking for more dexterity, such as construction [START_REF] Choi | An overview of drone applications in the construction industry[END_REF]. To successfully perform more complex missions, novel vehicle configurations are continually being explored. Winged UAVs [START_REF] Vourtsis | Wind defiant morphing drones[END_REF], capable of higher flight ranges, and different propulsion layouts [START_REF] Mancinelli | Dual-axis tilting rotor quad-plane design, simulation, flight and performance comparison with a conventional quad-plane design[END_REF], as the vehicles shown in Figure 1, capable of more dexterous flights, are often preferred. However, they create challenges with respect to accurate physics simulation. Precise aerodynamic modeling of airframes and propellers then becomes of paramount importance in order to correctly account for vehicle dynamics as it enables the as-sessment of vehicle stability, control effectiveness, and performance in various flight regimes. Such models are usually built with high order methods in the aeronautical engineering community, using computational fluid dynamics methods as in [START_REF] Koyuncuoglu | Cfd based multicomponent aerodynamic optimization for wing propeller coupling[END_REF] and [START_REF] Alvarez | High-fidelity modeling of multirotor aerodynamic interactions for aircraft design[END_REF]. However, the use of such computationally heavy tools remains impracticable at system simulation level for a given trajectory or full mission considering sensing and control systems, as a single point simulation can take hours to run. So, the challenge of accurately addressing aerodynamic behavior of unconventional UAVs within an acceptable computational cost for system simulation remains open.

In addition to the challenges associated with creating new vehicle configurations and defining accurate aerodynamic forces and moments, existing simulators often suffer from limited portability and adaptability. Many simulators are designed to run on specific operating systems or hardware, as NVIDIA GPUs, making them less accessible for researchers using different platforms. This can hinder collaboration and the ability to reproduce research findings across different research groups. In this paper, we present an extension to an existing simulation framework based on Bullet physics engine, which is operating system and hardware-agnostic. The simulator is intended to be used to model new vehicle designs, develop control and guidance laws and simulate mission scenarios while testing different vehicle configurations at the same time. The subsequent sections of this paper delve into the details of related work, proposed solution, and usecase scenarios, highlighting the benefits and potential applications of our approach. Contributions of the paper can be summarized as :

• Facilitate the simulation of new vehicle configurations • Easier and more intuitive definition of aerodynamic forces and moments

• Heterogeneous multi-vehicle simulation

• Improved portability through different operating systems

• Implementation into a well known open source autopilot system for rapid transition into hardware flights.

RELATED WORK

The desire for a simulation tool to evaluate different vehicle designs in [START_REF] Baskaya | A novel robust hexarotor capable of static hovering in presence of propeller failure[END_REF][START_REF] Luiz | Development of a missiontailored tail-sitter mav[END_REF], and novel control strategies for such vehicles, as in [START_REF] Barth | Model-free control algorithms for micro air vehicles with transitioning flight capabilities[END_REF], led to the search of available solutions. The interest in testing different guidance strategies [START_REF] Bilgin | Panel Method Based Guidance for Fixed Wing Micro Aerial Vehicles[END_REF] and even learning approaches [START_REF] Hovell | Learned multiagent real-time guidance with applications to quadrotor runway inspection[END_REF] increased even further the need for a reliable and flexible tool. The several available simulators for aerial vehicles are different in terms of objective, implementation, available vehicles, flexibility, and so on. So, the choice of evaluating existing solutions and selecting the most suitable one was made. Fadri et al. [START_REF] Furrer | Robot Operating System (ROS): The Complete Reference (Volume 1), chapter RotorS-A Modular Gazebo MAV Simulator Framework[END_REF] presented the Ro-torS, a simulator for different multicopter vehicles. The implementation in Robot Operating System (ROS) and Gazebo facilitates inclusion of sensors for odometry and an Inertial Measurement Unit (IMU). In the same vein, Baca et al. [START_REF] Baca | The mrs uav system: Pushing the frontiers of reproducible research, real-world deployment, and education with autonomous unmanned aerial vehicles[END_REF] presented the MRS UAV system, a ROS-based solution to simulate different vehicles while enabling real world deployment and facilitating reproducible research. Yunlong et al. [START_REF] Song | Flightmare: A flexible quadrotor simulator[END_REF] presented the Flighmare, a flexible quadrotor simulator based on Unity and a physics engine that allows for the inclusion of sensing system on the simulation, evaluation of reinforcement based strategies for control, and more. From the same research group, Foehn et al. [START_REF] Foehn | Agilicious: Open-source and open-hardware agile quadrotor for vision-based flight[END_REF] presented Agilicious, an open-source and open-hardware agile quadrotor for visionbased flights. This solution comes with agilib, a library with minimal dependencies to integrate Agilicious into a given infrastructure, and agiros, a ROS based library that allows for easy simulation and real flight. Fernandez-Cortizas et al. [START_REF] Fernandez-Cortizas | Aerostack2: A software framework for developing multi-robot aerial systems[END_REF] presented the Aerostack2, a ROS2 based simulator for multiaerial-robots systems. Keipour et al. [START_REF] Keipour | A simulator for fully-actuated uavs[END_REF] presented a Matlab based simulator for fully-actuated vehicles. With the main focus of studying learning based methods, Kulkarni et al. [START_REF] Kulkarni | Aerial gym -isaac gym simulator for aerial robots[END_REF] presented the Aerial Gym Simulator, based on the NVIDIA Isaac Gym. Panerati et al. [START_REF] Panerati | Learning to fly-a gym environment with pybullet physics for reinforcement learning of multi-agent quadcopter control[END_REF] also presented a Python based simulator adapted for learning methods. The authors introduced an environment called Gym-Pybullet-Drones, based on the Bullet physics engine, which supports multiple quadcopters and provides realistic collisions, aerodynamic effects, and reinforcement learning interfaces. With the objective of helping in education, Folk et al. [START_REF] Folk | Rotorpy: A python-based multirotor simulator with aerodynamics for education and research[END_REF] presented the RotorPy, a python based multirotor simulator that considers aerodynamic effects. More simulation frameworks and research directions can be found in the extensive review by Mairaj et al. [START_REF] Mairaj | Application specific drone simulators: Recent advances and challenges[END_REF]. Considering the needs of adding different vehicle configurations and aerodynamic and propeller modeling, as well as some aspects of the available simulators, such as maturity, ease of integration with the Paparazzi autopilot system and design environment, we opted to use the Gym-Pybullet-Drones, extending it according to our interest.

EXTENSIONS

Simulating heterogeneous fleets was one of the main objective for the improvement of the existing framework. However, a unified controller, suitable with most of the new and unconventional vehicle configurations (such as coaxial birotors, fully-actuated hexacopters, and tail-sitters), is also implemented. Some of these vehicles require higher-fidelity models to correctly represent aerodynamic forces and moments. Additionally, for the purpose of mission simulations, a more generic guidance algorithm becomes necessary. All of these issues have been addressed and presented below. This tool is under active development. So new vehicles, control strategies and improvements are added in an ongoing basis.

Vehicle Agnostic Controller

The controller is based on the work of Van Wijngaarden et al. [START_REF] Van Wijngaarden | Flight code convergence: fixedwing, rotorcraft, hybrid[END_REF] and revolves around the control of the angular accelerations in an incremental way. The attitude control law cal-culates the new command (u c ) as a function of the current command (u f ), current angular acceleration ( Ωf ) and thrust (T f ) as:

u c = u f + G + (ν -Ωf T f ) (1)
where G + is the Moore-Penrose pseudoinverse of the control effectiveness matrix G and ν is a virtual control law that uses a proportional feedback (k Ω ) to control the angular rates as:

ν = k Ω (Ω ref -Ω) T d (2)
where T d is the desired thrust, calculated by the outer loop, and Ω ref is calculated as a function of the quaternion attitude error (q err ) , obtained with the Kronecker product between the reference quaternion and the conjugate of the state quaternion:

Ω ref = k η q T err (3)
The velocity control law is:

v = v f + m [G T (η, T ) + G L (η, V )] -1 (a ref -a f ) (4)
where a ref and a f are the reference and current accelerations respectively, v f is the current velocity, and G T (η, T ) and G L (η, V ) are velocity effectiveness matrices for thrust and lift forces, obtained with the derivatives of such forces with respect to the control vector for the outer loop, defined as roll, pitch, and thrust ([ϕ, θ, T ]). The subscript [.] f represents the filtered signal for the current acceleration and velocity vectors. The Bullet physics engine runs at 240 Hz, while the controls run at 96 Hz.

Vector Field Guidance

A generic guidance law becomes beneficial when a wide variety of vehicle configurations needs to be simulated. So we implemented the Vector Field Guidance method from De Marina et al. [START_REF] Garcia | Guidance algorithm for smooth trajectory tracking of a fixed wing uav flying in wind flows[END_REF]. The basic guidance law uses the velocity error to generate the reference acceleration as:

a ref = [v ref -v] k v (5)
where, v is the current velocity, k v is decomposed as speed and heading gains for tangential and radial directions, and v ref is the reference velocity coming from a predefined vector field. The acceleration terms are neglected, and only velocity information is used. The implementation requires an analytical formulation of the trajectory that is used as a reference, e.g., a circle x 2 + y 2 + r 2 = 0. At any point in the environment, a "level-set" can be calculated which has a non-Euclidean distance metric with respect to the reference trajectory, and by driving this distance error towards zero, the vehicle converges to the reference trajectory. Figure 2 shows an example of the outcome of the vector field guidance. The desired velocity (direction and module) for a vehicle at each point in the space (represented by the arrows) is calculated in such a way that the desired trajectory, represented by the circle with radius 100 m, is correctly tracked.

Propeller Aerodynamics

To calculate in real time the thrust and moments generated by propellers assuming different geometries and flight including forward flight, four methodologies are implemented. Both low order data-driven methodologies proposed by Gill and D'Andrea [START_REF] Gill | Computationally efficient force and moment models for propellers in uav forward flight applications[END_REF] are implemented and indicated when experimental data for thrust and moments for the specific propeller are known. The hybrid blade element momentum (HBEM) by Davoudi and Duraisamy [START_REF] Davoudi | A hybrid blade element momentum model for flight simulation of rotary wing unmanned aerial vehicles[END_REF] is implemented to account for cases where no performance data is known, only the propeller geometry. In order to use the latter, a valid Matlab license is needed. Lastly, surrogate models can also be used to account for higher fidelity models, enabling a compromise between computational cost and accuracy. Regardless of the surrogate technique or generation strategy, any Pickle file containing the model can be read and used for the calculation. Surrogate models created using CCBlade [START_REF] Ning | Using blade element momentum methods with gradient-based design optimization[END_REF] were used alongside the simulator in [START_REF] Luiz | Assessment of closed loop dynamics in the multidisciplinary design and optimization of small uavs[END_REF] to simulate and optimize a tail-sitter in cruise flight. In [START_REF] Luiz | Assessment of methods for propeller performance calculation at high incidence angles[END_REF] seven different methodologies for propeller performance calculation were compared using the proposed simulator. The goal was to observe if, by changing the calculation method in a quadcopter trajectory simulation, the impact in the control input would be noticeable. Despite the fact that desired trajectory was correctly tracked by the vehicle independently of the propulsion methodology, we obtained a significant difference in motor rpm commands, as shown in Figure 3, that compares both methods by Gill and D'Andrea [START_REF] Gill | Computationally efficient force and moment models for propellers in uav forward flight applications[END_REF] and HBEM [START_REF] Davoudi | A hybrid blade element momentum model for flight simulation of rotary wing unmanned aerial vehicles[END_REF]. The HBEM led to smaller rpm values, which could lead to a smaller battery requirement during the design phase of a vehicle.

Vertical takeoff and landing UAVs

Modeling the aerodynamic forces and moments experienced by vertical takeoff and landing (VTOL) UAVs remains an open challenge. This kind of vehicle flies in a larger flight envelope than general aviation aircraft, including post-stall conditions. To account for that, the Φ-theory based model, as proposed by Lustosa et al. [START_REF] Lustosa | Global singularity-free aerodynamic model for algorithmic flight control of tail sitters[END_REF] was implemented. This model is algebraically simple and does not rely on angles of attack and sideslip, thus becoming singularity-free, while taking into account the influence of the wing, the deflection of wing control surfaces, and propeller-wing interaction. The same model has already being used by Tal et al. [START_REF] Tal | Global incremental flight control for agile maneuvering of a tailsitter fly-ing wing[END_REF] to develop a control law capable of enabling agile flight of tailsitters. In order to make this module applicable for any flying wing configuration, the contributions generated by the propeller-wing interaction in the forces and moments are not taken into account and left for future work, as the original formulation is specific for vehicles with two propellers. Below, we provide a brief explanation of how the aerodynamic forces and moments are generated. For further information, we refer the reader to [START_REF] Lustosa | Global singularity-free aerodynamic model for algorithmic flight control of tail sitters[END_REF]. Considering a flying wing UAV with 2 elevons and assuming the use of thin airfoils, three Φ-coefficients are needed to build the model using the Φtheory: the velocity-to-force coefficient Φ (f v) , the velocityto-moment coefficient Φ (mv) , and the angular-velocity-tomoment coefficient Φ (mω) :

Φ (f v) =   C D0 0 0 0 C Y0 0 0 0 2π + C D0   (6)
where C D0 and C Y0 are the minimum drag and side force coefficients, respectively.

Φ (mv) =   0 0 0 0 0 -c -1 ∆r(2π + C D0 ) 0 b -1 ∆ r C Y0 0   (7)
where b, c, and ∆ r are the wing span, wing mean aerodynamic chord, and flap chord percentage, respectively.

Φ (mω) =   C lp C lq C lr C mp C mq C mr C np C nq C nr   (8)
where each coefficient accounts for the variation of the rolling (l), pitching (m), or yawing (n) moment coefficient with dimensionless rate of change of roll p, pitch q, or yaw r rate [START_REF] Roskam | Airplane Flight Dynamics and Automatic Flight Controls. Number p. 1 in Airplane Flight Dynamics and Automatic Flight Controls[END_REF]. Additionally, it is necessary to model the so-called elevon effectiveness for both aileron and elevator configurations ζ f and ζ m ∈ R 3 . With that, the aerodynamic force on the body frame (F b ) at the aerodynamic center is calculated as:

F b = - 1 2 ρSΦ (f v) vv b + 1 4 ρSΦ (f v) [ζ f ×](δ 1 + δ 2 )vv b - 1 2 ρSΦ (mv) vBω b + 1 4 ρSΦ (mv) [ζ f ×](δ 1 + δ 2 )Bvv b (9)
where [ζ f ×] represents the vector product operation matrix of ζ f , ρ is the air density, S the wing area, δ 1 and δ 2 are the deflections of each elevon, v is the airspeed, v b and ω b are linear and angular speed in the body frame, respectively, and B a matrix defined as a function of wing span (b) and mean aerodynamic chord (c) as:

B =   b 0 0 0 c 0 0 0 b   (10)
The aerodynamic moment (M b ) acting on the body frame at the aerodynamic center is calculated as:

M b = - 1 2 ρSBΦ (mv) vv b - 1 2 ρSBΦ (mw) vω b + 1 4 ρSBΦ (mv) [ζ m ×](δ 1 + δ 2 )vv b + 1 4 ρSBΦ (mω) [ζ m ×](δ 1 + δ 2 )vBω b (11)
Figure 4 shows a visual comparison between a real and a simulated UAV.

Adding new extensions

Within the simulation environment, every robot is defined in a Unified Robot Description Format (URDF) file. So, defining a new vehicle and its "physics" consists in creating a new URDF file and defining a new function in the base aviary. In this new function, each force and moment acting on the vehicle has to be declared, as well as where (or which link) it is applied. Every existing physics can also be updated to account for new methods, such as a new methodology for propeller performance. Control and guidance strategies are implemented in the control aviary, so new solutions should be added to this Python script.

USE-CASE SCENARIOS

In the preceding sections, we have outlined the essential building blocks of the framework necessary for simulating a comprehensive use-case. In this section, we present potential applications for the tool. The primary aim is to demonstrate the framework's potential and provide guidance and motivation for other researchers by offering ready-to-use templates to initiate their own investigations.

Drone Racing : Time Optimal Trajectory Tracking

Drone racing has emerged as an exciting and rapidly growing sport, captivating both enthusiasts and researchers alike. With its unique blend of speed, agility, and precision, drone racing poses unique challenges that require advanced control strategies and trajectory tracking techniques. In the example shown in Figure 5, seven gates have been positioned in order to create a figure-eight maneuver. The drone starts from the middle section of the track, and calculates the fastest trajectory through the gates with a predefined sequence. The trajectory is generated by minimizing the snap of the polynomial trajectory, passing through the center of the gates. It can be seen that the desired trajectory calculated by the minimum snap (blue), and the actual flown trajectory (red) has a significant discrepancy. Other planning and control strategies can be implemented and evaluated in this context before real flights.

Vehicle Design and Control Law Optimization

The proposed simulation framework was used in a multidisciplinary optimization process applied to a tail-sitter vehicle. A UAV designed in [START_REF] Luiz | Development of a missiontailored tail-sitter mav[END_REF] for the IMAV20221 presented an undesired overshoot under high wind conditions, as shown in Figure 6. So, in [START_REF] Luiz | Assessment of closed loop dynamics in the multidisciplinary design and optimization of small uavs[END_REF], the entire vehicle simulation was embedded in an optimization problem, where the energy consumption for a circular lap was minimized, while reinforcing its wind robustness via constraints. The propulsion layout of the vehicle and the guidance control law gains were the designs variables for the problem. With the process, a more energy efficient and robust to wind vehicle was obtained. Figure 7 shows the differences in the flight path of the baseline and optimized vehicle.

IMPLEMENTATION TO PAPARAZZI AUTO PILOT SYSTEM

The integration of a simulator environment with a real drone autopilot system, such as Paparazzi [START_REF] Hattenberger | Using the Paparazzi UAV System for Scientific Research[END_REF], offers numerous advantages for the development and verification of control and guidance algorithms. By enabling the seamless trans- fer of code between simulation and real-world environments, this approach provides a powerful tool for iterative development, testing, and validation processes. The current develop- ment state allows for the use of Paparazzi's ground control station to simulate an entire mission. The bridge between simulation and Paparazzi enables the full replacement of the real system with the simulation, while assuming perfect sensing, for the moment. By harnessing the power of simulation, advancements in drone technology can be accelerated, leading to improved aerial capabilities and expanded applications across various industries.

CONCLUSION

We presented the extension of an existing framework for simulating aerial robots adapted to account for different vehicle configurations, control and guidance strategies. This open-source solution can be used during design phase, evaluating unconventional vehicle configurations and their perfor-mance with respect to a given mission. It can also be used for control and guidance laws development, testing, and benchmark. Its simple structure and the usage of Python allows for easy implementation of new features. This tool has already been used for different applications, such as drone racing and tail-sitter design and optimization. Ultimately, it will be used as well for simultaneous mission and heterogeneous swarm of vehicles design.

Figure 1 :

 1 Figure 1: Different vehicle configurations flying together are shown in the simulation environment.

Figure 2 :

 2 Figure 2: Illustration of vector field guidance, the arrows indicate the velocity direction at each point in the plane.

Figure 3 :

 3 Figure3: Discrepancies found in rpm commands for a quadcopter for different propeller methods: HBEM and the ones by Gill[START_REF] Luiz | Assessment of methods for propeller performance calculation at high incidence angles[END_REF].

Figure 4 :

 4 Figure 4: Simulated and real vehicle.

Figure 5 :

 5 Figure 5: Simulation of a drone passing through gates, representing an example drone racing scenario, with the desired trajectory in blue and the simulated one in red.

Figure 6 :

 6 Figure 6: Flight trajectory in the IMAV2022 competition: overshoot in one of the turns caused by strong winds.

Figure 7 :

 7 Figure 7: Flight trajectories of baseline and optimized vehicles under wind conditions, with error bounds defined as ±30[m].

Figure 8 :

 8 Figure 8: Simulation model running by Paparazzi Autopilot system.

https://www.imavs.org/2022

ACKNOWLEDGEMENTS

This work is part of the activities of ONERA-ISAE-ENAC joint research group. The authors would like to thank the ENAC Research Engineer Fabien Bonneval for writing the Paparazzi-pyBullet wrapper that is used in this work.


Benoît Dewandel

Sandra Lanini

Nicolas Frissant

Steady-state semi-analytical

Keywords: Underground dam, Subsurface dam, Steady-state semi-analytical solution, Fundamental solutions method, Managed Aquifer Recharge

teaching and research institutions in France or abroad, or from public or private research centers.

Steady-state semianalytical solutions for assessing the two-dimensional hydraulic head distribution induced by an underground dam in a sloping aquifer with artificial aquifer recharge and pumping Benoît Dewandel, Sandra Lanini, Nicolas Frissant To cite this version:

Benoît Dewandel, Sandra Lanini, Nicolas Frissant. Steady-state semianalytical solutions for assessing the two-dimensional hydraulic head distribution induced by an underground dam in a sloping aquifer with artificial aquifer recharge and pumping. Hydrogeology Journal, 2023, 10.1007/s10040-023-02734-2. hal-04314128

Introduction

Many countries face growing pressure on their water resources, with increasing demands for domestic, agricultural and industrial water supplies, coupled with a reduced resource due to climate change. This has created the need for artificially augmenting supplies from existing resources, and the development of groundwater technologies (e.g. [START_REF] Bouwer | Artificial recharge of groundwater: hydrogeology and engineering[END_REF][START_REF] Dillon | Future management of aquifer recharge[END_REF]. Underground (or subsurface) dams (Fig. 1) are one of such technologies [START_REF] Nishigaki | Underground dam technology in some parts of the world[END_REF][START_REF] Onder | Underground dam: A tool of sustainable development and management of groundwater resources[END_REF][START_REF] Izady | Optimal water allocation from subsurface dams: a risk-based optimization approach[END_REF]). An underground dam, usually-but not only (e.g., in karst aquifers, [START_REF] Stevanović | Damming underground flow to enhance recharge of karst aquifers in the arid and semi-arid worlds[END_REF][START_REF] Ahmed | Groundwater dams, General characteristics and historical development[END_REF][START_REF] Milanović | Engineering Karstology of dams and reservoirs[END_REF])-constructed in an alluvial area, often consists of a hydraulic barrier that may partially or completely penetrates the aquifer up to the impermeable bedrock. Such barriers may be semi-permeable or impermeable and may also not reach bedrock; the reader can refer to [START_REF] Kaleris | The effect of cutoff walls on saltwater intrusion and groundwater extraction in coastal aquifers[END_REF] and [START_REF] Wu | Preventing seawater intrusion and enhancing safe extraction using finitelength, impermeable subsurface barriers: 3D analysis[END_REF] for other dam configurations. The dam modifies groundwater flow and the groundwater level rises upstream and declines downstream of the dam. Upstream, the water can thus be abstracted and, ideally, used throughout the year. In coastal aquifers, many work has been done to investigate their use for alleviating seawater intrusion to increase fresh groundwater resources (e.g. [START_REF] Pool | Dynamics of negative barriers to prevent seawater intrusion[END_REF][START_REF] Kaleris | The effect of cutoff walls on saltwater intrusion and groundwater extraction in coastal aquifers[END_REF][START_REF] Abdoulhalik | How does layered heterogeneity affect the ability of subsurface dams to clean up coastal aquifers contaminated with seawater intrusion[END_REF], Sun et al. 2019[START_REF] Wu | Preventing seawater intrusion and enhancing safe extraction using finitelength, impermeable subsurface barriers: 3D analysis[END_REF][START_REF] Zheng | Insights of variable permeability full-section wall for enhanced control of seawater intrusion and nitrate contamination in unconfined aquifers[END_REF][START_REF] Zhang | Effects of a low-permeability layer on unstable flow pattern and land-sourced solute transport in coastal aquifers[END_REF]).

Underground dams have been used for a long time in arid and semi-arid zones, where groundwater is not available in sufficient quantity for part of the year, to retain water in alluvial aquifers [START_REF] Bceom | The underground dams, in French. Les barrages souterrains[END_REF][START_REF] Nilsson | Groundwater dams for small-scale water supply[END_REF][START_REF] Nishigaki | Underground dam technology in some parts of the world[END_REF][START_REF] Raju | Subsurface dams to harvest rainwater-a case study of the Swarnamukhi River basin, Southern India[END_REF][START_REF] Şen | Reservoirs for Water Supply Under Climate Change Impact-A Review[END_REF]. This technique was greatly improved and rendered more complex from the 1980s onward in China [START_REF] Sun | Impacts of artificial underground reservoir on groundwater environment in the reservoir and downstream area[END_REF], Japan [START_REF] Furukawa | The underground dam, a new concept in water resource development to meet the challenges of natural disaster and climate change[END_REF] and South Korea [START_REF] Ishida | Review: Sustainable use of groundwater with underground dams[END_REF], with increasingly long (several kilometres) and high (several tens of meters) dams.

Underground dams can be combined with managed aquifer recharge (MAR) systems, such as injection wells or percolation ponds, for improving the efficiency of the overall setup [START_REF] Bceom | The underground dams, in French. Les barrages souterrains[END_REF], Sun et al. 2019[START_REF] Choi | Effects of hybrid-type artificial groundwater recharge and underground barrier in a small basin[END_REF]. The vertical concealment of the aquifer by the underground dam is very variable, on the other hand, in most cases, the aquifer is concealed over its entire aquifer width and therefore over the entire width of the alluvial valley (e.g. [START_REF] Ishida | Construction of subsurface dams and these impacts on the Environment[END_REF], Sun et al. 2019[START_REF] Chang | Effect of subsurface dams on saltwater intrusion and fresh groundwater discharge[END_REF][START_REF] Zheng | Insights of variable permeability full-section wall for enhanced control of seawater intrusion and nitrate contamination in unconfined aquifers[END_REF]. A few cases of partial lateral occultation are apprehended by [START_REF] Wu | Preventing seawater intrusion and enhancing safe extraction using finitelength, impermeable subsurface barriers: 3D analysis[END_REF] and [START_REF] Frissant | État des connaissances sur le fonctionnement de deux systèmes aquifères du territoire du SAGE des Gardons et propositions d'amélioration de la connaissance : Les alluvions du Gardon et les molasses miocènes du bassin d'Uzès et les calcaires urgoniens de la Fontaine d'Eure[END_REF]..

Regardless of dam geometry, such subsurface structures act as a barrier to flow, modify the natural groundwater flow which may have an environmental, economic and/or social impact (basement flooding, soil salinization, seawater intrusion, etc., [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF]Yihdego 2016;[START_REF] Attard | Review: Impact of underground structures on the flow of urban groundwater[END_REF]). Any modification of natural groundwater flow has thus to be studied and anticipated, not only during construction, but also over the long term.

Much theoretical work was based on the analysis of numerical models for estimating their impact on groundwater flow (Merrick and Jewel 2003[START_REF] Paris | Hydrological effects of bounding the Venice (Italy) industrial harbour by a protection cut-off wall: a modeling study[END_REF][START_REF] Khairy | Analysis and design of a deep subsurface dam[END_REF][START_REF] Senthilkumar | Modelling the impact of a subsurface barrier on groundwater flow in the lower Palar River basin, southern India[END_REF][START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF][START_REF] Attard | Underground constructions impacts on groundwater flow in urban context[END_REF][START_REF] Kim | Validity evaluation of a groundwater dam in Oshipcheon River, eastern Korea using a SWAT-MODFLOW model[END_REF]. Others have explored various techniques to mitigate seawater intrusion to obstruct the inland movement of the saltwater wedge (e.g. [START_REF] Kaleris | The effect of cutoff walls on saltwater intrusion and groundwater extraction in coastal aquifers[END_REF][START_REF] Li | Laboratory tests and numerical simulations on the impact of subsurface barriers to saltwater intrusion[END_REF][START_REF] Wu | Preventing seawater intrusion and enhancing safe extraction using finitelength, impermeable subsurface barriers: 3D analysis[END_REF][START_REF] Zheng | Insights of variable permeability full-section wall for enhanced control of seawater intrusion and nitrate contamination in unconfined aquifers[END_REF][START_REF] Zhang | Effects of a low-permeability layer on unstable flow pattern and land-sourced solute transport in coastal aquifers[END_REF]. With regard to the impact on groundwater flow and therefore the distribution of induced hydraulic head -the subject of this study, some of this work has resulted in useful semiempirical equations on the long-term impact of disturbances caused by underground structures (tunnels, basements of buildings, deep foundations, etc.; [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF]. Few conventional analytical equations are available [START_REF] Marinos | Rise of the groundwater table when flow is obstructed by shallow tunnels[END_REF]Kavvadas 1997, Deveughèle and[START_REF] Deveughèle | Impact of an impervious shallow gallery on groundwater flow[END_REF], but, according to [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF], their validity remains limited. [START_REF] Fitts | Analytical modelling of impermeable and resistant barriers[END_REF] used analytical elements under steady-state condition for modelling impermeable and resistant barriers to groundwater flow, but they were considered with zero thickness (see also [START_REF] Strack | Principles of the analytic element method[END_REF]. [START_REF] Butler | Pumping test in non-uniform aquifers, the radially symmetric case[END_REF] described a transient analytical solution for the drawdown induced by a pumping well in a uniform matrix near a cylinder characterized by different hydrodynamic properties. This can be considered-if low hydraulic-conductivity values are used-as an obstacle to groundwater flow. El-Hames (2012) proposed a numerical solution for evaluating the transient of a water-table rise upstream of a dam, based on a general flow equation. However, the dam was of infinite-length and the solution cannot be used to assess the hydraulic head downstream of the dam. Therefore, no analytical solution has been proposed for evaluating the two-dimensional hydraulic head distribution induced by a rectangular underground dam, even for the simplest case, i.e. a finite-dimensional dam penetrating entirely into an infinite (or semi-infinite) aquifer. This is the subject of this study.

This study seeks alternative steady-state (semi-)analytical solutions based on the method of fundamental solutions (MFS), for computing hydraulic-head disturbance created by a fully penetrating impermeable dam of finite-size (rectangular) in a semi-infinite sloping aquifer. This technique is used in the field of physical and engineering problems to overcome the difficulty of finding an analytical solution in complex systems [START_REF] Chaiyo | The method of fundamental solutions for solving free boundary saturated seepage problem[END_REF][START_REF] Young | Method of fundamental solutions for multidimensional Stokes equations by the dual-potential formulation[END_REF]). In the field of groundwater hydraulics, MFS has been rarely used (Wang andZheng 2015, Neto et al. 2020). [START_REF] Kuo | Predicting pressure distribution in irregularly shaped reservoirs[END_REF], [START_REF] Yeung | flow reduction due to hydraulic barrier engineering structure: case of urban area flood, contamination and pollution risk assessment[END_REF] and [START_REF] Kuo | An image well method for predicting drawdown distribution in aquifers with irregularly shaped boundaries[END_REF] used this technique to propose drawdown solutions caused by wells pumping inside irregularly shaped reservoirs, and [START_REF] Huang | Estimating stream filtration from a meandering stream under the Robin condition[END_REF] for estimating flux filtration on a meandering stream caused by pumping wells. In this study, it is used to solve the problem of finite-dimensional obstacles to groundwater flow. MFS is a meshless method based on the fundamental solutions of governing equations [START_REF] Kupradze | The method of functional equations for the approximate solution of certain boundary value problems[END_REF]Aleksidze 1964, Golder 1995), which gives it certain advantages over conventional models (no model to build, no mesh generation, no discretization of physical boundaries). Its basic concept is to decompose the solutions of a partial differential equation by superposition of the fundamental solution with proper strengths, depending upon the location of boundary conditions. It assumes that a linear combination of the strengthened appropriate fundamental solution approximates the unknown solution sought. The accuracy of the method depends only on the distances between pairs of points where the fundamental solution is applied. In groundwater hydraulics, but not only, fundamental solutions are typically line-sink or point-source solutions.

Results of the proposed semi-analytical solution for an underground dam alone are analysed, compared to numerical modelling and discussed. Graphs with dimensionless parameters, relating the dam geometry to maximum hydraulic-head variations up-and downstream of the dam, and hydraulic-head variations compared to a condition without a dam are given. Analytical results are compared to two field examples that correspond to the theoretical situation, one of the Massillargues-Attuech site (France, Ladouche and Frissant, 2022) and the other from the literature [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF]. Theoretical formulation is given in case of a pumping/injecting well and/or a managed artificial recharge (MAR) near the underground dam, to examine the behaviour of composite devices (Fig. 1b).

The proposed solutions can be implemented as operational tools for engineers designing underground dams. It is not suggested that the proposed solutions replace existing analytical or numerical hydrogeological models, used for modelling disturbances induced by the building of underground structures. Rather, they are meant to supplement existing models for improving the design of such engineered systems.

Mathematical statements

Hydraulic head solution for an underground dam

This study considers an unconfined, sloping, homogeneous aquifer with a saturated thickness without dam D0 (in m; Fig. 2a). The aquifer dip is aligned with the x-axis direction and θ is the dip angle of the impermeable basement (in radian). The natural flow direction is parallel to the aquifer bottom, the natural (no dam) hydraulic gradient, i, being equal to tanθ. The aquifer is characterized by a constant and uniform hydraulic conductivity (K in m/s) and is limited in space downstream by a constant-head boundary at a distance 'd' (in m) from a Cartesian coordinate system (x,y), resulting in a semi-infinite aquifer. The underground dam has a rectangular shape, fully penetrates the aquifer and is impermeable, i.e. there is no flow through the dam. It is located at the centre of the coordinate system, and is defined by its length along the x-axis (Lx in m), and its length along the y-axis (Ly in m).

The governing steady-state groundwater flow equation is defined by the following non-linear Boussinesq equation: Assuming that the aquifer slope is small (i=tanθ ≈θ; θ≤0.02), which implies that groundwater flow is parallel to the substratum, a linearized form of the previous equation can be found by defining 𝑠𝑠(𝑥𝑥, 𝑦𝑦) = 𝐾𝐾(𝑥𝑥, 𝑦𝑦) -𝐾𝐾 0 and assuming that (𝐾𝐾 -𝐾𝐾 0 ) 𝐾𝐾 0 ⁄ is small (i.e. the first method of linearization is used; e.g. Zlotnick et al. 2017). Therefore, Eq.1 can be rewritten as:

𝐾𝐾𝐾𝐾 0 � 𝜕𝜕 2 𝑠𝑠 𝜕𝜕𝜕𝜕 2 + 𝜕𝜕 2 𝑠𝑠 𝜕𝜕𝜕𝜕 2 � -𝐾𝐾. 𝑖𝑖 𝜕𝜕𝑠𝑠 𝜕𝜕𝜕𝜕 = 0 (2)
with 𝑠𝑠(𝑑𝑑, 𝑦𝑦) = 0 and 𝑠𝑠(-∞, 𝑦𝑦) = 0.

The uncertainty introduced by the linearization method used is discussed later.

The no-flux condition through the four sides of the underground dam can be represented by applying the method of fundamental solutions. This assumes that the solution of the equation describing the hydraulic head in the sloping aquifer with an underground dam can be obtained from a linear combination of strengthened terms produced by application of the appropriate fundamental solution at a series of image wells. Along the dam sides, the superposition of the perturbations induced by each image well results in a zero flow.

The fundamental solution of the linearized equation (Eq.2) describing the steady-state hydraulic head in an inclined aquifer, bounded by a constant head boundary (𝐾𝐾 𝑜𝑜 𝑎𝑎𝑎𝑎 𝑥𝑥 = 𝑑𝑑) with a pumping well, can be found in Hantush (1964a & b). With the used coordinate system (Fig. 2a), the form of the fundamental solution is given by the saturated thickness D:

𝐾𝐾(𝑥𝑥, 𝑦𝑦) -𝐾𝐾 0 = 𝑄𝑄 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 e 𝛼𝛼𝜕𝜕 [K 0 (𝑟𝑟𝑟𝑟) -K 0 (𝑟𝑟 𝑛𝑛 𝑟𝑟)] = 𝑄𝑄. 𝐹𝐹(𝑥𝑥, 𝑦𝑦) (3)
where Q is the pumping rate (Q <0; in m 3 /s), K0(u) is the modified Bessel function of the second kind of the zero order, and e is the exponential function, 𝑟𝑟 = 𝑖𝑖 (2𝐾𝐾 0 ) ⁄ , 𝑟𝑟 = �𝑥𝑥 2 + 𝑦𝑦 2 and 𝑟𝑟 𝑛𝑛 = �(2𝑑𝑑 -𝑥𝑥) 2 + 𝑦𝑦 2 . Note that the same solution can be retrieved in [START_REF] Py | Theory of Ground Water Movement[END_REF]p.464-465).

According to the method of fundamental solutions, the hydraulic head in the sloping aquifer with an underground dam is given by:

ℎ(𝑥𝑥, 𝑦𝑦) = ∑ 𝛾𝛾 𝑗𝑗 . 𝐹𝐹 𝑗𝑗 �𝑥𝑥, 𝑦𝑦, 𝑥𝑥 𝑗𝑗 , 𝑦𝑦 𝑗𝑗 � 𝑛𝑛 𝑗𝑗=1 + 𝐾𝐾 0 + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) (4)
where x, y are the coordinates (in m) where the hydraulic head h has to be calculated, and xj and yj the coordinates of the j th image wells (up to n), the fundamental solution Fj is calculated to represent the underground dam, and γj are strengths applied to this solution that must be evaluated (Fig. 2b). If image wells are regularly located within the dam at a distance d' of the dam sides, the distance between them is: 𝜎𝜎 = 2 (𝐿𝐿𝑥𝑥 + 𝐿𝐿𝑦𝑦 -4𝑑𝑑′) 𝑛𝑛 ⁄ . Note that, when γj =0, Eq. 3 corresponds to the hydraulic head without dam, or the natural flow condition.

Incorporating Eq.3 in Eq.4 yields:

ℎ(𝑥𝑥, 𝑦𝑦) = ∑ 𝛾𝛾 𝑗𝑗 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 e 𝛼𝛼�𝜕𝜕-𝜕𝜕 𝑗𝑗 � �K 0 �𝑅𝑅 𝑗𝑗,1 𝑟𝑟� -K 0 �𝑅𝑅 𝑗𝑗,2 𝑟𝑟�� 𝑛𝑛 𝑗𝑗=1 + 𝐾𝐾 0 + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) (5) [---------------------------F-------------------------] [--------b---------]
where Rj,1 and Rj,2 are distances to the image wells (𝑅𝑅

𝑗𝑗,1 = � �𝑥𝑥 -𝑥𝑥 𝑗𝑗 � 2 + �𝑦𝑦 -𝑦𝑦 𝑗𝑗 � 2 and 𝑅𝑅 𝑗𝑗,2 = � �2𝑑𝑑 -�𝑥𝑥 -𝑥𝑥 𝑗𝑗 �� 2 + �𝑦𝑦 -𝑦𝑦 𝑗𝑗 � 2
). Eq.5 can be decomposed into two terms 'F', corresponding to the dam, and 'b' corresponding to the linear (i.e. natural) gradient.

The evaluation of strength coefficients γj is achieved by transposing the nil flow condition around the dam sides

( 𝜕𝜕ℎ(𝜕𝜕,𝜕𝜕) 𝜕𝜕𝜕𝜕 | 𝜕𝜕=±𝐿𝐿𝜕𝜕 2 ⁄ = 𝜕𝜕ℎ(𝜕𝜕,𝜕𝜕) 𝜕𝜕𝜕𝜕 | 𝜕𝜕=±𝐿𝐿𝜕𝜕 2 ⁄ = 0).
The following equations give the derivative forms of Eq.5 according to x and y.

According to x:

𝜕𝜕ℎ 𝜕𝜕𝜕𝜕 = ∑ 𝛾𝛾 𝑗𝑗 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 𝑟𝑟𝑒𝑒 𝛼𝛼�𝜕𝜕-𝜕𝜕 𝑗𝑗 � �K 0 �𝑅𝑅 𝑗𝑗,1 𝑟𝑟� -K 0 �𝑅𝑅 𝑗𝑗,2 𝑟𝑟� - �𝜕𝜕-𝜕𝜕 𝑗𝑗 � 𝑅𝑅 𝑗𝑗,1 K 1 �𝑅𝑅 𝑗𝑗,1 𝑟𝑟� - �2𝑑𝑑-�𝜕𝜕-𝜕𝜕 𝑗𝑗 �� 𝑅𝑅 𝑗𝑗.2 K 1 �𝑅𝑅 𝑗𝑗,2 𝑟𝑟�� 𝑛𝑛 𝑗𝑗=1 -𝑖𝑖 (6)
And according to y:

𝜕𝜕ℎ 𝜕𝜕𝜕𝜕 = ∑ 𝛾𝛾 𝑖𝑖 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 𝑟𝑟�𝑦𝑦 -𝑦𝑦 𝑗𝑗 �𝑒𝑒 𝛼𝛼�𝜕𝜕-𝜕𝜕 𝑗𝑗 � � K 1 �𝑅𝑅 𝑗𝑗,2 𝛼𝛼� 𝑅𝑅 𝑗𝑗,2 - K 1 �𝑅𝑅 𝑗𝑗,1 𝛼𝛼� 𝑅𝑅 𝑗𝑗,1 � 𝑛𝑛 𝑗𝑗=1 (7)
Where K1(x) is the modified Bessel function of the second kind of the first order.

Strength coefficients γj were evaluated on the dam sides at the closest distance (d') between the image wells and the dam sides (Fig. 2b). As the nil flow condition is computed in n locations, it results in a series of n linear equations that have to be solved simultaneously:

⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ 𝛾𝛾1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F11′+ 𝛾𝛾2 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F21′+ 𝛾𝛾3 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F31′+…. 𝛾𝛾n 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 Fn1′+b1′=0 𝛾𝛾1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F12′+ 𝛾𝛾2 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F22′+ 𝛾𝛾3 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F32′+…. 𝛾𝛾n 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 Fn2′+b2′=0 𝛾𝛾1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F13′+ 𝛾𝛾2 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F23′+ 𝛾𝛾3 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F33′+…. 𝛾𝛾n 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 Fan3′+b3′=0 … 𝛾𝛾1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F1j′+ 𝛾𝛾2 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F2j′+ 𝛾𝛾3 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F3j′+…. 𝛾𝛾n 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 Fni′+b j ′=0 … 𝛾𝛾1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F1n′+ 𝛾𝛾2 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F2n′+ 𝛾𝛾3 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 F3n′+…. 𝛾𝛾n 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 Fnn′+b n ′=0 ⎭ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎫ (8)
where superscript '′' denotes derivative terms of function 'F' and 'b' (Eq. 5), and subscripts the location where the calculation is performed and the number of the image well.

Equation 8 can be represented in matrix form: [𝐴𝐴] -1 ), and the hydraulic head is given by:

1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ F11′ F21′ F31′ … Fn1′ F12′ F22′ F32′ … Fn2′ F13′ F23′ F33′ … Fn3′ … . F1j′ F2j′ F3j′ … Fnj′ … F1n′ F2n′ F3n′ … Fnn′ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ γ 1 γ 2 γ 3 … 𝛾𝛾 𝑗𝑗 … γ 𝑛𝑛 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ -b 1 ′ -b 2 ′ -b 3 ′ … -b 𝑗𝑗 ′ … -b 𝑛𝑛 ′⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ (9) Equation 9 is of the form 1 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 [A]�γ 𝑗𝑗 �=[B],
ℎ(𝑥𝑥, 𝑦𝑦) = ∑ γ ȷ �𝑒𝑒 𝛼𝛼�𝜕𝜕-𝜕𝜕 𝑗𝑗 � �𝐾𝐾 0 �𝑅𝑅 𝑗𝑗,1 𝑟𝑟� -𝐾𝐾 0 �𝑅𝑅 𝑗𝑗,2 𝑟𝑟�� 𝑛𝑛 𝑖𝑖=1 + 𝐾𝐾 0 + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) (10)
Equation 10 shows that the steady-state hydraulic head induced by a fully penetrating underground dam does not depend on the hydraulic conductivity of the aquifer.

Hydraulic-head solution for an underground dam, a pumping or an injecting well, and an artificial recharge area

If other devices, like a pumping or injecting well and/or artificial recharge areas (MAR), are located near the underground dam (Fig. 2c), the same method can be used for evaluating the steady-state hydraulic head. Due to the linearity of the solution (Eqs. 2 and 4), the principle of superposition is used for adding such sources. It results that the terms corresponding to these other sources have to be added to the second part of Eq. 10 and thus to the known vector [B]. Therefore, [B], �γ ȷ ��, and thus hydraulic head, now also depend on hydraulic conductivity. The new system of linear equations obtained is also solved numerically.

For a well pumping (or injecting) at a rate Qp (in m 3 /s) located at xw, yw, the term to be added is identical to Eq. 3, except that x is replaced by x-xw, y by y-yw and

d by d-xw. If the well is pumping, Qp is <0; if it is injecting, Qp is >0.
For a rectangular recharge area (MAR), with uniform percolation rate R (m/s), located at xR, yR (centre of the recharge area) with lengths 2xL and 2yL (both in m) parallel to the axes of the coordinates systems (Fig. 2c), the hydraulic-head solution is the integral of Eq. 3 according to 2xL and 2yL [START_REF] Zlotnik | Estimating groundwater mounding in sloping aquifers for managed aquifer recharge[END_REF][START_REF] Dewandel | Artificial aquifer recharge and pumping: transient analytical solutions for hydraulic head and impact on streamflow rate based on the spatial superposition method[END_REF].

Therefore, the semi-analytical steady-state hydraulic-head solution combining an underground dam with a pumping or injecting well, and an artificial recharge area, is given by:

ℎ(𝑥𝑥, 𝑦𝑦) = ∑ γ ȷ �𝑒𝑒 𝛼𝛼�𝜕𝜕-𝜕𝜕 𝑗𝑗 � �K 0 �𝑅𝑅 𝑗𝑗,1 𝑟𝑟� -K 0 �𝑅𝑅 𝑗𝑗,2 𝑟𝑟�� 𝑛𝑛 𝑖𝑖=1 + 𝐾𝐾 0 + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) + 𝑄𝑄 𝑝𝑝 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 e 𝛼𝛼𝜕𝜕 𝑝𝑝 �K 0 �𝑟𝑟 𝑝𝑝 𝑟𝑟� -K 0 �𝑟𝑟 𝑝𝑝𝑛𝑛 𝑟𝑟�� + 𝑅𝑅 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 ∫ ∫ e 𝛼𝛼(𝜕𝜕-𝜕𝜕 𝑅𝑅 -𝜕𝜕′) [K 0 (𝑟𝑟′𝑟𝑟) -K 0 (𝑟𝑟 𝑛𝑛 ′𝑟𝑟)] +𝜕𝜕 𝐿𝐿 -𝜕𝜕 𝐿𝐿 +𝜕𝜕 𝐿𝐿 -𝜕𝜕 𝐿𝐿 𝑑𝑑𝑥𝑥′ 𝑑𝑑𝑦𝑦′ (11) with xp=x-xw; 𝑟𝑟 𝑝𝑝 = �(𝑥𝑥 -𝑥𝑥 𝑤𝑤 ) 2 + (𝑦𝑦 -𝑦𝑦 𝑤𝑤 ) 2 , 𝑟𝑟 𝑝𝑝𝑛𝑛 = �(2𝑑𝑑 -𝑥𝑥-𝑥𝑥 𝑤𝑤 ) 2 + (𝑦𝑦 -𝑦𝑦 𝑤𝑤 ) 2 , 𝑟𝑟′ = �(𝑥𝑥 -𝑥𝑥 𝑅𝑅 -𝑥𝑥 ′ ) 2 + (𝑦𝑦 -𝑦𝑦 𝑅𝑅 -𝑦𝑦 ′ ) 2 and 𝑟𝑟 𝑛𝑛 ′ = �(2𝑑𝑑 -𝑥𝑥 -𝑥𝑥 𝑅𝑅 -𝑥𝑥 ′ ) 2 + (𝑦𝑦 -𝑦𝑦 𝑅𝑅 -𝑦𝑦 ′ ) 2 .
The pumping or injecting well, as well as the recharge area, can be located up-or downstream from the underground dam. Examples of this solution are given later.

Numerical analysis of the proposed solutions

Accuracy of the semi-analytical model

The accuracy of the method used (MFS) depends on the distances between image wells (σ) and on the distance of these wells to the boundary of the underground dam (d', Fig. 2b). Examples here are given for the case of an underground dam alone (Eq.10), but results are similar if the solution includes a pumping or injecting well, or an artificial recharge area. Figure 3a gives examples of hydraulic gradient computation along the same underground dam for different d'/σ ratios. For ratios lower than 2 results are clearly not satisfactory and become acceptable for higher ratios. So the question is what ratio should be used to get acceptable results?

Figure 3b shows several sensitivity tests of the average hydraulic gradient calculation along the dam for different dam geometries (50≤Lx≤800, 2≤Ly≤20), different d' (0.1≤d'≤1.0) and number of image wells (100≤n≤4200). The hydraulic gradient is normalized according to the natural gradient (i.e., without a dam). The tests were also done with a different natural hydraulic gradient (i.e. no dam) and D0, without affecting the result (0.003≤i≤0.02 and 12≤D0≤42 m). The figure shows that the accuracy of the proposed semi-analytical solution, i.e. a nil hydraulic gradient along the dam, is clearly controlled by the number of image wells used, which is itself controlled by the dam geometry. Obviously, increasing n increases the result accuracy, but for reasonable computation times, from a few tens of seconds to a few minutes, the number of image wells can be optimized according to d'. Results show that for the ratio d'/σ ranging between 2.0 and 3.0, the semi-analytical solution gives satisfactory results.

Comparison with numerical modelling

The performance of the proposed solution has been compared with numerical modelling. The numerical modelling was developed with the MARTHE 7.4 © BRGM computer code [START_REF] Thiéry | Groundwater flow modelling in porous media using MARTHE[END_REF][START_REF] Thiéry | MARTHE calculation code -3D modelling of flows in hydro-systems -User manual for version 7.5[END_REF][START_REF] Thiéry | Modelling flow and heat transfer through unsaturated chalk -Validation with experimental data from the ground surface to the aquifer[END_REF]. MARTHE allows 2D or 3D modelling of flow and mass transfer in aquifer systems, including the effects of climatic and human influences. Groundwater flow is computed by a 3D finite-volume approach to solve the hydrodynamic equation based on Darcy's law and mass conservation. The grid geometry used in the numerical model (Fig. 4) is a 22.5×20 km rectangular grid (length on the x-axis: 22.5 km), with boundary conditions at constant head on the east and west sides, and no-flow boundary conditions on the north and south sides (length on the y-axis 20 km). The underground dam is located at the origin of the Cartesian coordinate system. In the x, y plane, cell size varies from 1×5 m at the origin up to1000×1000 m near the model boundaries. Vertically, the aquifer is a 42-m-thick layer subdivided into four 10-m-thick and one 2-m layers (i.e. D0=42 m) and is inclined (i=0.01). The dam is located at 2.5 km (d) from the eastern model limit (constant-head boundary). The aquifer is unconfined, the hydraulic head without dam is parallel to the aquifer substratum, and the model computes a steady-state hydraulic head, as in the analytical model. For numerical modelling purposes, hydraulic conductivity was set at 5.0×10 -4 m/s. However, its value does not affect the result.

Several rectangular underground dam geometries, Lx = 3, 5, 10 m and Ly = 15, 55, 105, 255, 405 and 505 m, were tested, but only some of them are presented below. However, comparisons are similar. Some of the results of the whole modelling are presented in the next section ( § 4.). Figures 5a,b and c compare the results of numerical modelling with analytically computed hydraulic-head values (Eq. 10). As shown on the figure, values computed with the proposed solution, for various 2D cross-sections along the x-axis, match the numerical modelling values.

The standardized root mean square error (RMSE D) was calculated for each numerical-analytical comparison of hydraulic-head values (see inserted tables on Fig. 5a,b&c). RMSE D values range between 0.14% and 0.81%, showing that the proposed semi-analytical solution (Eq.10) agree with the numerical modelling, a good indication of its accuracy. Note that Figure 5a corresponds to a case presented on Figure 3a (black curve with orange circles for d'/σ=2.5, n=2500).

Maximum hydraulic-head variation on the dam sides

The maximum hydraulic-head variation induced in the aquifer by the underground dam logically occurs in the middle of the dam, on its sides (i.e. h(-Lx/2,0) and h(Lx/2,0)). As the proposed semi-analytical solution (Eq.10) does not produce a clear mathematical expression of maxima, several tests tried establishing an empirical relationship between the dam geometry, the aquifer slope and thickness, and the maximum hydraulic-head variations up-and downstream of the dam. The distance to the constant-head boundary (d) was fixed at infinite (d=5000 m, in practice d>5Ly suffices); consequently, its influence is neglected. The tests were carried out only for a dam (no pumping, no artificial recharge area) to avoid overlying complex configurations.

Tests were done with 12≤D0≤42 m and 0.0001≤i≤ 0.02, 1≤Lx≤250 m and 3≤Ly≤1300 m, giving a Ly/Lx ratio ranging between 0.02 and 400. Differences in hydraulic head up-and downstream of the dam were normalized according to the natural hydraulic-head decrease that occurs at the dam location, i.e. (hup-hdown)/(iLx).

Figure 6a shows that the normalized hydraulic-head variation at the dam centre follows a typical characteristic curve as a function of the dam geometry ratio Ly/Lx. For Ly/Lx ratios over 5, this relation can be reasonably approximated by a linear relationship. For smaller ratios, the curve follows a power function starting with normalized hydraulic-head difference of 1 that corresponds to Ly/Lx →0 (absence of dam). As different aquifer thicknesses were used to construct this plot, it seems that its influence on the maximum hydraulic variations can be neglected (no correlation found). Similarly, a normalized maximum increase (hup-hnat.)/(iLx) and a maximum decrease (hnat.-hdown)/(iLx) of hydraulic head at the centre of the dam and on its sides, according to the natural hydraulic head (i.e. no dam, hnat.) were also plotted according to Ly/Lx ratios (Fig. 6b). The results show that differences in normalized hydraulic head are symmetrical and that they also follow near-linear relationships.

Results from the numerical modelling are plotted on graphs (Figs. 6a andb) and confirm the semi-analytical solution (Eq.10). Two field examples, corresponding to the scheme of the developed semi-analytical solution, are also presented. The Massillargues-Attuech subsurface dam is located about 40 km northwest of the city of Nîmes in France (Fig. 7). It features a 470 m long and 3 m wide underground dam (i.e. Ly/Lx=157), built in the 1990s to raise the groundwater level of the alluvial aquifer that became depleted due to gravel quarrying. The aquifer mainly consists of gravel with a hydraulic conductivity of about 5x10 -3 m/s and a thickness of 6 to 10 m, which is completely penetrated by the dam. The hydraulic gradient is about 3.5x10 -3 . The difference of about 1.7 m in hydraulic head, up-and downstream of the dam, was measured between the lake that is fully connected to the aquifer and an observation well located about 20 m downstream of the dam. The point corresponding to the difference in hydraulic head (normalized) falls almost on the theoretical curve (Fig. 6a). Using the defined empirical relationship, the difference in hydraulic head is evaluated at 1.67 m between both sides of the dam, a deviation of less than 2% according to the field measurements.

The second field example is given by [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF]. It corresponds to the basement of the El Prat de Llobregat railway station (Spain), where 29-m deep and impermeable diaphragm walls were used to make an enclosure, completely cutting an aquifer formed by sand and gravel deposits. Its thickness ranges from 12 to 20 m, with hydraulic conductivity ranging from 1 to 5x10 -4 m/s. The basement of the railway station can be assimilated to an underground dam that fully penetrates the aquifer. Its length and width are about 1800 m and 200 m, respectively (Ly/Lx=9). According to the given reference, the hydraulic gradient is 6.6x10 -4 , and the variation in hydraulic head from two observation wells, located on either side of the railway station, is about 1.25 to 1.3 m (1.275 m on average). This second field example agrees well with the theoretical curve (Fig. 6a), and the use of the defined empirical relationship gives a difference in hydraulic head of 1.31 m up-and downstream of the basement, thus a deviation of less than 3% according to field measurements.

Theoretical examples of underground dams with artificial recharge area and pumping well

This section presents several theoretical examples associating an underground dam and an artificial recharge area, as well as an underground dam with an artificial recharge area and a pumping well (Eq.10). The objective is to evaluate the effect of the hydraulic conductivity of the aquifer on the hydraulic head and to study how its distribution is impacted by the overall device.

For all simulations, the dam geometry is fixed at Lx= 5 m and Ly=500 m and the artificial recharge area is 2x10 4 m 2 (2xL=100 m and 2yL=200 m) with a uniform percolating flux (R) fixed at 5.56×10 -6 m/s. The natural hydraulic gradient (no dam, no pumping, no MAR) is i=0.005, D0=25 m, and d was set far away to overcome its influence (d=5 000 m). For the first example (Fig. 8), the dam is located at the eastern limit of the infiltrating area (see figure ). For the second one (Fig. 9), the pumping well is located in the middle between dam and recharge area, with a pumping rate fixed at 400 m 3 /h i.e. the same withdrawal rate as the injection one into the recharge area. For the purposes of this modelling, four hydraulic conductivity values were tested: 3×10 -3 , 5×10 -3 , 3×10 -2 and 10 -2 m/s, a range that generally represents aquifers formed of sand, gravel, or a mixture of the two [START_REF] Kruseman | Analysis and evaluation of pumping test data[END_REF].

Figure 8a shows the impact of an underground dam on the steady-state hydraulic head with an artificial recharge area; the black curve shows the steady-state hydraulic head with the dam only. Figure 8b shows the steady-state hydraulic head without a dam, showing that, with a dam alone, the hydraulic head rises by about 1.4 m upstream of the dam and decreases by the same amount downstream. With the addition of artificial recharge, the hydraulic head increases on both sides of the dam, up to +3.5 m upstream and even higher than the natural hydraulic head downstream for less permeable cases (K<10 -2 m/s). These increases are caused by the underground dam that upstream blocks the flux, but lets it pass around the sides, causing the downstream increase. Compared to the cases without a dam (Fig. 8b), the dam increases the hydraulic head induced by artificial recharge up to 1.6 m (K=3×10 - 3 ), locally multiplying by a factor of about two its effect. However, where the hydraulic conductivity is high (K>5×10 -3 m/s), the groundwater 'mound' induced by recharge area flattens out (Fig. 8b), and the effect induced by the dam becomes largely dominant (Fig. 8a).

Where a pumping well is present, it logically decreases the hydraulic head of the whole system (Figs. 9a &b). The dam increases the hydraulic head at the level of the pumping well by up to several tens of centimetres, depending on the case (Fig. 9a), and accentuates the 'performance' of the recharge area. Downstream, the hydraulic head computed for the different scenarios of hydraulic conductivity is very close to that with a dam alone, due to the same pumping rate as the injection rate into the recharge zone (Fig. 9a).

Discussion

Above, steady-state semi-analytical solutions were developed for evaluating the distribution of hydraulic head, caused by an underground dam, with or without a pumping or injection well or an artificial recharge area (Equations 10 and 11). They assumed that the rise or decline of the water table remains small in comparison to the initial saturated thickness of the aquifer (i.e. (𝐾𝐾 -𝐾𝐾 0 ) 𝐾𝐾 0 ⁄ is small). This method of linearization, commonly used for solving groundwater hydraulic problems (e.g. [START_REF] Theis | The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage[END_REF], was chosen to solve the standard flow equation in steady-state condition (Eq. 1). According to textbooks such as de [START_REF] Marsily | Quantitative Hydrology for Engineers[END_REF], it can be assumed that if (𝐾𝐾 -𝐾𝐾 0 ) 𝐾𝐾 0 ⁄ is <10%, the results can be considered as exact. Based on various theoretical and numerical modelling works for the groundwater mounding induced by artificial recharge in sloping aquifers, though extendable to other situations, [START_REF] Zlotnik | Estimating groundwater mounding in sloping aquifers for managed aquifer recharge[END_REF] proposed the following guidelines for this type of linearization. The relative differences are around a few percent if (D-D0)/D0≤0.5 and on about 10% if 0.5<(D-D0)/D0≤1. For diagnostic purposes, they also suggested that the results may be acceptable even at (D-D0)/D0>1. For example, the numerical modelling example given on Fig. 5a has a ratio of about 5.6% and comparison with the semi-analytical solution (Eq.10) is less than 1% (RMSE_D on curves is 0.7 to 0.8% and at the point of maximum difference in hydraulic head the difference is <0.6%). This again justifies the use of this type of linearization.

Concerning the maximum rise or depletion in hydraulic head caused by an underground dam (i.e. in the middle of the dam, on its sides), several tests with the semi-analytical solution (Eq.10) established an empirical relationship, depending only on the dam geometry (Ly/Lx ratio) and the natural hydraulic gradient, as was further confirmed by numerical modelling (Fig. 6). This empirical relationship assumed that the dam is rectangular and that its greatest dimension is perpendicular to the gradient; its validity is given for Ly/Lx ratios >5 and i≤0.02. These maxima, as mentioned earlier for the hydraulic-head distribution (Eq.10), do not depend on the hydraulic conductivity of the aquifer, since the flow condition without dam is only a function of the hydraulic gradient (i). This point was also mentioned by [START_REF] Pujades | Barrier effect of underground structures on aquifers[END_REF], who derived semi-empirical equations for quantifying the long-term effects of different types of hydraulic barriers (semi-permeable, partially or fully penetrating the aquifer) near the underground dam. The semi-analytical solutions proposed here complement Pujades et al.'s work, by evaluating the spatial hydraulic-head distribution caused by a dam. The proposed graph and the associated empirical relationship (Fig. 6) allows directly assessing the maximum hydraulic head variation on either side of the dam, knowing only the dam geometry and the regional hydraulic gradient (≤0.02). Its validity was tested with numerical modelling and at two field sites, obtaining an excellent agreement between the semi-analytical solution and, the modelled and observed variations. Also note that for the field examples, the validity of the proposed linearization of flow equation was not necessarily respected; indeed (D-D0)/D0 is 0.11 for the Massillargues-Attuech case and 0.04 for the El Prat de Llobregat railway station, thus reinforcing the linearization assumption used here.

As shown by the semi-analytical solution, the hydraulic head rise and decline may propagate far up-and downstream of the dam (see previous figures and Fig. 10 for a 3D view). According to the tests performed, a significant difference of a few centimetres from the natural hydraulic head occurs at about 4 to 5 times the length of the dam upstream and perpendicular to the dam, its effect downstream being felt up to 4 times the dam length. Parallel to the dam, the impact is much less, up to twice the dam length. This rise in hydraulic head may benefit upstream pumping wells, or it may reduce seawater intrusion, but it may also flood basements, or create new wetlands like the Gardon lake near Massillargues-Attuech (Fig. 7a). Conversely, downstream of the dam such depletion can decrease the yields of pumping wells -or even cause them to be abandoned, wetlands may fall dry and seawater can intrude coastal areas. It is thus necessary to carry out hydrogeological and pre-design studies for evaluating the environmental, economic and social impacts of all proposed underground dams (e.g. [START_REF] Ishida | Construction of subsurface dams and these impacts on the Environment[END_REF], Yihdego 2016[START_REF] Attard | Review: Impact of underground structures on the flow of urban groundwater[END_REF], Li et al 2018, Sun et al 2019…).

Several examples have shown the effect of hydraulic conductivity of the aquifer, when pumping, artificial recharge and an underground dam are coupled (Eq.11; Figs. 8 and9). Logically, a dam enhances the performance of artificial recharge by blocking flow. To some extent, the coupling with artificial recharge reduces the downstream depletion created by the dam, compared to a situation with a dam alone. With pumping wells, the dam prevents flow from propagating downstream, but they may accentuate the downstream depletion created by a dam. In the case of highly permeable aquifers (K>5×10 -3 m/s, Fig. 8a), the increase in hydraulic head created by the dam is much higher than that created by artificial recharge alone, despite a high injected flux (400 m 3 /h in the given examples). This suggests that, under such conditions, underground dams are more efficient for increasing the hydraulic head, and that, as mentioned by previous authors (BCEOM 1978, Choi and[START_REF] Choi | Effects of hybrid-type artificial groundwater recharge and underground barrier in a small basin[END_REF], such structures may be needed to prevent artificially recharged water from escaping. Thus, depending on the objective of an underground dam (store water for water supply, restore former water table conditions after quarrying activities or to mitigate seawater intrusion) such structures can be coupled with other devices, such as MAR, to improve their effectiveness.

The developed solutions assumed steady-state condition, and are thus suitable for evaluating long-term effects. However, during the building of an underground dam and after, it is of primary interest to model the transient behaviour of the hydraulic head. This allows, at least, to prevent a water-level rise or depletion near the dam and beyond, and to refine the design of the dam so as not to exceed a certain water level or to evaluate the time required to reach a quasi-steady state condition (El-Hames 2012). In addition, transient behaviour is needed for understanding how the resulting rise and/or depletion evolved under dry and/or wet conditions, and to assess, particularly in the dry season, how long the hydraulic head rise created by the dam will be effective for water supply wells, for example. Another point, not been explored here but which the authors want to emphasize, is the increase of flux at the dam edges. Figure 3a shows that the regional gradient near the dam walls is increased by a factor 4.5. Other simulations with a long dam (up to 1300 m, not shown here), increase the gradient by a factor of about 10. This can cause geotechnical issues to prevent edges from collapse or the appearance of voids in the ground.

Conclusions

This study proposes a steady-state semi-analytical solution for assessing the long-term spatial distribution of the hydraulic head caused by an underground dam whose extension concerns a limited part of the aquifer width Eq.10). The aquifer is assumed homogeneous, slanted and semi-finite. The dam fully penetrates the aquifer, is impermeable and is assumed to be rectangular with one of its length perpendicular to regional groundwater flow. The developed solution is based on the method of fundamental solutions and this work highlights that this tehcnique is a useful analytical technique to overcome the difficulty of finding an analytical solution in complex systems. However, the accuracy of the method depends upon the number of image wells used and their distance to the dam sides. Sensitivity tests show that the results are reliable for d'/σ ratios ranging between 2.0 and 3.0. The semi-analytical solution satisfactorily compares with numerical modelling. It is also demonstrated that the hydraulic head does not depend on aquifer hydraulic conductivity.

Dimensionless graphs relating the dam geometry to maximum hydraulic-head variations up-and downstream of the dam are given. An empirical linear relationship of the maximum difference in hydraulic head between the two sides of the underground dam, which depends only on the dam geometry (Ly/Lx ratio) and the regional natural hydraulic gradient (i) is established. This behaviour agrees with the results of previous studies and with the results of numerical modelling and at two field sites. The proposed relationship is valid for Ly/Lx ratios >5 and i≤0.02, and can be used for pre-design purposes.

In addition, the spatial extent of the rise and fall of the hydraulic head created by an underground dam was examined. Hydraulic-head changes propagate far upstream and downstream of the dam, about 4 to 5 times the length of the dam.

A Semi-analytical solution is also given when an artificial recharge area (MAR) and/or a pumping (or injecting) well are located near the underground dam (Eq.11). For these cases, hydraulic head further depends on aquifer hydraulic conductivity. It is shown that the dam enhances the performance of artificial recharge by reducing downgradient flow. Interestingly, in the case of highly permeable aquifers, the increase in hydraulic head created by the underground dam may be much higher than that created by artificial recharge alone, despite high injected flux. In these cases of highly permeable aquifers at least, combining a MAR device with an underground dam appears to be an attractive combination for preventing the escape of artificially recharged water.

Respecting the mathematical assumptions of the model, the proposed solutions can be used to assess the long-term spatial effects on groundwater to prevent possible environmental, economic and social impacts. They are intended to guide the design of such structures, especially to quickly test various configurations, which with conventional numerical modelling techniques can be time consuming. Future works should consider the transient behaviour of the hydraulic head not only during the building of the underground dam, but to evaluate how long is needed to reach quasi-steady state, and to assess during the dry season over how long the rise created by the dam will be effective for water supply wells and other abstraction needs. Other aquifer and underground dam geometries could also be tested.

Figure captions

  saturated thickness (in m) and ℎ(𝑥𝑥, 𝑦𝑦) = 𝐾𝐾(𝑥𝑥, 𝑦𝑦) + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) the hydraulic head (in m). The boundary conditions are: ℎ(𝑑𝑑, 𝑦𝑦) = 𝐾𝐾 0 , ℎ(-∞, 𝑦𝑦) = 𝐾𝐾 𝑜𝑜 + 𝑖𝑖(𝑑𝑑 -𝑥𝑥) and along the dam: is the saturated thickness without dam and i=tanθ.

  where [A] represents the equations matrix, [γj] is the unknown variable vector and [B] the known vector. With �𝛾𝛾 � 𝑗𝑗 � = �𝛾𝛾 𝑗𝑗 � 2𝜋𝜋𝐷𝐷 0 𝐾𝐾 , the dimensionless strength coefficients were evaluated numerically using matrix calculation (�𝛾𝛾 𝚥𝚥 � � = [𝐵𝐵].

Fig. 1 .

 1 Fig. 1. Conceptual models of the theoretical examples presented. a) An underground dam alone and b) an underground dam with a pumping well and a managed aquifer recharge device (MAR).

Fig. 2 .

 2 Fig.2. Sketch of the conceptual model of an underground dam. a) Section view (dam alone). b) Plan view (dam alone) with the displays of image wells; the equation on the left side corresponds to Eq. 5 and the one on the right corresponds to ones of Eq. 8. c) Plan view with an underground dam, a pumping or injecting well, and an artificial recharge area.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Sensitivity tests performed with Eq. 10 for an underground dam alone. a) hydraulic-gradient computation along the underground dam for different d'/σ ratios. Lx=3 m, Ly= 505 m, i= 0.01, D0=42 m, d'= 0.8m and 500≤n≤3000. b) d'/σ vs. normalized hydraulic-head gradient calculation along the dam (average of 1 𝑖𝑖

Fig. 5 .

 5 Fig. 5. Comparison of the analytical solution with numerical modelling (dam alone; Eq.10). a) Lx: 3 m, Ly: 505 m, b) Lx:5 m, Ly: 255 m, and c) Lx:3 m, Ly: 105 m. For all cases, D0=42 m, i: 0.01 and d: 2500 m. RMSE D: standardized residual mean square error.

Fig. 6 .

 6 Fig.6. Underground dam geometry ratio Ly/Lx vs. a) normalized hydraulic-head variation at the dam centre, i.e. h(-Lx/2, 0) and h(Lx/2, 0). The proposed equation is given for Ly/Lx>5. Also shown: results of the numerical modelling and the two field examples. b) Normalized maximum increase (hup-hnat.)/(iLx) and maximum decrease (hnat.-hdown)/(iLx) of hydraulic head at the dam centre and on its sides. On graph: h_dam=hup or hdown, hnat.: hydraulic head without dam (natural flow condition). Hydraulic-heads were computed with Eq.10.

Fig. 7 .

 7 Fig. 7. Massillarges-Attuech underground dam (470 m long and 3 m wide) in France. a) Location and water-table map in October 2022 (metres above sea level) on a 25x25 m cell grid; line contours: 0.5 m. b) Cross sectionshowing topography and water-table (m.a.s.l.); the underground dam creates, on its sides, a maximum hydraulichead variation of about 1.7 m (∆h).

Fig. 8 .

 8 Fig.8. Theoretical examples of an underground dam with a recharge area for different aquifer hydraulic conductivities (Eq.11 with Qp=0). Lx=5 m and Ly=500 m, 2xL=100 m and 2yL=200 m, R=5.56x10 -6 m/s, i=0.005, D0=25 m and d=5000 m. K=3×10 -3 , 5×10 -3 , 3×10 -2 and 10 -2 m/s. a) With the underground dam. b) Without dam, only the recharge area. Natural hydraulic head: no dam, no pumping, no recharge area.

Fig. 9 .

 9 Fig.9. Theoretical examples of an underground dam with recharge area and pumping well for different aquifer hydraulic conductivities (Eq.11). Lx=5 m and Ly=500 m, 2xL=100 m and 2yL=200 m, R=5.56×10 -6 m/s, i=0.005, D0=25 m and d=5000 m. K=3×10 -3 , 5x10 -3 , 3×10 -2 and 10 -2 m/s, Qp=400 m 3 /h. The pumping well is located 150 m from the dam and the MAR. a) With the underground dam. b) Without dam, only the recharge area and the pumping well. Natural hydraulic head: no dam, no pumping, no recharge area.

Fig. 10 .

 10 Fig. 10. 3D (half-space) graph of hydraulic-head disturbance created by an underground dam (Eq.10). Lx= 5 m and Ly=500 m, i=0.005, D0=25 m and d=5000 m. Black contour lines: hydraulic head and, blue contour lines (∆h): difference in hydraulic head compared to condition without dam. The arrows indicate the direction of groundwater flow.

  Fig.10

Acknowledgements

We warmly thank Mr. Marcel Masson who participated in constructing the experimental Massillargues-Attuech underground dam from 1989 to 1994, and who documented the design of the structure in 2018. The two anonymous Journal referees are thanked for their useful remarks and comments that improved the quality of the paper. The authors are grateful to Dr. H.M. Kluijver for revising the final version of the manuscript.

Funding information

This study was conducted by BRGM, and was funded by the internal BRGM CODHYDRO project and the 'Attuech Underground Dam' project funded by Etablissement Public Territorial de Bassin (EPTB) des Gardons, the French Water Agency Rhône Mediterranean & Corsica (AERMC), Occitanie Region and BRGM.

Conflict of interest statement

No, conflict of interest