Gabriel Marchesan Almeida

Remi Busseuil

Luciano Ost

Florent Bruguier

Gilles Sassatelli

Pascal Benoit

Lionel Torres

Michel Robert

est destinée au dépôt

PI and PID Regulation Approaches for Performance-Constrained Adaptive Multiprocessor

System-on-Chip Gabriel Marchesan Almeida, Student Member, IEEE, Rémi Busseuil, Luciano Ost, Florent Bruguier, Gilles Sassatelli, Member, IEEE, Pascal Benoit, Member, IEEE, Lionel Torres, and Michel Robert Abstract-Adaptive multiprocessor systems are appearing as a promising solution for dealing with complex and unpredictable scenarios. Given the large variety of possible use cases that these platforms must support and the resulting workload variability, offline approaches are no longer sufficient because they do not allow coping with time changing workloads. This letter presents a novel approach based on the utilization of PI and PID controllers, widely used in control automation, for optimizing resources utilization in Multiprocessor System-on-Chip (MPSoC). Several architecture characteristics such as response time during frequency changing, noise and perturbations are modeled and validated in a high-level model and results are compared to information obtained on a homogeneous MPSoC platform prototype. Power and energy consumption figures are discussed and two controllers are proposed: 1) PI-; and 2) PID-based controllers. Results show the system capability of adapting under disturbing conditions while ensuring application performance constraints and reducing energy consumption. Index Terms-Closed loop system, linear feedback control systems, multiprocessing systems, multitasking, system-on-a-chip.

I. INTRODUCTION

A DAPTATION is increasingly regarded as a promising technique to further optimizing a number of relevant parameters in future embedded systems, such as power consumption, quality of service, or even reliability. In multiprocessor systems, adaptation is often realized through dynamic task mapping, load balancing. Such techniques increase system efficiency through better use of processing and communication resources. These system-level adaptation strategies often rely on a centralized controller that performs the decision-making, and applies resulting decisions with typical latencies in the order of milliseconds. At lower-level, there also exist opportunities such as dynamic voltage and frequency scaling (DVFS). However, because of the distributed nature of the targeted systems, DVFS is either: 1) decided at design-time on a scenario-driven basis (frequency associated to a given system configuration); or 2) performed at run-time by means of distributed decision-making strategies as centralized approaches would incur a prohibitive latency (local monitoring and 2-way communication with the centralized unit).

DVFS techniques have been used in portable embedded systems in order to reduce the energy consumption and temperature of such systems. Most existing DVFS are defined at design time, where they are typically based on predefined profiling analysis (e.g., application data) that attempts at defining an optimal DVFS configuration [START_REF] Magklis | Profile-based dynamic voltage and frequency scaling for a multiple clock domain microprocessor[END_REF], [START_REF] Xie | Compile-time dynamic voltage scaling settings: Opportunities and limits[END_REF]. Puschini et al. [START_REF] Puschini | A game-theoretic approach for run-time distributed optimization on MP-SoC[END_REF] propose a scalable multiobjective approach based on game theory, which adjusts at run-time the frequency of each PE while reducing tile temperature and maintaining the synchronization between application tasks. Due to the dynamic variations in the workload of these systems and its impact on energy consumption, other adaptation techniques such as proportional-integral-derivative (PID)-based control have been used to dynamically scale the voltage and the frequency of processors [START_REF] Wu | Formal online methods for voltage/frequency control in multiple clock domain microprocessors[END_REF], [START_REF] Zhu | Feedback edf scheduling exploiting hardwareassisted asynchronous dynamic voltage scaling[END_REF], and recently, of networks-on-chip (NoCs) [START_REF] Ogras | Variation-adaptive feedback control for networks-on-chip with multiple clock domains[END_REF], [START_REF] Sharifi | Feedback control for providing qos in NoC based multicores[END_REF]. These techniques differ in terms of adopted control parameters (e.g., task deadline, temperature) and response times (period necessary to stabilize a new voltage/frequency). In [START_REF] Zhu | Feedback edf scheduling exploiting hardwareassisted asynchronous dynamic voltage scaling[END_REF], authors propose an analytic (queuedomain network) approach that explores the use of a PID controller, which adjusts the multiple clock domain (MCD) processor frequency according to the workload change. In turn, Zhu et al. [START_REF] Zhu | Feedback edf scheduling exploiting hardwareassisted asynchronous dynamic voltage scaling[END_REF] combine a PID-based controller with a DVS scheduler in order to make the system adaptive to varying workload, while reducing the power consumption.

Ogras et al. [START_REF] Ogras | Variation-adaptive feedback control for networks-on-chip with multiple clock domains[END_REF] propose an adaptive control technique for a multiple clock domain NoC, which considers the dynamic workload variation aiming at decreasing the power consumption. Besides, the proposed control technique ensures the operation speed and the frequency limits of each island (clock domain, defined according to a given temperature). The effectiveness of the proposed adaptive control technique was evaluated considering different scenarios (e.g., comparison with a PID controller) for a MPEG-2 encoder. In [START_REF] Sharifi | Feedback control for providing qos in NoC based multicores[END_REF], a PID controller is used to provide a predetermined throughput for multiple voltage-frequency/voltage-clock NoC architectures. The proposed PID-based controller sets the voltage and frequency of each NoC island by reserving virtual channels weights (primary control parameter) in order to provide the necessary throughput for different application communications, while saving energy. In [START_REF] Almeida | Predictive dynamic frequency scaling for multi-processor systems-on-chip)[END_REF], we have proposed a PID-based approach for ensuring application performance requirements. However, only PI components were used and no consideration about power and energy consumption is taken into account. The present letter presents the following contributions: 1) power and energy consumption are considered when tuning processor frequency; and 2) a PI-only controller is proposed and compared to a PID-controller;

The work presented in this letter relies on an open-scale MPSoC system that has been proposed in [START_REF] Almeida | An adaptive message passing MPSoC framework[END_REF]. It aims at devising a smart distributed frequency scaling strategy that makes it possible to meet real-time application requirements in the presence of perturbations originated from various sources among which the application-related time-changing workloads and especially the continuous adaptation process that the system undergoes, such as task migrations. This letter is organized as follows. Section II illustrates the proposed strategy while Section III discusses experimental results and conclusions are drawn in Section IV.

II. PROPOSED STRATEGY

Fig. 1 gives an overview of the proposed approach. As it can be observed, one PID controller is devoted to each task in the system that must ensure soft-real time constraints. In this example, there is one task per network processing unit (NPU), so one PID controller for each processor is required. In the case where multiple tasks are executed in the same NPU, a system with multiple PID controllers in the same NPU could be proposed. In this case, a simple heuristic would choose the highest frequency calculated by the different PID controllers in order to avoid deadline misses. The PID controller is implemented as a service in the operating system and it represents an overhead of less than 1% in terms of total occupied memory footprint.

Monitoring information such as task throughput is fed into the PID controller module that will match the actual throughput with the desired throughput (setpoint) and will then calculate an error value . This value is used for calculating and components and as result the PID controller will indicate a frequency such as to reach a given setpoint according to reactiveness factor initially set. It is important to observe that the PID management and calculation is performed at run-time. In cases where actual throughput is lower than the setpoint, the PID controller will indicate a frequency superior to the actual one in order to reach the setpoint and to satisfy application performance constraints. On the other hand, when the actual throughput is much higher than the expected one, the controller will sign to a lower frequency in order to save platform resources and consequently energy. The proposed strategy consists in deciding controller components on a task basis. To this purpose, a SystemC/TLM-based MPSoC simulation is executed in order to obtain the step-response of the process (see right box in Fig. 1) which corresponds to the NPU executing the task it hosts. In this scenario, processor frequency is changed from 55 MHz to 1005 MHz and task throughput is monitored. The system is linear, that means the system's behavior remains the same under such frequency change. Based on the high-level model, a number of different configurations of controllers can be explored. Each one exhibits different features such as speed, overshoot and static error. Once the process is modeled, PID components are fine tuned by using Simulink and the values of and are fed as input to the MPSoC platform. Fig. 2 presents the cycle-accurate simulation results showing both high-level model and platform behaviors. It is important to observe a very close match between the two models. In this scenario a PI-based model was used and components have been set to and . Fig. 3 illustrates the network processing unit represented in a high-level model. The system is mainly composed of:

1) process: represented by the DISCRETE FILTER in Fig. 3. It characterizes the system behavior based on the step-response process previously described; 2) noise generator: it represents system throughput variations of a given application; Fig. 3. High-level model of a network processing unit with PID controller, system characterization, and perturbation generator.

TABLE I REPRESENTATION OF THREE DIFFERENT PERTURBATION SCENARIOS WITH THEIR RESPECTIVE IMPACT

3) nerturbation generator: it reproduces the system behavior in the occurrence of perturbation in the system such as task migration, application's workload changes, etc; 4) PID controller: represented by the DISCRETE PID CON- TROLLER in Fig. 3. It symbolizes both PI and PID controllers used for tuning processor frequency according to application constraints. In the perturbation scenario, an application task is migrated from a different processor generating a perturbation in the system. Once our operating system is implemented in a thread-based way where CPU processing resources are shared among applications, whenever there is a new task running on it, the task will consume CPU power and then the application performance of the tasks previously running on it will decrease. The process modeling is composed of three steps: 1) perturbation scenario is implemented and validated in a real platform and performance figures are obtained; 2) based on this information, platform behaviors are modeled in a high-level model and perturbation amplitude in terms of performance are considered; and the last step consists in 3) matching both high-level model results in order to calibrate and fine adjust high-level model components.

III. CASE STUDY AND RESULTS

The following scenarios are based on a MJPEG video decoder application composed of five tasks running on a processor and a different application (a synthetic application implemented in such way that it requires a considerable amount of processing power) is migrated to the processor, disturbing the MJPEG application performance. Three different perturbation scenarios ( and ) are created based on different perturbation factors defined respectively as follows: 1) low; 2) medium; and 3) high, represented in Table I. . Application performance is analyzed for two approaches: 1) no adaptation is done and processor frequency keeps constant; and 2) processor frequency is tuned at real-time according to PID controller output values. The perturbation is configured as follows: 1) perturbation start: 300 ms; 2) perturbation period: 700 ms; and 3) application performance impact factor: 58%.

It is possible to observe that when there is no PI controller, the task throughput is reduced to approximately 440 KB/s during the perturbation period. If we consider soft-real time applications where performance constraints must be ensured, this scenario could imply on, for instance, freezing a video during the decoding process or having frame losses during the perturbation period. By adopting our proposed strategy it is possible to notice that the PI controller takes only 100 ms to react and reach the initial setpoint, that represents the minimal performance requirements of a given application. When perturbation is over, we can clearly observe that task throughput increases due to the fact that the processor is running in a high frequency in order to compensate the previous existing perturbation. Thus, the PI controller starts decreasing the processor frequency in order to save platform resources, getting back close to the initial setpoint.

In Table II, the application performance considering three different perturbation scenarios is presented. For each scenario we have analyzed three different configurations: 1) no PI-/PID- based controller is used and processor frequency keeps constant; 2) PI-and; 3) PID-based controllers are implemented. We can observe that in scenario , where perturbation impact factor is higher, the average task throughput when using PI-/PID-based controllers is of 1,000 KB/s, resulting on getting the desired performance, represented by the setpoint. When using the approach that does not apply PI-/PID-based controller, we can observe that the average throughput is of 1,675 KB/s, representing a waste of platform resources. This is due to the fact that it is necessary to configure the processor to run in a higher frequency to be capable of dealing with scenarios under perturbation conditions. The metric used for choosing the processor frequency consists in selecting the lowest possible frequency where processor can guarantee soft real-time application requirements under perturbation conditions.

In order to obtain the network processing unit power consumption as precisely as possible, we have synthesized the processor on a 65 nm CMOS technology. The design kit is provided by STMicroelectronics and two libraries are used: high voltage threshold (HVT) is used on the critical path of the circuit and low VT (LVT) is used on the remaining of the circuit. The power consumption measurements were obtained from an execution of an embedded operating system's boot loader on a processor running at 50 MHz. Note that, the memory consumption is not considered due to the fact that its operating frequency does not change with the processor frequency. The network processing unit power consumption is presented in Table III.

Table IV presents the power and energy consumption for three different perturbation scenarios. We can clearly see that power and energy consumption are significantly reduced when using PI-and PID-based controllers. In , the energy consumption is reduced by 32% while power consumption is reduced by 42% in . It is important to observe that the energy/power consumption saving can be much higher in scenarios with longer perturbation periods and where impact factor of perturbations over applications is higher. We emphasize that by using our proposed approach it is possible not only complying with application performance constraints and saving energy/power but also it is possible to reduce waste of platform resources in order to reach an optimal solution.

IV. CONCLUSION

This letter has presented a novel approach based on the utilization of PID controllers for energy consumption reduction while ensuring application performance constraints. The strategy does not rely only on preventing application deadline-misses but also attempts to save energy, once the processor frequency is adjusted at real-time according to application requirements. We claim that our approach can be easily integrated to linear systems and, as result, platform resources utilization can be optimized. Moreover, by using PID controllers it is possible to save up to 42% in terms of power by tuning processor frequency according to application needs. The letter also presents: 1) PI-and; 2) PID-based controllers. As can be observed, PID-controllers are intended to react faster under disturbing conditions when compared to PI-only controllers, however its power consumption is higher. The good choice of which controller to use in multiprocessor system-on-chip platforms will be a trade-off between power consumption and the desired system's reactiveness.

Fig. 1 .

 1 Fig. 1. Architectural overview of the MPSoC platform with PID controllers.

Fig. 2 .

 2 Fig. 2. System behavior-simulink versus MPSoC platform.

Fig. 4 .

 4 Fig. 4. Perturbation representation-PI versus no PI.

Fig. 4

 4 Fig.4depicts the system response under perturbation condition . Application performance is analyzed for two approaches: 1) no adaptation is done and processor frequency keeps constant; and 2) processor frequency is tuned at real-time according to PID controller output values. The perturbation is configured as follows: 1) perturbation start: 300 ms; 2) perturbation period: 700 ms; and 3) application performance impact factor: 58%.It is possible to observe that when there is no PI controller, the task throughput is reduced to approximately 440 KB/s during the perturbation period. If we consider soft-real time applications where performance constraints must be ensured, this scenario could imply on, for instance, freezing a video during the decoding process or having frame losses during the perturbation period. By adopting our proposed strategy it is possible to notice that the PI controller takes only 100 ms to react and reach the initial setpoint, that represents the minimal performance requirements of a given application. When perturbation is over, we can clearly observe that task throughput increases due to the fact that the processor is running in a high frequency in order to compensate the previous existing perturbation. Thus, the PI controller starts decreasing the processor frequency in order to save platform resources, getting back close to the initial setpoint.In TableII, the application performance considering three different perturbation scenarios is presented. For each scenario we have analyzed three different configurations: 1) no PI-/PID-

TABLE II APPLICATION

 II PERFORMANCE FOR THREE DIFFERENT SCENARIOS TABLE III NETWORK PROCESSING UNIT (NPU) POWER CONSUMPTION USING STMICROELECTRONICS DESIGN KIT FOR 65 NM TECHNOLOGY

TABLE IV NPU

 IV POWER AND ENERGY CONSUMPTION REPRESENTATION FOR THREE DIFFERENT SCENARIOS USING 65 NM TECHNOLOGY

The authors are with Laboratory of Informatics, Robotics, and Microelectronics of Montpellier (LIRMM), University of Montpellier II, Montpellier 34095, France (e-mail: gmalmeida@ieee.org; remi.busseuil@lirmm.fr; ost@lirmm.fr; florent.bruguier@lirmm.fr; gilles.sassatelli@lirmm.fr; pascal.benoit@lirmm.fr; lionel.torres@lirmm.fr; michel.robert@lirmm.fr).


Amine Dehbaoui
email: dehbaoui@lirmm.fr

Thomas Ordas

Victor Lomné

Philippe Maurine

Lionel Torres

Michel Robert

Incoherence Analysis and its Application to Time Domain EM Analysis of Secure Circuits

Since several years, Electromagnetic (EM) radiations of ICs are considered as source of noise and interferences from an EM Compatibility (EMC) point-of-view. Electromagnetic interferences (EMI) are unwanted disturbances that affect digital circuits due to electromagnetic conduction or electromagnetic radiation emitted from an internal or external source. Today electromagnetic radiations of ICs also represent vulnerabilities for hardware security modules like smartcards. Recently, from a hardware security point of view, researchers, industrials and governmental agencies are focusing with strong interest on these radiations. This paper aims at introducing a pragmatic technique WGMSI (Weighted Global Magnitude Squared Incoherence) allowing localizing hardware security modules within the whole EM noise generated by the environment or other modules of the circuits.

I. INTRODUCTION

In the last century, modern cryptology has mainly focused on defining cryptosystems resistant against theoretical attacks. However, with the increasing use of secure embedded systems like smartcards, researchers focused on exploiting the physical syndromes leaking from secure devices during a cryptographic operation to disclose the key. As a result, a new kind of attack called Side-Channel Attack (SCA) has appeared. Among the known attacks, some exploit the timing behavior of Integrated Circuits (IC) [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF], while others exploit the global power consumed by IC such as the well known Differential Power Analysis (DPA) [START_REF] Kocher | Differential Power Analysis[END_REF]. Recently, the Electro-Magnetic (EM) emanations of embedded systems have been identified as a major threat [START_REF] Gandolfi | Electromagnetic analysis: Concrete results[END_REF] [START_REF] Quisquater | Electromagnetic analysis (ema): Measures and countermeasures for smart cards[END_REF]. The efficiency of the EM channel is mainly due to the inner properties of EM emissions. Their ability to propagate through different materials is the most interesting one since it allows an attacker targeting the bounded hardware area to integrate the cryptographic algorithm under attack or part of it. This is all the more interesting since it also allows getting round global hardware countermeasures against power analyses such as the use of detached power supplies [START_REF] Shamir | Protecting smart cards from passive power analysis with detached power supplies[END_REF] by focusing the analysis on reduced die areas. However, this requires the use of small magnetic sensors to localize the leaking spots and thus, implies a quadratic increase (with the square of the ratio between the length of the measured chip and the size of the sensor) of the number of points to be attacked using Differential Electro-Magnetic Analyses (DEMA) [3][4]. Thus, according to the Firstly, it requires only few EM measurements to be efficiently applied. Secondly, it remains efficient (a) even in presence of data-independent EM emanations such as the ones generated by the clock tree or any 'always on' analogue blocks but also (b) in presence fully data-independent parasitic emissions such as noise. Thirdly, it allows finding positions where DEMA might be successful with a reduced set of EM traces. Finally, as last advantage, this non invasive and contactless technique can be applied with success, as demonstrated in Section 3, even if the circuit under attack is encapsulated. The remainder of this paper is organized as follows. Section 2 introduces theoretical explanations about the proposed localization technique called W eighted Global M agnitude Squared Incoherence (W GM SI)technique and how to couple W GM SI technique with EM near field scanning systems. Section 3 presents some concrete results related to the application of W GM SI localization technique to different mappings on FPGA of a same design. This section also demonstrates the efficiency of W GM SI cartography to lead a reverse engineering of the IC. Finally, conclusions are drawn in Section 4.

II. GLOBAL MAGNITUDE SQUARED INCOHERENCE

A. Problem definition

Let us consider, for simplicity, that there are only two local sources of EM emissions within the chip: the source CB (corresponding to the cryptographic block) which is datadependent, and a source S (such as clock, or an always on analogue block), which is data-independent. In that case, if a probe placed close to the IC surface but far from CB and close to S collects: only a small fraction of the data-dependent emissions radiated by CB (since the magnetic field amplitude decreases rapidly as the square or the cube of the distance [START_REF] Dhia | Electromagnetic Compatibility of Integrated Circuits: Techniques for low emission and susceptibility[END_REF]) and large portion of the emissions of S. On the contrary, if the probe is placed really close to CB, the probe collects a large fraction of the data-dependent emissions radiated by CB during its operation and a small portion of the emissions of S. From the considerations above, one may conclude that positioning the probe close to the coprocessor CB, results in collecting time domain traces differing significantly one from another, and therefore, that it is easy to localize the cryptographic module. However, these last claims do not hold ! Indeed, if the EM emissions radiated by the source S are significantly greater than the emissions of the CB or, if S and CB sources are really close one from the other, it is still extremely difficult to localize CB, even if the probe is placed just above, since most of the signal collected by the probe is generated by S. As a result, any localization technique based on time domain amplitude analyses may not work in presence of large EM data-independent emanations sources such as clock generators, PLL or I/O interfaces, or even in presence of large environmental electromagnetic noise sources. As an illustration, Fig. 1c shows a map revealing the maximum amplitude of the EM emissions measured at several coordinates of the FPGA package surface during a DES ciphering (this EM emission maps has been obtained with the experimental setup introduced Section 3). Fig. 1b is an X-Ray photography of the package containing the circuit under attack. Finally, Fig. 1a discloses the routing (obtained with Xilinx ISE tool suite) of the considered circuit, running at 50MHz and integrating a DES module, a finite state machine and a RS232 interface for communication purpose. As shown, it appears impossible to correlate Fig. 1a with Fig. 1c even if the die area is roughly known thanks to the X-Ray photography Fig. 1b. It thus appears all the more difficult to identify the DES module on this kind of EM cartography and thus to decide where to position the magnetic sensors above the package to perform a successful CEMA only on the DES module, in order to avoid potential global hardware countermeasures such as [START_REF] Shamir | Protecting smart cards from passive power analysis with detached power supplies[END_REF].

B. About the behavior of EM emissions of IC

DPA exploits by statistical means the data-dependent behavior of the switching current consumed by circuits during a computation of a cryptographic module. This behavior is due to inner properties of the CMOS logic which consumes energy (much more than in the idle state) only to switch from one logical state to another [START_REF] Yeap | Practical Low Power Digital VLSI Design[END_REF]. EM emissions of a circuit are mainly generated by flows of electrical charges through the different metal wires connecting logic gates but also trough wires supplying the circuit [START_REF] Dhia | Electromagnetic Compatibility of Integrated Circuits: Techniques for low emission and susceptibility[END_REF] [START_REF] Ordas | Near-Field Mapping System to Scan in Time Domain the Magnetic Emissions of Integrated Circuits[END_REF]. Since the switching of gates generates a current flow through the circuit interconnect, we may conclude that these switching generate some data-dependent EM emissions at different points in the circuit according to the power distribution network [START_REF] Dhia | Electromagnetic Compatibility of Integrated Circuits: Techniques for low emission and susceptibility[END_REF] [START_REF] Ordas | Near-Field Mapping System to Scan in Time Domain the Magnetic Emissions of Integrated Circuits[END_REF]. These data-dependent behaviors may be exploited by statistical means, using for example DEMA [3][4], to retrieve the secret key. Even if the magnitudes of both power consumption and computation time of logic cells are roughly known, it is extremely difficult to deduce any characteristic about the EM emissions generated by gates due to the complexity of the power distribution grid of actual IC. As a result, the only conclusion we may draw and consider in the remainder of the paper is that gates generates some EM perturbations and more precisely generates some data-dependent harmonics located somewhere in the whole EM emission spectrum. Within this context, the proposed technique allows disclosing the data-dependent behavior of EM emissions in the frequency domain without making any assumption on the EM emission characteristics. It is based on spectral incoherence analysis of two time domain signals as detailed below. The only observation on which is based the method is the following: considering two successive hardware operations, we are sure that some gates switch during one computation and does not switch during the other, while some gates switch during both operations. This leads to the following intuitive conclusion that guides the development of our proposal: between two cryptographic operations some characteristics of the EM emissions remain constant (coherent) from one operation to another, while some character-ristics completely change (are incoherent). Such a behavior is said partially data-dependent in the rest of the paper, and the proposed W GM SI technique [START_REF] Dehbaoui | Enhancing Electromagnetic Attacks using Spectral Coherence based Cartography[END_REF] aims at disclosing circuit areas characterized by this behavior.

C. Coupling W GM SI criterion with EM near field scanning

Coupling W GM SI analysis with Near Field Scanning system to localize the points characterized by partially datadependent EM emissions and thus leaking spots is straightforward. The basic idea is to collect for each (X,Y) coordinates above the integrated circuit at least two different time domain traces of the magnetic field corresponding to two different data processing. Finally, W GM SI values are computed for all (X,Y) positions. This provides W GM SI cartographies revealing positions characterized by partially data-dependent EM emissions. Note however, that computing W GM SI values for more than two data and averaging the results is not theoretically required but may lead to better results in practice.

III. EXPERIMENTAL RESULTS

To validate the effectiveness of the W GM SI analysis in localizing spots with partially data-dependent EM emissions, 2 kinds of validation were performed. The one aimed at correlating the obtained W GM SI cartographies with the design floorplans while the second aimed at demonstrating that spots characterized by the highest W GM SI values are good candidates for CEMA.

A. Analyzed Design

The two aforementioned validation steps were performed considering a design mapped onto a FPGA circuit and more precisely a Spartan3-1000 Digilent board. Note the Spartan die is encapsulated in a cavity-up Ball Grid Array (BGA) package. The mapped design integrates a RS232 block to communicate with the PC, a finite state machine that manages the communications and the behavior of the chip. Three different floorplans of this design were elaborated with ISE tool suite as shown Fig. 4. This was done to definitively validate the efficiency of W GM SI map in disclosing area with partially data-dependent radiations.

B. W GM SI maps vs. design floorplans

W GM SI maps aim at disclosing (X,Y) coordinates at which the EM emissions captured by the probes are partially data-dependent. As a first validation, we therefore scanned the whole package surface during a DES ciphering operations and computed the W GM SI maps for the 3 considered mappings. The expected result was, at least, to localize the DES module and more precisely to define an area (characterized by higher W GM SI values) containing the DES module or part of it since (a) the probe was placed at roughly 500µm from the die and (b) since EM waves are dispersive. Note that the distance from the probe to the die was estimated by removing the package of an identical FPGA and measuring the package thickness. This procedure also allowed measuring the die dimensions (incorporating the IO pads): roughly 8mm by 8mm. The whole package was scanned with a displacement step of 500µm. This resulted in acquiring EM emissions at 1156 coordinates for each mapping. It took 2 hours, for one mapping, to collect the EM emissions corresponding to 5 different encryptions. Note that to increase the Signal to Noise Ratio, each ciphering was performed 20 times and the average computed. Fig. 3 shows the 3 obtained W GM SI maps. The lower left map of Fig. 3 is to be compared to Fig. 1c. This comparison demonstrates that W GM SI map provides more valuable information than a maximum time domain amplitude analysis. Moreover and as expected the comparison of these maps with the corresponding floorplans demonstrates, considering the accuracy of the die size measure (0.5mm), that W GM SI criterion allows disclosing the hardware area corresponding to the DES module. One surprising point is that areas with partially data-dependent EM emissions have different size even if the DES module occupied roughly the same number of slices. This is especially true for W GM SI maps corresponding to mappings 1 and 2 that disclose significantly larger and smaller areas than that effectively occupied by the DES module. This is probably due to the power/ground network specificities effectively involved in the supply of the DES module in each mapping; supply rails being recognized as important sources of EM emanations [6][8].

C. W GM SI and Frequency Bandwidth

As generally recognized CMOS gate dissipates power, this power dissipation is strongly affected by non-zero input signal transition times τ in . Depending on F ast or Slow input signal settling times [START_REF] Maurine | Transition time modeling in deep submicron cmos[END_REF], and according to the Biot-Savart Law, the gate radiation will be spread in the frequency domain, depending on τ in (1),(2) :

|Env B | f ast = µ0
4π .I f ast max .τ in .sinc 2 a (ω. τin 2 ). (1 + 2.cos(ω.τ in )).(sinθ 2sinθ 1 )

(1)

|Env B | slow = µ0 4π .I slow max .τ in .sinc 2 a (ω. τin 2 ). (sinθ 2 -sinθ 1 ) (2)
Where µ 0 , τ in , θ 1 and θ 2 , are respectively the permeability of free space, the ramp duration, and angles involved in Biot-Savart law. As shown in [START_REF] Maurine | Transition time modeling in deep submicron cmos[END_REF], gates obey to two types of switching processes, fast and slow ones. Although both processes may exist within the same circuit, the fast one is much more common within typical circuit [START_REF] Lasbouygues | Logical effort model extension to propagation delay representation[END_REF]. Such switching mechanisms, are characterized at first order by a current consumed taking the shape of a trapezoid for fast processes, and the shape of a triangle for slow ones.

Considering F ast and Slow input signal settling times, Fig. 3 gives the power spectral density of the power radiated by a CMOS inverter, and the measured one using a perfect differentiator (magnetic loop). As observed, for slow switching process, amplitudes measured with a differentiator sensor are more significant at lower frequencies, while, for fast switching process, the highest amplitudes are localized at higher frequencies. 4, where T ck, Slack and I block max are respectively the clock period, the timing margin design and the maximum switching current consumed by the block. The resulting spectral density is given by (3) and represented in Fig. 5.

|Env B | block = µ 0 4π .I block max .τ.sinc 2 a (ω. τ 2 ).(sinθ 2 -sinθ 1 ) (3)
However, the values involved in this expression are very different from those considered for a gate. Indeed, for a combinational block I max and τ are much greater. So to go further in our reverse engineering analysis, we performed static WGMSI maps for different frequency bandwidths. Fig. 6 shows the considered design placement and the corresponding peak to peak map.W GM SI maps for the considered bandwidths ( IV. CONCLUSION

A new technique to localize Electro-Magnetic hot-spots in hardware cryptographic module was presented in this paper. It is based on the assumption according to which: from a data processing to another one, EM emissions radiated by an integrated circuit have some coherent characteristics and some incoherent characteristics. This claimed property, called partially data-dependence of EM emissions, has been first validated experimentally by verifying its correctness with some measured traces. Then, we deduced from this observation, a localization technique. This technique allows (a) localizing cryptographic modules and more precisely leaking points thanks to EM near field mapping and (b) selecting rationally a reduced set of points of interests for electromagnetic analyses. Finally, concrete results have been given on an iterative DES mapped on a FPGA. These results have demonstrated the interest of using incoherence analysis of EM emissions.

Fig. 1 .

 1 Fig. 1. (a) Design floorplan (b) X-Ray (c) Peak to Peak map

Fig. 2 .

 2 Fig. 2. Circuit floorplans and related full package W GM SI maps obtained for the three considered mappings of the design

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. First-Order Power Spectral Densities a) radiated by a CMOS inverter b) collected with a perfect differentiator

  [200MHz -400MHz] and [800MHz -900MHz]) are shown in Fig. 7. As expected, supply rails and pads appear clearly in the frequency bandwidth [200MHz -400MHz], while, logical blocks appear in the frequency bandwidth [800MHz -900MHz].

Fig. 6 .Fig. 7 .

 67 Fig. 6. a) Design Placement b) Peak to Peak map


Nejmeddine Ouhichi
email: nejmeddinouhichi@gmail.com

Radhouane Hamdi

Fethi Lachaal

Hakim Gabtni

Olivier Grünberger

Nejmeddine Ouhich

Olivier Grunberger

Geophysical and hydrogeological investigations of water leakage from a reservoir dam to a coastal aquifer: the Lebna Case Study (Northeastern Tunisia)

Keywords: Recharge, surface water-groundwater interaction, hydrological balance, resistivity, Electrical conductivity, semiarid zone

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In Mediterranean countries, dams play a key role in development. Dams are generally built to provide supplementary water that can satisfy the needs of the population, mainly including urban, water flood control and rural irrigation demands. Dams favor the regulation of the hydrologic regimes of rivers and retain sediments [START_REF] Altinbilek | The role of dams in development[END_REF]. Nevertheless, they may also cause considerable water losses by evaporation, thereby increasing the salinity of the impounded water and lowering the volume of available water needed for groundwater recharge and alternative uses. Simultaneously, such structures, in turn, may promote the recharge of aquifers. According to de Vries and [START_REF] De Vries | Groundwater recharge: An overview of processes and challenges[END_REF], the aquifer recharge is achieved via two key processes: the direct infiltration of precipitation through the soil in outcrop areas and the localized recharge arising from surface water accumulation and sideways percolation (as is the case in this study).

Localized recharge in reservoir lakes depends on various processes and mechanisms [START_REF] Simmers | Groundwater recharge: An overview of estimation 'problems' and recent developments[END_REF], and several authors have considered this process to be an indirect supply process [START_REF] Nazoumou | Impact des barrages sur la recharge des nappes en zone aride: étude par modélisation numérique sur le cas de Kairouan (Tunisie centrale[END_REF]. Thus, recharge can occur in natural lakes or be caused by dams.

Therefore, the recharge process relies on many different phenomena, such as dam releases, infiltration wells, and irrigation returns. Leakage through geological embankments that enclose reservoirs can add to silting and evaporation and may result in considerably high water losses [START_REF] Mesney M | Un grand barrage à haut risque. Saddam dam en Irak[END_REF] and pose severe threats to the stability and security of dam structures [START_REF] Benfetta | Les fuites d ' eau dan s les barrages dans le monde: Quelques exemples algeriens water leaks in dams in the world[END_REF]). The natural recharge by dams is disrupted in some way by climate change which is increasing from year to year and has short and long-term negative effects on surface water and groundwater in the northern part of the African continent. Since among these consequences the irregularity of precipitation, water degradation, stagnation of sludge in dams, dead slice, drought, and increased extraction of groundwater, etc. [START_REF] Hamed | Climate impact on surface and groundwater in North Africa: a global synthesis of findings and recommendations[END_REF] Leaks occur through materials geological boundaries with a degree of porosity, permeability or discontinuity sufficient to allow the occurrence of water movement [START_REF] Kim | Assessment of water seepage through a geologic barrier surrounding a large reservoir using groundwater levels, soil condition, and a numerical model[END_REF]. Leakage is considered a negative effect when it promotes the instability of a dike and/or a significant loss of water; however, leakage can also be considered a positive effect when a nearby aquifer is recharged. The amount of available research on this subject is limited, but examples of reservoir leaks, such as leaks in the such as the Negratín Dam (Spain), Gauting (China), Karoun (Lebanon), and Conqueyrac (France), were cited by [START_REF] Benfetta | Les fuites d ' eau dan s les barrages dans le monde: Quelques exemples algeriens water leaks in dams in the world[END_REF]. Recent works have focused on studying dam leakage processes by quantifying, qualifying or even spatially and temporally delimiting leakage terms (Al-Fares 2011; [START_REF] Berhane | Geological and geophysical investigation of water leakage from two micro-dam reservoirs: Implications for future site selection, northern Ethiopia[END_REF][START_REF] Berhane | Water leakage investigation of micro-dam reservoirs in Mesozoic sedimentary sequences in Northern Ethiopia[END_REF][START_REF] Gutiérrez | Leakage problems in dams built on evaporites. The case of La Loteta Dam (NE Spain), a reservoir in a large karstic depression generated by interstratal salt dissolution[END_REF][START_REF] Kim | Assessment of water seepage through a geologic barrier surrounding a large reservoir using groundwater levels, soil condition, and a numerical model[END_REF][START_REF] Kim | Hydrogeochemical interpretation of water seepage through a geological barrier at a reservoir boundary[END_REF][START_REF] Metwaly | A hydrogeophysical study to estimate water seepage from northwestern Lake Nasser, Egypt[END_REF][START_REF] Sevil | Sinkhole investigation in an urban area by trenching in combination with GPR, ERT and high-precision leveling. Mantled evaporite karst of Zaragoza city, NE Spain[END_REF]. Research on the processes associated with leakage through geologic barriers or natural banks is complex due to limited knowledge of the geometry and properties of subsurface rocks. Moreover, because different causes cause leakage in different dams [START_REF] Benfetta | Les fuites d ' eau dan s les barrages dans le monde: Quelques exemples algeriens water leaks in dams in the world[END_REF], methods must be adapted to the target case study.

Various techniques have been used to study leakage water entering the boundaries of lakes or dams. Some authors have used isotopic tracers [START_REF] Alazard | Etude des relations surface -souterrain du système aquifère d ' El Haouareb ( Tunisie Centrale ) sous contraintes climatiques et anthropiques[END_REF][START_REF] Grünberger | Quantification of water exchange between a hill reservoir and groundwater using hydrological and isotopic modelling (El Gouazine, Tunisia)[END_REF][START_REF] Yi | Detecting the leakage source of a reservoir using isotopes[END_REF] to detect the contribution of lake water to underground water. Leakage estimates have been computed using surface water budgets, Darcy's law or even by correlating reservoir lake levels with underground water piezometric levels [START_REF] Bouteffeha | A water balance approach for quantifying subsurface exchange fluxes and associated errors in hill reservoirs in semiarid regions[END_REF]. Likewise, researchers have used geophysical techniques such as vertical electrical sounding (VES) or electrical tomography [START_REF] Bièvre | Improvement of 2D ERT measurements conducted along a small earth-filled dyke using 3D topographic data and 3D computation of geometric factors[END_REF] to study flow path extents and to delineate leakage areas.

In Tunisia, the government constructed numerous small earthen reservoirs for irrigation, artificial groundwater recharge, and sediment retention. According to the "General Direction of Dams and Large Hydraulic Works" (DGBGTH) the state office in charge, numerous dams in the country are suspected to experience leakage (e.g., Beni M'Tir, Kasseb, Barbra, Bou Heurtma, and El Aroussia). Nevertheless, only a few scientific studies have been performed in this region (e.g., [START_REF] Alazard | Etude des relations surface -souterrain du système aquifère d ' El Haouareb ( Tunisie Centrale ) sous contraintes climatiques et anthropiques[END_REF][START_REF] Dassi | Identifying sources of groundwater recharge in the Merguellil basin (Tunisia) using isotopic methods: Implication of dam reservoir water accounting[END_REF]. The particular case of the Lebna Dam located in Cap Bon (Tunisia) draws attention. Built to deliver surface water to large coastal irrigation regions to diminish the overpumping of the coastal aquifer [START_REF] Gaubi | A GIS-based soil erosion prediction using the Revised Universal Soil Loss Equation (RUSLE) (Lebna watershed, Cap Bon, Tunisia)[END_REF][START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF], the dam quickly showed evidence of noncontrolled recharge occurring through geological material detected through the modification of piezometric levels. Since its establishment, the Lebna reservoir has never dried up, indicating that the leakage does not endanger the surface water supply. Nevertheless, in light of the available evidence representing this situation, the effectiveness of this recharging system controlled by natural spills in front of an underground salinization region experiencing sea water intrusion is uncertain.

The main objectives of this study are to depict the geological barrier formation by employing fieldwork observations and disclosing the hydrological features of the infiltration processes across the dam bank. To determine the infiltration rate and the ratio between the increase in the Lebna reservoir water level and leakage, data collected from the reservoir were used to characterize groundwater flow through field sampling data. Finally, the aquifer geometry was established, and a surface water-groundwater relationship was defined in the study area.

STUDY AREA

The Lebna watershed situated in Tunisia (Fig. 1a) has an area of approximately 210 km 2 and is located in the eastern coastal plain in northeastern Tunisia (Fig. 1b). The watershed is bounded on the northwest by the Abdurrahman Mountains, which reach elevations of 650 m.a.s.l, to the northwest by the town of Menzel Temime, to the southwest by the town of Tafelloune, and to the southeast by a slightly marked relief, occupied mainly by annual crops and orchards at elevations close to sea level (Fig. 1c).

[Insert Figure 1]

In the southern part of the Lebna watershed, a dam was built at the junction of the El Ouedyenne and Bou Dokhane rivers. The Lebna Dam is the largest constructed dam in Cap Bon; it has an area of approximately 2.822 km 2 , and its capacity is almost 30.4 Mm 3 ; its average yearly water supply reaches 20 Mm 3 .

As in many coastal regions, overpumping has caused the water table to be overexploited, thus leading to marine intrusion even before dam construction [START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF]. Several studies have aimed to detect, delineate, or even model this marine intrusion affecting the Korba coastal aquifer. Some of these works include Ziadi rt al. (2019), [START_REF] Ziadi | Use of time-domain electromagnetic (TDEM) method to investigate seawater intrusion in the Lebna coastal aquifer of eastern Cap Bon, Tunisia[END_REF], [START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF], [START_REF] Kouzana | Seawater intrusion and associated processes: Case of the Korba aquifer (Cap-Bon, Tunisia)[END_REF], [START_REF] Kouzana | Geophysical and hydrochemical study of the seawater intrusion in Mediterranean semi arid zones. Case of the Korba coastal aquifer (Cap-Bon, Tunisia)[END_REF], [START_REF] Zghibi | Understanding groundwater chemistry using multivariate statistics techniques to the study of contamination in the Korba unconfined aquifer system of Cap-Bon (North-east of Tunisia)[END_REF]; [START_REF] Zghibi | Assessment of seawater intrusion and nitrate contamination on the groundwater quality in the Korba coastal plain of Cap-Bon (North-east of Tunisia)[END_REF]; [START_REF] Zghibi | Groundwater modelling and marine intrusion in the semi-arid systems (Cap-Bon, Tunisia)[END_REF][START_REF] Paniconi | Modeling and analysis of seawater intrusion in the coastal aquifer of Eastern Cap-Bon, Tunisia[END_REF].

MATERIALS AND METHODS

Hydrological Balance

Surface water balance [START_REF] Thornthwaite | Instructions and tables for computing potential evapotranspiration and the water balance[END_REF] methods are not frequently used because, at the end, the recharge rate is considered to be equal to the difference between the terms characterizing the decrease in the reservoir volume and evaporation effects, and these terms are often uncertain. However, in our case, this method can lead to satisfactory results, as we have a well-calibrated station on the Lebna dam site at which regular measurements are taken.

In this study, the determination of the hydrological balance of the Lebna Dam is based on the daily water volume conservation method [START_REF] Albergel | Les lacs collinaires en Tunisie: enjeux, contraintes et perspectives[END_REF]. This calculation must consider all the water inputs and outputs from the reservoir [START_REF] Bouteffeha | A water balance approach for quantifying subsurface exchange fluxes and associated errors in hill reservoirs in semiarid regions[END_REF][START_REF] Hatch | Spatial and temporal variations in streambed hydraulic conductivity quantified with time-series thermal methods[END_REF][START_REF] Schmadel | Deducing the spatial variability of exchange within a longitudinal channel water balance[END_REF]. Data on precipitation, evaporation, pumping, dike leakage, spilling, draining, and pumping for the National Water Supply and Distribution Company (SONEDE) have been collected from the Lebna Dam and used with water level measurements taken in the reservoir to assess the water balance. We assume that when no rain or runoff occurs (most days in the Mediterranean climate), the lake experiences no inflow, and the decreases are explained only by evaporation, the recharge (or natural leakage) volume through the geological materials, and other withdrawals. Then, an estimate of the daily leakage [m 3 day -1 ] can be calculated as follows:

V leak = ∆V-V ev -[V spill +V Pump +V SONEDE +V gate + V Seepage ] [1]
Where:

∆V is the stock variation in the dam reservoir from day n to day n+1 during a dry period [m 3 day -1 ];

V ev is the volume of water evaporated from the reservoir surface [m 3 day -1 ];

V spill is the volume of water passing over the spillway (when the reservoir is full) [m 3 day -1 ];

V Pump is the volume of water directly pumped for irrigation [m 3 day -1 ];

V SONEDE is the volume of water pumped by the National Water Supply and Distribution

Company [m 3 day -1 ];

V Leak is the daily subsurface exchange volume from the dam reservoir to the aquifer (dam reservoir-groundwater interaction) [m 3 day -1 ];

V gate is the output drainage volume [m 3 day -1 ] (the dam operator prepares for floods by opening the outlet gate to prevent the spillway from working); and

V seepage is the seepage volume [m 3 day -1 ] under the dam's dike. This value is computed as a power law fitted on a reserved flow plume and as a function of the water body change in the reservoir as follows:

V seepage = 86400*(Q day(n+1) +Q day(n))/2 [2] Q = a. Wl b
where a is equal to 0.00196702780619; b is equal to 1.06336220122;

T indicates time;

Q indicates the operating flow rate (m 3 .s -1 ); and

Wl is the water level in cm.

All the variables available for observation were measured during the period from 1990 to 2018 (V, V E , V spill, V Pump , V SONEDE , V gate , and V Seepage ). Rainfall data were used only to determine dry periods. The rainfall records were obtained from a rain gauge located at the dam site.

Two slug tests were performed in the two rivers, El Ouedyenne and Bou Dokhane, at the entrance of the reservoir dam to assess the permeability and the water exchanges between the rivers and groundwater upstream of the Lebna reservoir.

Water balance assessments can be applied on different scales depending on the target resolution. For the Lebna Dam, analysis of the water balance process begins using daily water [Insert Table 1]

Geological Survey

To complement the existing geological maps, a field mapping survey of the Lebna Basin was performed in August 2019 at a scale of 1:10000 to determine the levels permitting natural leakage from the Lebna Dam. The survey was based on previous regional works regarding sedimentary deposits in the Cap Bon region (northeastern Tunisia) and included sedimentological aspects [START_REF] Colleuil | Etude stratigraphique et néotectonique des formations néogènes et quaternaires de la région Nabeul-Hammamet[END_REF][START_REF] Salem | Contribution à la connaissance de la géologie du Cap Bon: stratigraphie, tectonique et sédimentologie[END_REF], biostratigraphic aspects (Hooybergs 1995), tectonic aspects [START_REF] Salem | Contribution à la connaissance de la géologie du Cap Bon: stratigraphie, tectonique et sédimentologie[END_REF]) and stratigraphic and biostratigraphic continuity aspects up to the Hammamet region [START_REF] Temani | Les Ostracodes plio-pléistocènes des séries sédimentaires de la bordure orientale du Cap Bon (coupe de l'Oued Lebna, Tunisie orientale)[END_REF]. Detailed observations were performed on dam bank outcrops and of two series in quarry pit facies (Q1 and Q2) located at approximately 300 and 600 m from the reservoir ridge, respectively, representing useful outcrops for cartography and calibration.

Applied Geophysical Survey

Vertical Electrical Sounding VES

The VES surveys describing vertical variations in the apparent electrical resistivity of geological formations [START_REF] Samouëlian | Electrical resistivity survey in soil science: A review[END_REF] were performed in April 2019 at 68 sites in the study zone (Fig. 1b) using the Schlumberger electrode configuration, which is less time consuming and provides the subsurface geologic layers at a higher resolution compared to the other configurations [START_REF] Loke | A comparison of smooth and blocky inversion methods in 2D electrical imaging surveys[END_REF]. The acquisition was performed using an ABEM TERRAMET~R SAS 4000 along the downstream part of the Lebna reservoir body. The principal objective was to determine the extent of leakage through the geological bank, obtain a global idea of the extension of the recharge zone, and define the geometrical structure of the region downstream of the dam. The AB length extended up to 150 m and reached 300 m in some VES regions during data acquisition, depending on field access.

The determination of layers 'thicknesses and resistivity is manually obtained by relating the VES field curves to the theoretical master curves and auxiliary point charts. The interpretation has been carried out confronting the 1D inversion technique software (IPI2WIN, Russia) and the WINSEV program (Jenny and [START_REF] Borreguerro | WINSEV, programme d'interprétation des sondages électriques verticaux réalisés selon le dispositif Schlumberger[END_REF]). An example of the VES treatment is shown on (ESM_2)

Electrical Resistivity Tomography (ERT)

Geophysical investigations were conducted using ERT profiles at two sites on the downstream edge of the reservoir (Fig. 1b). The locations of these profiles were chosen in the vicinity of the recharge input zone and according to the monitoring wells available for calibration purposes. A grid of four profiles was formed at each of the two sites. A Wenneralpha array with 10-m spacing between electrodes was constructed to collect the resistivity data, reaching an investigation depth of approximately 70 m. The electrical resistivity tomography data are presented in the form of inverted resistivity model including topography obtained using Res2DINV Software. The inversion algorithm is based on the standard smoothness-constrained least-squares inversion algorithm (Gauss-Newton method) (De Groot-Hedlin and [START_REF] Groot-Hedlin | Occam's inversion to generate smooth, two dimensional models from magneto telluric data[END_REF][START_REF] Sasaki | Two new and two resurrected species of the sciaenid genus Johnius (Johnius) from the West Pacific[END_REF]. For all models, the root mean square (RMS) error of the inversion in the resistivity model is less than 5%, which is well acceptable and the models are tolerable (ABEM, 2010)

Piezometry and Water Quality

Knowledge of piezometric data is of sharp interest when assessing groundwater flow directions and identifying aquifer recharge zones [START_REF] Chenini | Piezometric levels as possible indicator of aquifer structure: Analysis of the data from Maknassy basin aquifer system (Central Tunisia)[END_REF].

Four piezometric field campaigns were conducted in the Lebna Basin region in the wet (February) and dry (August) seasons of two successive years, 2019 and 2020. Phreatic groundwater levels were measured in 71 available wells distributed in the study zone (Fig. 1b). The resulting isopiezometric lines were interpolated with the kriging method using Surfer 14 software (Golden Software, Golden, CO, USA). The water table piezometric level represented on maps was measured relative to sea level (meters above sea level: m.a.s.l) based on data obtained in the Shuttle Radar Topography Mission (SRTM), on a differential global positioning system (DGPS) survey conducted at some points and from a 1:25000 topographic map of the region. The relations between groundwater and the Lebna reservoir were considered in the interpolation. The dam water level was considered a piezometric level representing the reservoir shoreline.

For the water quality electrical conductivity measurements, 71 wells were monitored over the same periods. This monitoring campaign also concerned the Lebna Dam water, and the resulting data were used to produce spatialized maps to delimit areas affected by the groundwater recharge induced by the impoundment of the dam.

RESULTS

Lebna Basin Lithological Description

Since 1971, numerous investigations have been conducted in this area, including geological surveys, drilling, geological logging, permeability tests, and geophysical investigations.

Specifically, many authors have conducted detailed reviews of the geology of the Cap Bon and Korba regions [START_REF] Abdeljaouad | L'Eocène continental en Afrique du Nord: Essai de corrélations biostratigraphique et paléoclimatique. Impact régional périméditerranéen[END_REF][START_REF] Ennabli | Etude hydrogéologique des aquifères du Nord-Est de la Tunisie pour une gestion intégrée des ressources en eau[END_REF][START_REF] Hooyberghs | Synthèse sur la stratigraphie de l'Oligocène, Miocène et Pliocène de Tunisie[END_REF][START_REF] Salem | Contribution à la connaissance de la géologie du Cap Bon: stratigraphie, tectonique et sédimentologie[END_REF][START_REF] Hooyberghs | Synthèse sur la stratigraphie de l'Oligocène, Miocène et Pliocène de Tunisie[END_REF]. We focused our research essentially on the Lebna Dam and its downstream section (Fig. 1b). The formations outcropping in the study area include Mio-Plio-Quaternary deposits (Ben Salem, 1998). The lower part of the middle Miocene deposits corresponds to detrital deposits, while the upper part is composed of lenticular sandstones and marls with some lignite levels, referred to as the Saouaf Formation. In the study area, deposits from the Upper Miocene are missing due to the widespread erosion that occurred during the Miocene orogeny [START_REF] Ennabli | Etude hydrogéologique des aquifères du Nord-Est de la Tunisie pour une gestion intégrée des ressources en eau[END_REF]. Transgressional marine Pliocene sediments cover unconforming older folded and eroded formations and constitute the main aquifer corresponding to Astian sandstone sediments, called the Porto Farina Formation. These deposits are mainly composed of sandstone, sand, and marl alternations topped by sandstones and sands [START_REF] Damak-Derbel | Le Pliocène marin du Cap Bon (Tunisie): exemple de dépôts gravitaires[END_REF]. This lithology changes laterally to contain argillaceous sands or consolidated sandstones that are more or less argillaceous. Pliocene formations largely outcrop in the northern reach of the Lebna River but are completely covered by Quaternary deposits in the Tafelloune and Diar El Hajjej areas [START_REF] Ennabli | Etude hydrogéologique des aquifères du Nord-Est de la Tunisie pour une gestion intégrée des ressources en eau[END_REF]. These Quaternary deposits are composed of two units. The lower marine facies unit corresponds to sandy limestones containing mollusks. The upper unit is composed of continental sediments [START_REF] Ozer | Essai de corrélation du Pléistocène supérieur de la Sardaigne et de la Tunisie[END_REF]) containing oolitic limestones and coprolites or pelloïdes. These deposits currently form a coastal consolidated wind dune built following marine regression [START_REF] Chakroun | Découverte d'une faune de mammifères du Pléistocène supérieur dans la grotte d'El Geffel (cap Bon, Tunisie)[END_REF]. Old, consolidated dunes cover the Tyrrhenian deposits [START_REF] Ozer | Essai de corrélation du Pléistocène supérieur de la Sardaigne et de la Tunisie[END_REF], forming an essential feature of the littoral topography of the Cap Bon peninsula. This geological investigation of the study area allowed us to establish a Lebna Basin map series (Fig. 2), containing seven units, as follows, from the base to the top (1 to 7):

1. A gypsum gray marl layer;

2. alternating sandstone layers corresponding to a continuity of the Souaf Formation from the upper Miocene, which is thickest in this zone and forms several depositional sequences;

3. a series of 5-to 6-m-thick plastic gray clays, known in the published literature as the Sidi Barka Formation [START_REF] Colleuil | Etude stratigraphique et néotectonique des formations néogènes et quaternaires de la région Nabeul-Hammamet[END_REF]) of the Lower Pliocene; these clays can be considered the floor of the aquifer;

4. upon these clays lies a sedimentary unit with a total thickness of 8 to 10 m, composed of bioclastic sandy levels alternating with lumachelic sandstones and organized in several layers. The base of each layer shows bioclastic sands with fine and homogeneous granulometry and an ochre-yellow color, some of which contain sandy lenses. The fauna appears to be reduced, showing fragments of lamellibranch valves.

This sedimentary unit ends with a very fossiliferous sandstone term containing large Pectinidae and Ostreidae fossils. This stratigraphic unit is known in the literature as the Hammamet Formation consisting of sand and sandstone [START_REF] Colleuil | Etude stratigraphique et néotectonique des formations néogènes et quaternaires de la région Nabeul-Hammamet[END_REF]);

5. the section ends with coarse reddish sand that becomes yellowish following breakage and is fine-grained and beige to pink;

6. the latter end contains magnesian limestone crusts that outcrop extensively in the area.

These crusts likely correspond to the continental Villafranchian described in the literature; and 7. current alluvium sourced from upstream watersheds during flooding events constitutes the reservoir lake floor.

[Insert Figure 2]

Water Balance Analysis

A statistical analysis performed at different lag-time intervals allowed the selection of an optimal interval time of 10 days; this lag time showed the lowest dispersion in the leakage

estimates. An analysis of the curve plotted after calculations showing the average infiltration versus the dam water level (Fig. 3) revealed the following information: the infiltration value rose as the water level increases, and two different slopes exist. One slope was found between the 12.5 and 16-m levels, where infiltration rises very slowly from 0.015 mm.day -1 (3.6×10 - 5 m 3 .s -1 ) to 0.025 mm.day -1 (6.3×10 -5 m 3 .s -1 ); under these conditions, the contribution of subsurface flow is restricted. A second slope was found in which concentrated recharge increases at 16 m and reaches a maximum of 0.09 mm.day -1 (2.1 m 3 . s -1 ) at a water level of 18.04 m.

[Insert Figure 3]

The average effect can be summarized with a curve expression that links the groundwater recharge to the water level in the lake. This expression, applied to the time series of the dam Notably, low annual recharge corresponds to dry years and high annual recharge corresponds to rainy years (Fig. 4). Most likely, a linear correlation exists between the total annual rainfall and the annual recharge. As more rainfall occurs, the reservoir is filled, and more recharge is thus generated. However, the infiltration tendency is in accordance with the reservoir filling behavior. Indeed, the maximum infiltration value (5.65 Mm 3 .y -1 ) was recorded during the 1999-2000 year, when maximum rainfall occurred, and the minimum infiltration value (1.6

Mm 3 .y -1 ) corresponded to 2010-2011, when minimum rainfall occurred. In addition, March, April and May presented the highest dam infiltration values. Based on the infiltration law and by interpolating these values, the infiltrated volume was estimated to be 17.9407 mm/year over the entire watershed. The global computed leakage rate represented approximately 16% of the inflow entering the reservoir. This amount was then compared with the evaporation value of 29.54%. A schematic figure was drawn to visualize each term separately using annual averages (Fig. 5).

[Insert Figure 4]

[Insert Figure 5]

Geophysical Results

Electrical Resistivity Tomography (ERT)

The obtained ERT profiles were inverted to electric tomography resistivity models and highlighted three distinguished zones (Fig. 6 and7). The main objective was to detect the preferential infiltration zones, and the edges of the dam were found to allow more infiltration than the other zones.

[Insert Figure 6]

[Insert Figure 7]
For the first layer, 2 to 3 meters from the surface, resistivity values between approximately 20 and 60 Ω.m were detected, with a clear change from one profile to another due to the topography effect. Some heterogeneity existed due to the presence of the clays or marl beds that had already been mapped around the dam bank. This level was not completely continuous; thus, it was absent in some zones and formed small, convex shapes on some profiles. The second layer showed high resistivity values between 140 and 700 Ω.m, contrasting with the deepest part of the ERT profile records corresponding to the clay bedrock (with resistivity values lower than 20 Ω.m). The described layers had the same order of resistivity as the highest topographic part of the profiles, although their depths varied. This layer also showed a highly variable thickness.

Lateral changes from high to low resistivity, especially near the surface, may have been caused by discontinuities and/or abrupt facies changes. The integrated interpretation of the data aimed to distinguish among the infiltration zones at the two different investigated sites to identify the zone that was the most susceptible to leakage. The water levels on the day's measurements were taken were well marked on the resistivity profiles. The area above the water level traced on the profiles showed very high resistivity values, belonging to the damaquifer contact zone with the sandy sandstone formation of the upper Pliocene. The contact areas permitting aquifer recharge at the first site (site 1) were more frequent than those at the second site (site 2).

Vertical Electrical Soundings (VES) Interpretation

The determination of layers 'thicknesses and resistivity is manually obtained by relating the VES field curves to the theoretical master curves and auxiliary point charts. The interpretation has been carried out confronting the 1D inversion technique software (IPI2WIN, Russia) and the WINSEV program (Jenny and [START_REF] Borreguerro | WINSEV, programme d'interprétation des sondages électriques verticaux réalisés selon le dispositif Schlumberger[END_REF]. The VES inversion provided an electrical image of the underground structure from the surface up to a depth of 100 m. six maps were selected among 14 maps produced (ESM_2). The maps represented the apparent resistivity variation over the whole zone covered by the Lebna Basin along a section of subsoil varying in thickness between 3 and 75 m. The resistivity range denoted on the maps varied between 0 and 700 Ω.m, and high resistivity values [> 90 Ω.m] were concentrated downstream of the reservoir lake. The apparent resistivity showed circular to subcircular shapes ranging from green to red on the scale and indicated either the existence of permeable layers allowing recharge near the dam or even zones farther away that were convenient for the implementation of water drilling.

In the area near the dam, on the right bank side, the maps show high-resistivity zones (>70 Ω.m) at different AB/2 locations in contact with the reservoir. Consequently, the dam bank is surrounded by high-value contours, allowing water to flow and partly explaining the water exchange with the Korba aquifer.

The existing borehole BH1 allowed the calibration of the proximal VES, S2-5 (ESM_3).

However, the available lithological columns are reconstructed from field observations and calibration based also on two quarries pit facies Q1 and Q2. Interpretation of the rest of the VES will be based on the relationship between resistivity and facies obtained from the BH1 calibration providing the necessary data to develop isobaths and isopach maps of the aquifer delivering information on aquifer geometry.

As result, the top isobath map of the shallow aquifer (Fig. 8.a) represents deformations that affect the Plio-Quaternary deposits in the study site. This map shows a NE-SW structure located approximately 500 m south of the lake bank, corresponding to a high structure bordered by a lower area (''H2'' on the diagram). Faults or monoclines related to the late Quaternary tectonic extension could have generated this structure; these structures could be studied by seismic refraction, as was in the study by [START_REF] Burollet | Structures and tectonics of Tunisia[END_REF]. The resulting isobath map revealed an important element for consideration: the high structure existing in the subsurface on the southern side of the dam. This structure must also be considered to characterize the hydrogeological system. ''H1'' and ''H2'' on the map also denote high structures, but their elevations result from topography effects, as they are located farther upstream with higher inclinations. The isopach map of the shallow aquifer (Fig. 8.b) shows that the thickness of the sandy aquifer ranged between 5 and 10 m in large areas but was thicker in certain places. A thick zone was located directly south of the dam, enabling a preferential flow path toward the SE direction in accordance with the topographic slope toward the sea.

[Insert Figure 8]

Piezometric and Electrical Conductivity Results

The general behavior of the piezometric levels of free aquifers is mostly related to topographic variations [START_REF] Chenini | Piezometric levels as possible indicator of aquifer structure: Analysis of the data from Maknassy basin aquifer system (Central Tunisia)[END_REF]. This is the case for the piezometric levels of the sandy aquifer studied herein, the piezometric levels of which are mainly oriented parallel to the seashore. By reviewing piezometric maps of 2020 (Fig. 9), we detected a southeast general flow direction exiting the southeast bank of the reservoir. The flow path directions show discrepancies in their general directions that are in accordance with the aquifer geometry heterogeneity determined in previous research. Indeed, the average piezometric gradients range from approximately 3‰ and can drop to inferior levels of 5‰ in areas where the aquifer is particularly thick (>20 m); this trend was observed in all 4 seasons, even in 2019, (ESM_4). Near the sea shore (<2 km distance), piezometric levels below sea level illustrate the influence of overpumping above sea water intrusion, as previously studied by researchers [START_REF] Kingumbi | Évaluation des transferts entre barrage et aquifères par la méthode de bilan d'une retenue en zone semi-aride[END_REF][START_REF] Kouzana | Seawater intrusion and associated processes: Case of the Korba aquifer (Cap-Bon, Tunisia)[END_REF][START_REF] Zghibi | Assessment of seawater intrusion and nitrate contamination on the groundwater quality in the Korba coastal plain of Cap-Bon (North-east of Tunisia)[END_REF][START_REF] Ziadi | Mineralization and pollution sources in the coastal aquifer of Lebna, Cap Bon, Tunisia[END_REF]).

[Insert Figure 9]

The underground mapping of electrical conductivity (EC) (Figure 10,ESM_5) illustrated EC values ranging from 0.4 to 7.5 mS.cm -1 . The lowest EC values identified in the groundwater were close to those observed in the dam water: 0.8 mS.cm -1 in winter and 1.3 mS.cm -1 in the dry season. Around the reservoir banks, the ECs were usually within these values up to a distance of 1 km during winter and a distance of 500 m during the dry season. Excluding the August 2020 situation, the EC values identified at some distances from the reservoir seemed to stabilize between 1.8 and 2.6 mS.cm -1 (corresponding to total dissolved solids (TDS) values of ~0.9 to ~1.3 g/L) in large areas. These values are characteristic of mixing between the reservoir water and underground water, corresponding to a TDS of 2 g/L (EC ~4 mS.cm -1 ) [START_REF] Hamouda | Origine de la minéralisation dans l'aquifère plio-quaternaire de la côte orientale du cap Bon (Tunisie)[END_REF]). In the southeastern part of the study area, the EC values were generally high (> 3.5 mS.cm -1 ) due to saline intrusion. Similarly, 2 wells in the northwestern part of the study area showed EC measurements greater than 3.5 mS.cm -1 during the 4 seasons sampled (wells W60 & W61). In August 2019, a clear and continuous dilution plume was observed coming from the dam side, while in August 2020, the plume seemed to have been displaced toward the south, as pronounced increases were identified in EC in the wells located at distances between 1 and 2 km from the Lebna Dam.

[Insert Figure 10]

DISCUSSIONS

A geological survey was performed to identify and map a sandy formation outcrop known as the Hammamet Formation of upper Pliocene age; this material is present on the southern ridges of the studied reservoir. The position of the formation lies on a layer of Sidi Barka clay (Lower Pliocene), favoring infiltration when water levels are higher than 15 m.a.s.l. The rise of the water level to the dam-dumping height of 18 m.a.s.l should increase at the same time the infiltration surface area and the hydraulic gradient increase. Dam balance calculations indicated that an average water volume of 3.76 Mm 3 .y -1 +/-0.02 infiltrated into the water table downstream of the dam throughout the 1990-2017 period, with high variations occurring between wet and dry years. As mentioned by [START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF] to be coherent with his modeling of the Korba aquifer, the estimated volume was approximately 3.73 Mm 3 .y -1 in 2004. [START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF] used a unique value as a permanent input to represent the coastal aquifer, whereas we were able to compute varying recharge estimates that could be entered to represent scenarios dedicated to forecasting the evolution of saline intrusion under climate change conditions. The extension and limits of the aquifer downstream of the reservoir were precisely determined by geophysical means combining VES and ETR methods. Tectonic feature and facies changes influenced the aquifer layer flow paths from the reservoir to the sea; for example, the Hamamet Formation presented large thickness variations with clear piezometric patterns. Areas with large aquifer thicknesses implied low piezometric gradients.

This description of the aquifer limits constitutes added value that creates opportunities for the simulation and understanding of the responses of the aquifer to recharge and pumping. From the EC and piezometric evolution data, 4 different areas of interest were discriminated:

1-Near the reservoir (0-500 m from the south banks), low EC (<1.8 mS.cm -1 ) and piezometric levels were connected to the corresponding lake levels with low gradients, indicating the direct influence of the recharge zone in the Hammamet Formation.

2-Downstream (0.5-1 km), tectonic structures narrow the flow paths, and a higher-gradient zone is found with higher EC values [1.8 to 2.6 mS.cm -1 ]. The salinity in this region is typical of mixing between underground water and reservoir water 3-Downstream (1-2 km), a thicker aquifer zone with low piezometric gradients and similar EC values can be found.

4-

At approximately the sea shoreline (distance<3 km), the piezometric gradients are disturbed, and the EC values are elevated (>2.6 mS.cm -1 ), indicating overpumping and salt intrusion.

A combination of VES and ERT techniques was used in this study, but as Kayode, Odukoya, Adagunodo, and Adeniji (2018) they found that a combination of seismic refraction and resistivity methods is adequate for obtaining robust information regarding seepage pathways in and around dams. Similarly, the self-potential (SP) and induced polarization (IP) methods are known to be sound for monitoring subsurface seepage [START_REF] Panthulu | Detection of seepage paths in earth dams using self-potential and electrical resistivity methods[END_REF]; thus, we considered applying these methods in the future near the lake zone, but these methods may be inadequate for sites with piezometric levels deeper than 4 to 8 meters. In this study, the geophysical campaign was carried out during the dry season. It is always possible to perform regular follow-up measurements at different times to detect the resistivity variations on the banks of the Lebna reservoir.

The uncertainties in our piezometric maps are difficult to reduce; few points were measured by DGPS methods (9 out of 71), and for the rest of the analyzed points, SRTM 90 and topographic map data were used, contributing an accuracy of 2 meters. In our experience incorporating DGPS values with topographic map data, the errors were under 50 cm. Nearly all piezometric data are collected in traditional wells that are very often affected by pumping and/or stand as abandoned wells with suspicious water renewal processes. Pumping for irrigation and water returns were not evaluated in the study. [START_REF] Kerrou | Status of the Korba groundwater resources (Tunisia): Observations and three-dimensional modelling of seawater intrusion[END_REF] estimated an annual pumping value of 100 mm for irrigated zones. We assumed that this value was insufficient for inducing a large water table perturbation. Nevertheless, we recommend complementary temporal monitoring with a pressure probe to document these effects. In most previous works that explored the dam leakage phenomenon [START_REF] Dreybrodt | Karstification below dam sites: A model of increasing leakage from reservoirs[END_REF][START_REF] Zhang | Impacts of the sanmenxia dam on the interaction between surface water and ground water in the Lower Weihe River of Yellow River Watershed[END_REF][START_REF] Zhao | Interactions between dam-regulated river flow and riparian groundwater: A case study from the Yellow River, China[END_REF][START_REF] Zhou | Surface water-groundwater interactions of Xiluodu Reservoir based on the dynamic evolution of seepage, temperature, and hydrochemistry due to impoundment[END_REF], piezometric data were used as the main evidence. In several cases, the data interpretations were based on punctual measurements recorded downstream of the studied dam and correlated to the water level in the reservoir [START_REF] Zhou | Surface water-groundwater interactions of Xiluodu Reservoir based on the dynamic evolution of seepage, temperature, and hydrochemistry due to impoundment[END_REF]. Whenever piezometric maps were presented in the works mentioned above, they were used to describe the flow in general and not the flow perturbations resulting from leakage flow contributions. This consideration resulted in our study being more focused on the detection of zones under the direct leakage reservoir dam effect.

In this study, a relation between the water level and infiltration was established, and the resulting curve included slope changes that corresponded to the limits of sandy layers located in the reservoir bank, as described by geology field surveys. This fact illustrated that an increase in recharge occurs when the lake water level reaches the level of the sandy formations on the southeastern reservoir ridges. The uncertainty in the estimated infiltration volume, obtained through the calculation of some terms in the balance equation, resulted from the manual centimetric measurement accuracy and the simplified working hypothesis regarding the water budget. In contrast, a leakage-water level relationship was established by [START_REF] García-López | Conjunctive use of water resources as an alternative to a leaky reservoir in a mountainous, semiarid area (Adra River basin, SE Spain)[END_REF] in a karstic system. They evidenced a progressive decrease in leakage rates over a certain threshold level and explained this phenomenon by the progressive "saturation of the rock" (the filling of a hydrogeological compartment via a karstic sinkhole).

In this case, the resulting equilibrium between the piezometric level and the water level of the reservoir would explain the decrease in leakage flows. In our case, the reservoir level can be considered a piezometrically imposed limit that conditions recharge through a combination of the nearly impermeable alluvial sediments and highly permeable aquifer levels that are reached when the reservoir water level is higher than 16 m.a.s.l. [START_REF] Boulton | The functional significance of the hyporheic zone in streamsand rivers[END_REF] proposed three spatial scales to represent these interactions: the basin/watershed scale, contact zone scale, and sediment scale; these scales were also considered in this paper. Despite previous studies being conducted at the basin scale [START_REF] Kouzana | Geophysical and hydrochemical study of the seawater intrusion in Mediterranean semi arid zones. Case of the Korba coastal aquifer (Cap-Bon, Tunisia)[END_REF][START_REF] Ziadi | Use of time-domain electromagnetic (TDEM) method to investigate seawater intrusion in the Lebna coastal aquifer of eastern Cap Bon, Tunisia[END_REF], few studies analyzed the recharge evolution over different time periods for dam sites located in basins similar to the Lebna Basin. In addition, most previous research has analyzed surface water-groundwater interactions using single methods rather than combinations of multiple field methods; the use of a combination of methods thus represents the added value of the present work.

CONCLUSION

Water leakage should be considered when constructing of hydraulic structures such as dams, reservoirs or artificial lakes. This phenomenon can contribute to the study region either negatively (through water losses or dam breakage) or positively (by contributing to the recharge of the groundwater), as in our case study. Investigating the geological nature of the bank layers prior to construction is a major step that can be taken to determine whether leakage will occur and, if so, the order of magnitude. In cases such as our study, the leakage volume mentioned (17% of the dam inflow) was not as critical as that reported in [START_REF] Somers | Quantifying groundwater-surface water interactions in a proglacial valley, Cordillera Blanca, Peru[END_REF], who suggested a 49% exchange with the subsurface. Nevertheless, this phenomenon was welcomed, as it represented enhanced contribution to the recharge of an overexploited aquifer affected by marine intrusions.

Recharge through the percolation of dam water into natural beds represents an artificial input source to groundwater in arid areas. The choice of method for estimating recharge must be guided by the research objectives, the data availability and the possibility of acquiring additional data and the scale of the study area. In this paper, we estimated the groundwater recharge volume over a 27-year period in a Mediterranean catchment area using a water balance approach. Electrical resistivity measurements were applied downstream of the Lebna Dam reservoir to precisely detect the locations where leakage occurred. Based on vertical electrical surveys and tomography, resistivity methods allowed the detection of highresistivity regions corresponding to leakage areas.

The current findings are very useful for stakeholders, especially, to enhance the importance of the Lebna dam as a source of aquifer recharge for the upstream region. Characterization of the Lebna dam / aquifer relationship, improved by a non-controlled leakage of 3.7 Mm 3 . y -1 from Lebna dam to upstream aquifer and leakage mapping are calculated and presented for the first time.

The characterizations of water leakage and aquifer (groundwater balance, water level, geometry) are useful for: i) local and national water management authorities (DGRE and CRDA Nabeul), to improve sustainable groundwater management strategy, and ii) groundwater uses specially the farms in upstream dam, to identify the depth of drilling and wells to reach the aquifer and to know the available groundwater resource.

ACKNOWLEDGMENT

This work has been supported by the International Laboratory NAILA (http://lmi-naila.com/)

and ALTOS-PRIMA Project (https://www.altos-project.org). Authors are grateful for the Lebna dam authority (DGBGTH) and staff that provided site access and shared their operational data set.

  level measurements (ESM_1.a). According to the objective of estimating leakage from Lebna reservoir, only dry-season periods were selected, and the daily negative balance values were maintained (ESM_2.b). Due to the difference in scale between the water level measurements (measured on a centimeter scale) and the estimated infiltration corresponding to a few millimeters, it was necessary to conduct several aggregations at 5-, 10-, 15-and 20-day intervals. The balance terms were calculated in the same way. However, the value assigned to the date represented a sum of the 5-, 10-, 15-or 20-day increments. From the water level graph derived as a function of time at the chosen interval (ESM_1.c), we classified years, according to the water level intervals, into eight groups (Tab.1). Each group was assigned an instantaneous average water level value attributed to an average infiltration rate. May, June, July and August represent the recession period. When estimating the infiltrated water volume, this period is preferred. The infiltration volume was estimated according to the surface area of the Lebna watershed.

  water level, allowed us to estimate a volume of 3.76 Mm 3 .y -1 over a period of 27 years. The examination of the reservoir levels from 1990 to 2017 (ESM_1a) showed many periods of water storage and destocking. The water volume stored in the reservoir reflected rapid and frequent increases inFebruary of 1996February of , 1999February of , and 2015 and in November of 1999; these increases were related to high rainfall during the corresponding period. Less-frequent flood events occurred inMarch 2007 and April 2009.

Figure 1 :

 1 Figure 1: Location and geologic map of the study area. (1) Tunisian map; (b) geographical situation and digital terrain model of the Lebna watershed relative to the Cap Bon peninsula and to the country of Tunisia; (c) geological setting based on maps (Sheet N° 22: Realized by M. Arnould (1951), actualized by H. Ben salem (1997) and Sheet N° 15: by Bensalem 1997), with well-sampled, borehole, quarry and VES locations; and (c) ERT profile locations.

Figure 2 :

 2 Figure 2: Geological map of the Lebna Dam basin, illustrated by pictures captured in 2019 depicting the sedimentology and nature of outcropping units in the study area and the AB cross section realized in the downstream region of the Lebna Basin.

Figure 3 :

 3 Figure 3: Infiltration from the Lebna reservoir as a function of the water level.

Figure 5 :

 5 Figure 5: Schematic illustration of the water balance terms of the Lebna Dam.

Figure 4 :

 4 Figure 4: Average monthly infiltration rate from the Lebna Dam reservoir over 27 years (1990-2017).

Figure 6 :

 6 Figure 6: Resistivity models of site 1 including the topography derived by the ERT method and the Wenner array of site 1.

Figure 7 :Figure 8 :Figure 10 :

 7810 Figure 7: Resistivity models including topography obtained using the Wenner ERT method of site 2.

Table 1 .

 1 Annual classification according to the water level interval (measured during the dry

	75
	76	season)
		Group:	1	2	3	4	5	6	7	8
	Interval in m
		[from-to]:

DATA AVAILABILITY

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

CONFLICT OF INTEREST

We are not aware of any conflicts of interest related to this publication.