Industry-Driven versus Natural Groundwater Flow Regime at the Dead Sea Coastal Aquifer

: The coexistence of nature and anthropogenic development requires continuous monitoring and research to address and respond to unforeseen threatening processes that occur with time. This is particularly relevant to the groundwater ﬂow regime in the coastal aquifer adjacent to the Dead Sea, the level of which is dropping, and the industrial evaporation ponds, whose levels are rising. The increasing hydraulic gradient between the two water bodies has produced severe leakage through the pond embankments. To prevent this leakage, a vertical deep sealing wall was built along the embankment. In this study, the overall leakage is calculated by mass balance, and the subsurface leakage component is numerically simulated, based on the mass balance and hydrological observations. Some of the leakage discharges into surface canals and some at the Dead Sea. The leakage volume increased from 20 mcm/year in the 1980s to 100 mcm/year before the sealing wall was built (in 2012), and from 60 mcm/year once the wall was established to 80 mcm/year today. Using the calibrated model, the leakage volume is predicted to increase in the next few decades, mainly through the Ye’elim alluvial fan. Further research effort is needed to come up with new preventive measures.


Introduction
Worldwide exploitation of water resources increased rapidly during the 20th century; today, many water resources are being fully utilized, and even overexploited [1,2]. The world is, therefore, facing a very severe water crisis, affecting developed and developing countries alike. Overexploitation of surface-water and groundwater sources, and pollution of the remaining water sources, has led to ecological disasters in many areas worldwide [3,4]. Control over these processes has only been achieved in a few places, allowing man and nature to coexist [5][6][7]. The Dead Sea basin of Israel and Jordan is a classic example of coexistence between settlement, agriculture, industry, tourism and natural reservations, despite the fact that its water resources are fully utilized; there, for the moment, the conflicts between anthropogenic development and nature is under control.
Utilization of water in the Dead Sea drainage basin expanded dramatically in the latter half of the 20th century and, so far, the Dead Sea level has dropped by approximately 40 m [8,9]. The fresh water that flows to the Dead Sea is entirely used throughout its tributaries by about 25 million people, living in Israel, Jordan, Syria and Lebanon. The water is used for urban consumption, agriculture, industry and tourism. The ultimate solution to the lake's desiccation is to build a sea canal that will stabilize the lake at a constant level [10]. Several alternatives have been proposed for this sea canal ("Med-Dead" from the Mediterranean Sea, and "Red-Dead" from the Red Sea), and one of them is likely to be implemented sooner or later to allow the continued harmonious existence of people and nature in the Dead Sea watershed [11].
Concurrently, on the micro-scale along the Dead Sea shore, several collisions between anthropogenic development and nature nearly spun out of control, but scientific research regained control by implementing appropriate measures toward environmental preservation. For example, stream incision and sinkhole development around the Dead Sea began to damage roads, agriculture, tourism and other civilian infrastructures [12,13]. However, understanding the geomorphological mechanisms [14] and mapping the threatened zones [15] enabled the continuation of human existence by shifting the infrastructure to non-threatened zones.
In the present study, we focus on the conflict between industry and nature in the southern part of the Dead Sea coastal aquifer. The study addresses the groundwater flow regime and demonstrates how the flow field has changed several times due to human intervention in the environment. Two groups of chemical plants were established in this area, one in Israel and the other in the Kingdom of Jordan, for the production of potash and other minerals from the lake brine ( Figure 1). These factories supply agricultural fertilizers to many countries, allowing them to increase food production. In addition, they provide a respectable and stable source of livelihood for the residents of the area. Furthermore, the evaporation ponds of the factories form the infrastructure for many hotels on their shores. Nevertheless, these plants have induced a dramatic change in the lake-aquifer interaction, diverting the groundwater flow field, and at certain points, the evaporation ponds have been in danger of collapse. However, preventive actions were taken at the right time, based on scientific understanding. Similarly, today, changes in the groundwater flow field are once again endangering the industry and tourism and must be properly confronted.

The Dead Sea Region
The Dead Sea basin is a rhomb-shaped pull-apart basin, and the deepest structure along the Dead Sea rift valley [16]. Its marginal faults create a steep cliff above the Dead Sea shore. The Dead Sea (Figure 1) is the lowest land location on Earth; its level was 434 This study demonstrates how hydrological-environmental scientific research and collaboration between all stakeholders will enable maintaining coexistence in the region. In this way, coordination is likely to continue in the future as well, on both the small scale of the chemical plants and the large scale of the entire basin.

The Dead Sea Region
The Dead Sea basin is a rhomb-shaped pull-apart basin, and the deepest structure along the Dead Sea rift valley [16]. Its marginal faults create a steep cliff above the Dead Sea shore. The Dead Sea (Figure 1) is the lowest land location on Earth; its level was 434 m below sea level (mbsl) in 2020. It is a hypersaline and deep terminal lake with a maximum depth of about 300 m and a salinity of~342 g/l TDS (total dissolved solids). In the last few decades, the Dead Sea level has been dropping at an average rate of~1.1 m per year (Figure 2), mainly due to damming of its feeding rivers and exploitation of upstream aquifers [13,17]. This drop is predicted to continue in the next few decades until a new limnological equilibrium is achieved, probably at an elevation of 500 to 550 mbsl [8,18,19]. The southern part of the Dead Sea is much shallower (~10 m deep when the lake level is 400 mbsl), flooding only when the lake level is above 402 mbsl [18]. During most of the 20th century, the Dead Sea level was higher than 402 mbsl and the lake extended southward ( Figure 1B). Since the 1980s, however, the Dead Sea level has dropped below this threshold and the southern basin has naturally dried [14,19]. Instead, artificial saline water reservoirs were established which serve as the evaporation ponds for the Dead Sea Works ( Figure 1C). The water to these ponds is supplied via pumping from the northern basin of the Dead Sea ( Figure 1D). The water level in these ponds is actually rising with time, at a rate of 0.2 m per year (Figure 2), due to halite precipitation on the pond floor and raising of the pond embankment to maintain a water depth of about 2 m for efficient evaporation (The Dead Sea Preservation Government Company data). After~40 years of anthropogenic activity, there are now two saline water bodies,~10 km apart, with a 40-m difference in elevation-the Dead Sea to the north and the evaporation ponds to the south.  Two main aquifers are found west of the Dead Sea: the regional Upper Cretaceous Judea Group, and the local alluvial Dead Sea Group sediments along the Dead Sea coast [20][21][22]. The Dead Sea Group is composed of rift-filling sediments that are relatively young, namely Pliocene and Holocene. The marginal faults of the rift place the Judea Group carbonates (Albian-Turonian) in lateral contact with the Dead Sea Group sediments exposed along the Dead Sea shore. These sediments are roughly divided into two hydrogeological facies [23]: permeable clastic sediments deposited in fan deltas (conglomerate, gravel and sand), and impermeable lacustrine sediments (clay, marl, aragonite, gypsum and salt) deposited between the deltas. Two main aquifers are found west of the Dead Sea: the regional Upper Cretaceous Judea Group, and the local alluvial Dead Sea Group sediments along the Dead Sea coast [20][21][22]. The Dead Sea Group is composed of rift-filling sediments that are relatively young, namely Pliocene and Holocene. The marginal faults of the rift place the Judea Group carbonates (Albian-Turonian) in lateral contact with the Dead Sea Group sediments exposed along the Dead Sea shore. These sediments are roughly divided into two hydrogeological facies [23]: permeable clastic sediments deposited in fan deltas (conglomerate, gravel and sand), and impermeable lacustrine sediments (clay, marl, aragonite, gypsum and salt) deposited between the deltas.

The Evaporation Ponds of the Dead Sea Works
The water that is exploited from the Dead Sea to the evaporation ponds is pumped at station P88, at the southwest corner of the Dead Sea ( Figure 1D). From there, the water flows gravitationally southward in the open "feeding canal" ( Figure 1D) to station P5 at the northwest corner of the evaporation ponds. There, the water is raised to the northern evaporation ponds. The evaporation ponds are roughly divided into two types: ponds for halite precipitation, which is, in fact, a waste material, and ponds for carnallite precipitation, which is the raw material for potash production. The northern ponds are mostly halite ponds ( Figure 1D). In the northern ponds, the water flows slowly by gravitation from pumping station P5 to stations P11 and P33, from which the water is pumped to the southernmost ponds. There, the water again flows gravitationally northward to the end brine discharge point. At this point, the end brine flows out from the ponds through the Arava stream back to the Dead Sea ( Figure 1D,E). As the brine flows along this route, the solutes become more concentrated, the density increases, and the Na/Cl ratio decreases as well.

Leakage of Brine from the Ponds
The Arava stream divides the Israeli evaporation ponds on the west and the Jordanian evaporation ponds on the east ( Figure 1D). Due to the dropping Dead Sea level in the last decades, the Arava stream incision has accelerated (Dente et al. 2015), threatening the stability of the pond embankments. Indeed, significant leaks were detected along the foot of the pond embankments due to the increasing hydraulic gradient between the evaporation ponds and the lake. These threats were, in fact, the direct result of the changing hydrological system, induced by anthropogenic development.
The Dead Sea Works prevented embankment collapse in 2012 using several measures. They enforced the embankments by dynamic compaction and constructed an underground vertical sealing wall-0.5 m thick and 20 m deep-at the center of the embankments. These measures were successful and significantly reduced the leaks. In other words, the natural groundwater flow regime, which was controlled by the increasing hydraulic gradient, was changed again by the industry through the underground seal. In preventing the collapse of the evaporation pond embankments, the coexistence between man and nature was retained.
Nevertheless, in the northwest corner of the ponds, no seal was installed because of the alluvial fan of the Ye'elim wadi ( Figure 1E). This alluvial fan is composed of gravel and conglomerate, contains large boulders to a depth of more than 60 m [24], and does not allow construction of any deep cutoff. This fan is the Achilles' heel in terms of leakage. Indeed, a few hundred meters north of the pond embankment, within this alluvial fan, leaking brine discharges into the marine canal and flows through the Arava stream to the Dead Sea ( Figure 1D,E). The chemical composition of the brine in the marine canal was found to be similar to the chemical composition of the northern evaporation ponds (Supplementary  Materials Table S1), and this leak is increasing with time, again endangering the industry.
Additional leakage from the evaporation ponds occurs through the eastern embankments toward the Arava stream. Certainly, the constructed cutoff stopped most of the leakage, but not all of it. The current small leak is detected as halite precipitation along the Arava stream banks (Supplementary Materials Figure S1; [25]) and as decreasing density along the Arava stream (Supplementary Materials Figure S2; [26]).
The total quantity of the leakage, through the Ye'elim fan delta and beneath the enforced embankments, is calculated accurately by mass balances conducted every year by the Dead Sea Works. However, the leaking brine flow field has never been quantitatively computed using a hydrological model. In this study, the overall leakage over time was recalculated and the underground component of the leakage volumes and pathways were numerically computed, under past, present and future conditions.

Hydrological and Geochemical Characteristics of the Leakage Zone
There is not enough hydrogeological data on the leakage zone (the area between the lake and the ponds) because it was covered by the Dead Sea until the 1980s and, since then, access to this region has been difficult and limited (due to swamp mud and sinkholes). The existing data are as follows: pumping tests in three shallow (~20-30 m) boreholes, which were drilled in the rift-filling permeable gravel (Ye'elim 10, Ye'elim 14 and Hz2; Figure 3), and in one borehole drilled to impermeable clay ( Figure 3). These pumping tests show hydraulic conductivity of about 100 m/day for the gravel infrastructure and about 3 m/day for the clay units [24]. The Na/Cl ratio (equ. ratio) in the Ye'elim 12 borehole (in the gravel unit; Figure 3) was 0.18 in 2011, lower than the typical ratio of 0.2 for the Dead Sea water [27], verifying the hypothesis of leakage from the ponds. Unfortunately, these boreholes are not accessible today.
A Na/Cl ratio of 0.5 was measured by the Dead Sea Works in 2020 in the Rahaf borehole, located 6 km north of the Ye'elim boreholes and also open in the gravel unit [28]. This higher ratio indicates that leakage from the ponds has not reached this site ( Figure 3). A new series of shallow boreholes (~30 m) was recently drilled by the Dead Sea Works for engineering purposes (BHCR boreholes, Figure 3). The hydraulic heads measured in these boreholes show an almost flat groundwater level of about 10 m below the pond level (due to the embankment). These parameters were used for the model calibration (Section 3.3).

Surface Flow Discharge Measurements
The surface water discharge at the marine canal was measured by a SEBA Hydrometrie calibrated propeller-type current meter (Supplementary Materials Figure S3B).

Surface Flow Discharge Measurements
The surface water discharge at the marine canal was measured by a SEBA Hydrometrie calibrated propeller-type current meter (Supplementary Materials Figure S3B). The measurements were conducted in a convenient section that was built in the canal. The whole flow section was divided into small subsections of 0.3-0.4 m width, in which the water velocity was measured at three depths (0.2, 0.4 and 0.8 m from the canal floor). All of the velocity measurements were integrated and multiplied by the section area to get the discharge. The measurement of 24 Mart 2020 was conducted from the canal bank using a rod. The other measurements were conducted using a crane (Supplementary Materials Figure S3A), and are, therefore, more accurate. The measurements were conducted by Yama Company.

Mass Balance of the Northern Evaporation Ponds
A mass balance of the brine in the northern evaporation ponds was calculated using independent data from the Dead Sea Works. The balance was based on measured inflows and measured and calculated outflows. Since the pond volume has remained stable throughout the years, the inflow and outflow are expected to be equal. The inflows consisted of the water pumped from the Dead Sea measured at station P88 ( Figure 1D), the rainfall on the ponds measured by a nearby gauge (in Sedom), and a component of fresh groundwater flow from the west that is assumed to be negligible. The outflows included the evaporated water (calculated by the Dead Sea Works), halite precipitation on the pond floor (measured) and pumping from the northern to southern evaporation ponds at stations P33 and P11 (also measured) ( Figure 1D). Because each component has a different density, we have conducted mass balance, not volume balance, calculations. The volume of each component (except the halite precipitation that is given as a mass) was multiplied by its density (all densities were measured). It was then possible to calculate the mass balance. The gap between the inflows and the outflows was the leakage (Equation (1)): The data for the mass balance calculations for the period 2003-2018 were submitted by the Dead Sea Works. It should be noted that the Dead Sea Works calculates separate mass balances for each of the chemical elements, for production and economic purposes; thus, they are reliable. Nevertheless, all values include some errors or uncertainties.

Flow and Transport Numerical Modeling
A three-dimensional numerical groundwater flow and solute transport model of the aquifer between the evaporation ponds and the Dead Sea was introduced and calibrated. This model was based on a previous model [21,29] introduced for the entire Eastern Mountain Aquifer (EMA) of Israel, the eastern boundary of which is formed by the Dead Sea rift. The results of the previous model were used to define the boundary conditions for the current model, which was included within the area of the previous model ( Figure 4A,C). The model's eastern boundary extended along 10 km of the Dead Sea's southwest coast, the Arava stream and 10 km of the evaporation ponds on the northwest coast. The western boundary was set along an arbitrary line at a distance of~15 km from these coasts. The northern and southern boundaries were two parallel straight lines connecting the eastern and western boundaries ( Figure 4A). The model extended to a depth of 100 m from the Dead Sea coast (the lowest elevation). The sediments within the model were roughly divided into three lithological units: (1) carbonates of the Judea group, having relatively low permeability on the western side of the rift marginal fault, (2) high-permeability rift-filling gravel, and (3) low-permeability rift-filling clay ( Figure 4B). The numerical groundwater flow and solute transport model was introduced by the FEFLOW code [31], using the finite element method for solving the groundwater flow equation (Equation (2)) and the solute transport equation (Equation (3)). The numerical groundwater flow and solute transport model was introduced by the FEFLOW code [31], using the finite element method for solving the groundwater flow equation (Equation (2)) and the solute transport equation (Equation (3)).  [32]) and constant prescribed hydraulic head on the western boundary, with the hydraulic heads derived from the EMA model ( Figure 4C); saline water (c = 342,000 mg/l [33]) and constant prescribed hydraulic head on the eastern boundary. The eastern boundary was divided into three parts, with a suitable hydraulic head set for each: the evaporation ponds, the Arava stream and the Dead Sea ( Figure 5). The Arava stream elevation is well known in the northeast corner of the evaporation ponds [14] and a linear slope was assumed between this corner and the Dead Sea.
The hydraulic head boundary condition along the saline lakes is a saltwater boundary condition. Namely, the equivalent freshwater hydraulic head increases with depth according to the density ratio of the fresh and saline waters (Equations (4) and (5)): The flow model parameters, namely the hydraulic conductivities, were calibrated under the current levels of the Dead Sea (434 mbsl) and evaporation ponds (390 mbsl). The calibration was conducted by fitting the model results to three data systems: observed hydraulic heads in wells, the calculated mass balance of the northern evaporation ponds and pumping tests in some wells. The hydraulic conductivities were adjusted to four units: (1) the Judea group carbonates (this value was derived from the EMA regional model); (2) rift-filling gravel; (3) rift-filling clay; (4) the enforced pond embankment ( Figure 4B). Units 1-3 were calibrated using the current conditions (2020), whereas unit 4 was calibrated according to the reconstructed ponds' mass balance, namely, the leakage decrease due to the enforcing of the pond's embankment in 2012 (Section 4.2). Longitudinal and transverse dispersivities of 200 and 20 m, respectively, were chosen for all of the four units, which are reasonable for the scale involved [32,34].
The calibrated model was used to reconstruct the leakage from the evaporation ponds to the Dead Sea, starting from 1980 (when the Dead Sea was divided into the northern basin and the southern evaporation ponds) until the present (2020), and to predict the leakage in the next 20 years. For the reconstruction and the prediction, the model was run in several period steps: 1981-1990, 1991-2000, 2001-2007, 2008-

Premodeling: Evaporation Pond Mass Balance and Marine Canal Discharge Monitoring
Here, we present the results of the preliminary process, where we defined the boundary conditions for running the numerical model. These results refer to the total brine mass that has leaked from the ponds over the past 40 years (since the splitting of the Dead Sea into two basins), and to the partial brine mass that emerged as a spring north of the pond embankment and was transported through the marine canal to the Dead Sea. The total mass leakage minus the leakage flowing in the marine canal gives the brine mass that leaked from the ponds and was transported in the aquifer, and for which the flow and transport numerical model was introduced and calibrated.
The leakage amounts from the northern salt ponds ( Figure 1D) were obtained from the Dead Sea Works for the early years , but were accurately recalculated for the later years (2003-2018). The leakage volumes are shown in Figure 6 and the mass balance calculations are detailed in Supplementary Materials Table S2. The general trend was an increase in leakage over time. Nevertheless, a large increase occurred in 2008-2012, when the leakage volume "jumped" from 40-60 million cubic meter (mcm) per year to ~130 mcm/year, due to technical failures in the embankments, and a large decrease occurred in 2013, from ~130 mcm/year to ~65 mcm/year, due to the embankments' enforcement. It should be noted that these occurrences well illustrate the relationship between industrial activity and natural forces, namely the formation and elimination of threats to the man-nature coexistence.
As described in Section 3.2, all of the balance components included some uncertainty and errors in both the calculated (evaporation and halite precipitation) and measured (pumping) components. The error of the fluxes at the pumping stations is about 1.5%; the error of the brine densities varies between 0.7-2.2%; the error of the halite mass precipitation is 1.5%; all are relatively small. However, the evaporation rate error is relatively large; specifically, 12-30%. Therefore, for the sake of numerical modeling, we divided the leakage period into steps, each with an approximated constant leakage value (Figure 6, Sup-

Premodeling: Evaporation Pond Mass Balance and Marine Canal Discharge Monitoring
Here, we present the results of the preliminary process, where we defined the boundary conditions for running the numerical model. These results refer to the total brine mass that has leaked from the ponds over the past 40 years (since the splitting of the Dead Sea into two basins), and to the partial brine mass that emerged as a spring north of the pond embankment and was transported through the marine canal to the Dead Sea. The total mass leakage minus the leakage flowing in the marine canal gives the brine mass that leaked from the ponds and was transported in the aquifer, and for which the flow and transport numerical model was introduced and calibrated.
The leakage amounts from the northern salt ponds ( Figure 1D) were obtained from the Dead Sea Works for the early years (1980-2002), but were accurately recalculated for the later years (2003-2018). The leakage volumes are shown in Figure 6 and the mass balance calculations are detailed in Supplementary Materials Table S2. The general trend was an increase in leakage over time. Nevertheless, a large increase occurred in 2008-2012, when the leakage volume "jumped" from 40-60 million cubic meter (mcm) per year tõ 130 mcm/year, due to technical failures in the embankments, and a large decrease occurred in 2013, from~130 mcm/year to~65 mcm/year, due to the embankments' enforcement. It should be noted that these occurrences well illustrate the relationship between industrial activity and natural forces, namely the formation and elimination of threats to the mannature coexistence.
As described in Section 3.2, all of the balance components included some uncertainty and errors in both the calculated (evaporation and halite precipitation) and measured (pumping) components. The error of the fluxes at the pumping stations is about 1.5%; the error of the brine densities varies between 0.7-2.2%; the error of the halite mass precip-itation is 1.5%; all are relatively small. However, the evaporation rate error is relatively large; specifically, 12-30%. Therefore, for the sake of numerical modeling, we divided the leakage period into steps, each with an approximated constant leakage value ( Figure 6, Supplementary Materials Table S3). The constant value of leakage for each period step is obviously an estimation, but it is better for the numerical simulation compared to the yearly fluctuations.
The surface discharge rate in the marine canal was monitored for the first time in 2016 and then three times in 2020 ( Table 1). The momentary discharge in 2016 was equal to 39 mcm/year (assuming constant yearly discharge). In 2020, the momentary discharge was found to be in the range of 82-58 mcm/year. As mentioned in Section 3.1, the March 2020 measurement was less accurate. Thus, the 2020 leak was estimated to be~60 mcm/year.
As expected, the overall leakage from the ponds increased with time, as did the discharge to the marine canal. For the current conditions, on which the numerical calibration was conducted, namely the 2017-2020 period step, a value of 55 mcm/year was assumed for the discharge to the marine canal out of the overall leakage of 80 mcm/year ( Figure 6, Supplementary Materials Table S3). Therefore, 25 mcm/year leaks from the ponds and continues in the aquifer toward the Dead Sea as groundwater.  Table S3). The constant value of leakage for each period step is obviously an estimation, but it is better for the numerical simulation compared to the yearly fluctuations.
The surface discharge rate in the marine canal was monitored for the first time in 2016 and then three times in 2020 ( Table 1). The momentary discharge in 2016 was equal to 39 mcm/year (assuming constant yearly discharge). In 2020, the momentary discharge was found to be in the range of 82-58 mcm/year. As mentioned in Section 3.1, the March 2020 measurement was less accurate. Thus, the 2020 leak was estimated to be ~60 mcm/year.
As expected, the overall leakage from the ponds increased with time, as did the discharge to the marine canal. For the current conditions, on which the numerical calibration was conducted, namely the 2017-2020 period step, a value of 55 mcm/year was assumed for the discharge to the marine canal out of the overall leakage of 80 mcm/year ( Figure 6, Supplementary Materials Table S3). Therefore, 25 mcm/year leaks from the ponds and continues in the aquifer toward the Dead Sea as groundwater.  Table S3). The dashed black line represents the time-step approximations used for the numerical model.

Model Calibration and Initial Conditions
The calibrated hydraulic conductivities for the different geological units were as follows: 0.06 m/day for the Judea Group carbonates (derived from the whole EMA model, Levy et al., 2019Levy et al., , 2020; 70 m/day for the rift-filling gravel; 4.5 m/day for the rift-filling clay; 0.03 m/day for the compacted embankment ( Figure 7A). These parameters gave the best fit between observed and calculated values related to the following four sets of parameters: (1) 25 mcm/year for the net subsurface leakage from the ponds (in addition to  Table S3). The dashed black line represents the time-step approximations used for the numerical model.

Model Calibration and Initial Conditions
The calibrated hydraulic conductivities for the different geological units were as follows: 0.06 m/day for the Judea Group carbonates (derived from the whole EMA model, Levy et al., 2019Levy et al., , 2020; 70 m/day for the rift-filling gravel; 4.5 m/day for the rift-filling clay; 0.03 m/day for the compacted embankment ( Figure 7A). These parameters gave the best fit between observed and calculated values related to the following four sets of parameters: (1) 25 mcm/year for the net subsurface leakage from the ponds (in addition to the 55 mcm/year at the canal for the 2017-2020 period step); (2) the relationship between the two hydraulic conductivity values of the two rift-filling units adequately fit the pumping tests (3 m/day for the clays and~100 m/day for the gravel [24]); (3) the errors between the observed and calculated hydraulic heads in the boreholes were minimized (Figure 7B), namely mean absolute error (MAE) was 2.4 m; (4) the hydraulic conductivity of the embankment enabled a decrease in the leakage between the 2009-2012 and 2013-2016 period steps bỹ 50 mcm/year. It should be noted that these hydraulic conductivities are the horizontal ones; the vertical hydraulic conductivities were 0.1 of the horizontal ones [23]. The specific storage was 0.0004 [23].
The model parameters that were calibrated for the current conditions were used for a forward simulation, starting in 1980 (prior to the splitting of the single Dead Sea basin into two saline lakes) to the present, and for a prediction of the next 20 years. The simulations for the past years were conducted in the same period steps as those approximated for reconstructing the evaporation pond balance (Section 4.1, Figure 6 and Supplementary Materials Table S3). Each simulation period step adopted constant hydraulic head boundary conditions for the Dead Sea, evaporation ponds and the Arava stream ( Figure 8, Supplementary Materials Table S4).
The initial conditions were defined for the conditions of a single saline lake, with water level of 395 mbsl. The groundwater hydraulic head distribution and the salinity field of these initial conditions are shown in Figure 9. It should be noted that these hydraulic conductivities are the horizontal ones; the vertical hydraulic conductivities were 0.1 of the horizontal ones [23]. The specific storage was 0.0004 [23]. The model parameters that were calibrated for the current conditions were used for a forward simulation, starting in 1980 (prior to the splitting of the single Dead Sea basin into two saline lakes) to the present, and for a prediction of the next 20 years. The simulations for the past years were conducted in the same period steps as those approximated for reconstructing the evaporation pond balance (Section 4.1, Figure 6 and Supplementary Materials Table S3). Each simulation period step adopted constant hydraulic head boundary conditions for the Dead Sea, evaporation ponds and the Arava stream (Figure 8, Supplementary Materials Table S4).
The initial conditions were defined for the conditions of a single saline lake, with water level of 395 mbsl. The groundwater hydraulic head distribution and the salinity field of these initial conditions are shown in Figure 9.

Flow and Transport Simulations between the Two Saline Lakes
The calibrated model was used to simulate the groundwater flow regime in the coastal aquifer between the evaporation ponds and the Dead Sea in the past, from 1980 to 2020, and in the future, from 2021 to 2040. These forward simulations were conducted separately for each period step, taking the results of the previous period step as initial conditions for the next one. Although the "real" boundary conditions are continuously changing, we simulated the process in period steps, each with constant boundary conditions, for two reasons: (1) simplifying the model simulations; (2) simulating in accordance with the mass balance results shown in Figure 6 (Section 4.1).
The results of the flow and transport simulations of all period steps, including maps of hydraulic heads and maps of salinity distribution and particle tracking, are shown in Supplementary Materials Figure S4. The results of three important period steps, namely the first one (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990), the present one (2017-2020) and the last one (2031-2040), are shown in Figure 10. The maps in the upper row show the temporal increase in the hydraulic gradient between the evaporation ponds and the Dead Sea. The maps in the bottom row show that the aquifer near the Dead Sea was saturated with saline water after 1980, because it served as the lake floor prior to its drying, but it was washed thereafter (the washing also was observed in previous study [35]). At the same time, brine from the evaporation ponds leaked northward to the aquifer. According to the particle-tracking analysis, the brine that started to leak from the evaporation ponds 40 years ago still (as of 2020) has not arrived to the Rahaf borehole, as indeed detected in its chemical composition (see Section 2.4). Nevertheless, in the next period step (2021-2040) the leaked brine is expected to reach this borehole.
From 1980-2012, leakage from the evaporation ponds occurred throughout the northern and eastern embankments. After the construction of the underground sealing wall in 2012, most of the leakage stopped. However, at the northwestern corner of the ponds, at the Ye'elim alluvial fan ( Figure 1E), brine leakage continued, and since then it has dramatically increased. The high leakage through the fan is reflected in the velocity vectors, which are several orders of magnitude greater than for its surroundings ( Figure 11).
Based on the numerical simulations and the mass balance, the total leakage to the aquifer was disassembled into five leakage components, three of which were numerically calculated and two estimated. The three calculated components were: through the Ye'elim fan, through the northern embankment, and below the northern embankment. The estimated components were: through the marine canal (measured only in the two last period steps) and through the eastern embankment ( Figure 12, Supplementary Materials Table S5).
The overall leakage accumulates with time mainly because of the increasing hydraulic gradient between the evaporation ponds and the Dead Sea. In the 2009-2012 period step, the leakage showed an extreme increase, mainly by the "through embankments" components, probably due to technical failures in the embankments (personal communication, Dead Sea Works). In 2012, the high leakage through the embankments stopped following their compaction and construction of the underground vertical sealing wall. In the modeling simulations, the hydraulic conductivity of the northern embankment changed accordingly (Figure 7). Since 2012, the leakage through the Ye'elim alluvial fan and below the embankment has been showing a gradual increase ( Figure 12) and is expected to continue growing as long as the hydraulic gradient increases. The Achilles' heel is the Ye'elim alluvial fan where most of the leakage is occurring. Therefore, in the next stage of the "industrial vs. natural flow regime" conflict, this site must be focused on.  Note the difference in the high salinity between (D-F) indicating brine washing with time on the newly exposed land. In (D-F), The black-white lines are particle tracking of brine leakage from the pond's northern embankment progressing toward the Rahaf borehole (location is shown in (E,F)). The dashed rectangular area on map A represents the area that is enlarged in maps (D-F). The dashed rectangular area on map E is the area shown in Figure 11.
Water 2021, 13, x FOR PEER REVIEW 15 of 18 Figure 10. Maps of groundwater level (meters above sea level, (A-C) and salinity distributions (D-F) for three period steps: first, present and future. The head contour difference is 50 m in general (values on yellow background) but 10 m between the ponds and the lake, where the increasing gradient with time is notable (values in white). Note the difference in the high salinity between (D-F) indicating brine washing with time on the newly exposed land. In (D-F), The black-white lines are particle tracking of brine leakage from the pond's northern embankment progressing toward the Rahaf borehole (location is shown in (E,F)). The dashed rectangular area on map A represents the area that is enlarged in maps (D-F). The dashed rectangular area on map E is the area shown in Figure 11. Figure 11. Groundwater velocity vectors (ranked by colors), demonstrating the massive leak through the Ye'elim fan delta (see Figure 10E). Figure 11. Groundwater velocity vectors (ranked by colors), demonstrating the massive leak through the Ye'elim fan delta (see Figure 10E).

Summary and Conclusions
The groundwater flow regime in the coastal aquifer adjacent to the Dead Sea's southern edge is significantly affected by both industrial activity and natural forces. Since 1980, the one unified lake has been divided into two basins: the northern deep Dead Sea, and the evaporation ponds of the Dead Sea Works in the south. With time, the evaporation pond level is increasing and the Dead Sea level is decreasing; thus, the hydraulic gradient between the two is increasing. The embankments of the evaporation ponds were built before the hydraulic gradient formed; they were not designed to withstand the developing loads. As a result, technical failures have developed with time in the embankments, and in the years 2008-2012, there was a large "jump" (more than double) in the leakage volume. Based on continuous monitoring and scientific research, the Dead Sea Works drew the correct conclusions, and compacted the embankments and constructed the vertical sealing wall inside the embankments. Consequently, the groundwater flow regime was completely changed, blocking most of the leakage. Nevertheless, it is not easy to confront nature and to withstand the increasing hydraulic gradient. Therefore, since 2012, the leakage through the Achilles' heel-the Ye'elim alluvial fan, where no seal was placedhas been continuously increasing. Since the hydraulic gradient is expected to increase in the foreseeable future, the brine leak will increase as well, and a new threat to the embankment may arise. This paper demonstrates how the conflict between anthropogenic development and natural forces in the Dead Sea region has evolved with time. There is no other alternative for maintaining the coexistence between man and nature but to continue monitoring, conduct research, and be flexible enough to respond in time to the unexpected. The coexistence between settlement, agriculture, industry and tourism on one hand, and nature on

Summary and Conclusions
The groundwater flow regime in the coastal aquifer adjacent to the Dead Sea's southern edge is significantly affected by both industrial activity and natural forces. Since 1980, the one unified lake has been divided into two basins: the northern deep Dead Sea, and the evaporation ponds of the Dead Sea Works in the south. With time, the evaporation pond level is increasing and the Dead Sea level is decreasing; thus, the hydraulic gradient between the two is increasing. The embankments of the evaporation ponds were built before the hydraulic gradient formed; they were not designed to withstand the developing loads. As a result, technical failures have developed with time in the embankments, and in the years 2008-2012, there was a large "jump" (more than double) in the leakage volume. Based on continuous monitoring and scientific research, the Dead Sea Works drew the correct conclusions, and compacted the embankments and constructed the vertical sealing wall inside the embankments. Consequently, the groundwater flow regime was completely changed, blocking most of the leakage. Nevertheless, it is not easy to confront nature and to withstand the increasing hydraulic gradient. Therefore, since 2012, the leakage through the Achilles' heel-the Ye'elim alluvial fan, where no seal was placed-has been continuously increasing. Since the hydraulic gradient is expected to increase in the foreseeable future, the brine leak will increase as well, and a new threat to the embankment may arise. This paper demonstrates how the conflict between anthropogenic development and natural forces in the Dead Sea region has evolved with time. There is no other alternative for maintaining the coexistence between man and nature but to continue monitoring, conduct research, and be flexible enough to respond in time to the unexpected. The coexistence between settlement, agriculture, industry and tourism on one hand, and nature on the other, could quickly succumb to disaster. Loss of control over such processes is, unfortunately, a normal scenario, which raises concerns about the fate of the area.
On the macro-scale, overexploitation of fresh water in the Dead Sea catchment basin has accelerated the lake's drying. This will probably stop in the near future, by supplying sea water through a sea canal, or in the far future, by natural stabilization of the lake level at about 100 m below the present level. Stream incisions and sinkhole formation threaten the infrastructure around the lake, but these processes are still under control. As seen in the present study, the groundwater flow regime has changed back and forth due to actions taken by the industry, and collapse of the evaporation pond embankments has been avoided.  Figure S2. Evidence of brine leakage along the Arava stream based on changes in water density (based on Sheffer, 2013 [24]). (A) Map of the sampling points, marked 1 to 7 from south to north. (B) Measured water densities show an increase between points 1 and 4 due to the entrance of end brines from the bromine factory (BEB) and from the Israeli (IEB) and Jordanian (JEB) potash factories, and a decrease between points 4 and 7 due to leakage from the evaporation ponds. Figure S3. (A) Measuring water discharge in the marine canal at several points throughout the flow section. (B) Using a propeller-type current meter. Figure S4. Maps of groundwater level (meters above sea level, top row) and salinity distributions (bottom row) for each of the time steps (past reconstruction and future predictions). The head contour difference is 50 m in general but 10 m be-tween the ponds and the lake, where the increasing gradient with time is notable. Note the brine washing with time on the newly exposed land. The black-white lines are particle tracking of brine leakage from the pond's northern embankment progressing toward Rahaf borehole. The dashed rectangular area on the left top map represents the enlarged area in the bottom row maps. The dashed rectangular area on the 2017-2020 period map (bottom row) is the area shown in Figure 11 in the main text. Table S1. Geochemical parameters of the brines at several measurement points (based on Farber, 2016 [30]). The similarity between the first four values and the difference between these and the last two values indicate leakage from the evapo-ration ponds. Table S2. Mass balance calculations at the northern evaporation ponds for 2003-2018. Table S3. Representative leakage volumes from the evaporation ponds for the different time periods (in accordance with Figure 6 in the main text). Table  S4. Constant head boundary conditions and total representative leakage volume for the different time periods, serving as input for the numerical modeling (in accordance with Figure 8 in the main text). Table S5. Reconstructed and predicted volumes of brine leaked from the evaporation ponds through different routes at each time period (in accordance with Figure 12 in the main text). Colors match those in Figure 12 in the main text.