Impacts of Large-Scale Groundwater Exploitation Based on Long-Term Evolution of Hydraulic Heads in Dhaka City, Bangladesh

: Dhaka city has emerged as the fastest-growing megacity, having more than 20 million in-habitants, with a growth rate of 3.62%. Unplanned and rapid urbanization, coupled with exponential population growth, has signiﬁcantly altered the groundwater dynamics in Dhaka city. This study concentrates on the evolution of long-term piezometric heads of the Upper Dupi Tila aquifer (UDA) and the Middle Dupi Tila aquifer (MDA) based on long-term hydrographs, piezometric maps and synthetic graphical overviews of piezometric trends. Due to over-exploitation, the piezometric level (PL) has declined deeper than − 85 and − 65 m PWD (Public Works Department reference datum) in UDA and MDA, respectively. The highest rate of decline was observed in the south-central to southeastern parts of the city both in UDA (4.0 m/year) and MDA (5.74 m/year). The results clearly show that the rates of decline in PL vary from 2.25 to 5.74 m/year in both aquifers of the city, and urban expansion has greatly affected the shape and extent of the depression cone over the past four decades. The magnitudes of the depression cones in both aquifers seem to pose a considerable threat to groundwater resources, indicating that the current exploitation is not sustainable at all.


Introduction
Groundwater overexploitation due to urban expansion and population growth is a critical issue in Dhaka city [1], as it is in many other developing or developed cities in the world such as in Beijing [2,3], Central India [4], the Indus Basin [5], Northwest India [6], Northern Chile [7], California's Central Valley in the U.S. [8], the U.S. High Plains [9], the Middle East [10], the Hai River plain in China [11] and Belgium [12], among others. A continuous decline in groundwater levels (GWLs) has been observed in many places in the past half-century [13,14] but there is no efficient and adequate approach to groundwater sustainability [13]. Mexico, Melbourne, Jakarta, São Paulo, Central Iran, Chennai and major cities in India have been critically affected by massive groundwater abstraction and have already experienced land subsidence and the drying up of rivers and wetlands [6,[15][16][17][18][19][20][21][22][23][24].
Dhaka, the capital of Bangladesh, is densely populated and is the largest city in the country. It is the world's 6th most populated city [25]. It started to expand exponentially after the independence of Bangladesh in 1971. This exponential population growth accelerated groundwater abstraction in order to meet the water demand. In 2018, 78% of the supply water came from underground sources and 22% from treating surface water [26].

History of the Expansion of Dhaka City
The pattern of expansion and urbanization of Dhaka city has mainly been shaped by the landscape and elevation of the land [55]. The historic evolution of Dhaka city generally followed two patterns: old Dhaka or the historic core, and new Dhaka or the northern expansion [56]. The city reached its present status, with its fast-growing population and increasing size, through a series of progressive changes. The history, growth and development of Dhaka can be categorized into six periods: the pre-Mughal period (1205-1610), the Mughal period (1620-1757), the East India Company period (1758-1858), the British colonial period , the Pakistan period  and the Bangladesh period (from 1971 onwards) [45,57] (Figure 2).
In the pre-Mughal period, about 400 years ago, Dhaka was limited to an area of 1.5 km 2 around the river Buriganga [45,57]. According to Chowdhury and Faruqui (1989) [58], Dhaka attracted European traders-the Portuguese, the Dutch, the English, the French and the Armenians-because of its commercial importance and they established their trading posts in the 18th century. In the British period, Dhaka expanded to 22 km 2 and in the Pakistan period to 50 km 2 . The city extended and developed to the west and the north during the Mughal and British periods. Dhaka city was established as the capital of East Pakistan in the Pakistan period. This city progressively became the central point of attraction and development towards the north. A large number of people from other districts started migrating to Dhaka from the Pakistan period onward [59].
After the independence of Bangladesh in 1971, the main urban development and rapid expansion started in Dhaka. The development expanded from the old area towards the west and north in the period between 1972 and 1980 [60,61]. From 1980 to 1990, the city expanded up to Mirpur, Airport and further north to the banks of the Turag river. The expansion to the east and northwest was limited due to the presence of lowlands and depressions, and therefore major development occurred along the north-south corridor

History of the Expansion of Dhaka City
The pattern of expansion and urbanization of Dhaka city has mainly been shaped by the landscape and elevation of the land [55]. The historic evolution of Dhaka city generally followed two patterns: old Dhaka or the historic core, and new Dhaka or the northern expansion [56]. The city reached its present status, with its fast-growing population and increasing size, through a series of progressive changes. The history, growth and development of Dhaka can be categorized into six periods: the pre-Mughal period (1205-1610), the Mughal period (1620-1757), the East India Company period (1758-1858), the British colonial period , the Pakistan period  and the Bangladesh period (from 1971 onwards) [45,57] (Figure 2).
In the pre-Mughal period, about 400 years ago, Dhaka was limited to an area of 1.5 km 2 around the river Buriganga [45,57]. According to Chowdhury and Faruqui (1989) [58], Dhaka attracted European traders-the Portuguese, the Dutch, the English, the French and the Armenians-because of its commercial importance and they established their trading posts in the 18th century. In the British period, Dhaka expanded to 22 km 2 and in the Pakistan period to 50 km 2 . The city extended and developed to the west and the north during the Mughal and British periods. Dhaka city was established as the capital of East Pakistan in the Pakistan period. This city progressively became the central point of attraction and development towards the north. A large number of people from other districts started migrating to Dhaka from the Pakistan period onward [59].
After the independence of Bangladesh in 1971, the main urban development and rapid expansion started in Dhaka. The development expanded from the old area towards the west and north in the period between 1972 and 1980 [60,61]. From 1980 to 1990, the city expanded up to Mirpur, Airport and further north to the banks of the Turag river. The expansion to the east and northwest was limited due to the presence of lowlands and depressions, and therefore major development occurred along the north-south corridor of Dhaka city [55]. Initially, the progress of development was limited to the available highlands. But later on, development also occurred in the low-lying areas. The major land filling areas were mainly in the eastern and south-eastern parts of the city [60].
Until 1990, the main city was restricted between the Balu and Turag rivers in the east and west, respectively, due to the area beyond this limit being low, swampy and floodaffected [62]. In 1996, the Dhaka Metropolitan Development Plan (DMDP) administered by Rajdhani Unnayan Kartripakkha (RAJUK) provided a long-term strategic plan up to 2015 for the development of a greater Dhaka, which includes Dhaka city and its surrounding six municipalities (1528 km 2 ) [63]. In the 1990s and 2000s, the new Dhaka expanded north to its present boundaries. Between 2000 and 2015, the city expanded continuously in all directions and occupied all the vacant spaces [60,62,64]. The city did not yet have any definite and well-planned outline for homogenous and cohesive growth; rather, different private organizations occupied and filled up major ditches, swamps and marshes in an unplanned manner [65,66]. of Dhaka city [55]. Initially, the progress of development was limited to the available highlands. But later on, development also occurred in the low-lying areas. The major land filling areas were mainly in the eastern and south-eastern parts of the city [60]. Until 1990, the main city was restricted between the Balu and Turag rivers in the east and west, respectively, due to the area beyond this limit being low, swampy and floodaffected [62]. In 1996, the Dhaka Metropolitan Development Plan (DMDP) administered by Rajdhani Unnayan Kartripakkha (RAJUK) provided a long-term strategic plan up to 2015 for the development of a greater Dhaka, which includes Dhaka city and its surrounding six municipalities (1528 km 2 ) [63]. In the 1990s and 2000s, the new Dhaka expanded north to its present boundaries. Between 2000 and 2015, the city expanded continuously in all directions and occupied all the vacant spaces [60,62,64]. The city did not yet have any definite and well-planned outline for homogenous and cohesive growth; rather, different private organizations occupied and filled up major ditches, swamps and marshes in an unplanned manner [65,66].  [60,67]). Figure 2. The historical development of Dhaka City at different stages of its growth (modified from [60,67]).

Climate
Dhaka city has a humid tropical monsoon-type climate with an annual average rainfall of 1955 mm. About 80% of the annual average rainfall occurs between the months of May to October. The average potential evaporation is 836 mm per annum. The average relative humidity ranges from 55% to 79%. The maximum monthly average temperature is 30 • C in June and the minimum monthly average temperature is 19 • C in January (Figure 3). , x FOR PEER REVIEW 5 of 29

Climate
Dhaka city has a humid tropical monsoon-type climate with an annual average rainfall of 1955 mm. About 80% of the annual average rainfall occurs between the months of May to October. The average potential evaporation is 836 mm per annum. The average relative humidity ranges from 55% to 79%. The maximum monthly average temperature is 30 °C in June and the minimum monthly average temperature is 19 °C in January (Figure 3).

Geomorphology, Surface Geology and Topography
Dhaka is situated in the southern extremity of the Madhupur tract, a north-south elongated Pleistocene terrace [54] having an undulating topography, elevated from the Recent floodplain [53]. Fifteen geomorphic units have been noted, but centrally most of the parts are covered by the Madhupur terrace unit (Figure 4a). Most of the north-eastern, middle eastern and northwestern margins remarkably cover swamp or depression units. All low-lying geomorphic units, e.g., alluvial valley, active flood plain, natural levee, point bar, abandoned channel bed, swamp and depression area, have already been filled up due to rapid urban growth.
The surface geology mainly follows the geomorphic expression of the area ( Figure  4b) and is classified into four broad categories: Madhupur Clay residuum, Alluvium silt and clay, Alluvium silt, and Marshy peat and clay. From a geological and geomorphological perspective, the city mainly covers two broad geomorphological units. One is high lands, which correspond to Madhupur terrace, and the other is lowlands, which correspond to alluvial deposits, mainly comprising swamp or depression, back swamp, channel bar, flood plain and valley fill units [68,69]. The relatively flat topography also perfectly resembles the geomorphic features of the area. The western part is marked by an abrupt change in elevation, whereas relatively gentle slopes are maintained in the eastern part. The surface elevation of gullies, valleys and bars is somewhat higher (>9 m) in the southwestern and northeastern parts than in the floodplain, depressions and abandoned channels, but the undulating topography was mainly created or developed by the rivers Buriganga, Balu and Turag. The general slope of the land surface is from north to south, more specifically towards the southeast, in which ground surfaces merge with the Recent floodplain of the rivers Buriganga, Balu and Shitalakhay. The elevation differences are remarkably reflected by distinct landforms: a high elevation is observed in the Madhupur terraces region, normally ranging from 15 m to >20 m above mean sea level (amsl), and a

Geomorphology, Surface Geology and Topography
Dhaka is situated in the southern extremity of the Madhupur tract, a north-south elongated Pleistocene terrace [54] having an undulating topography, elevated from the Recent floodplain [53]. Fifteen geomorphic units have been noted, but centrally most of the parts are covered by the Madhupur terrace unit (Figure 4a). Most of the north-eastern, middle eastern and northwestern margins remarkably cover swamp or depression units. All low-lying geomorphic units, e.g., alluvial valley, active flood plain, natural levee, point bar, abandoned channel bed, swamp and depression area, have already been filled up due to rapid urban growth.
The surface geology mainly follows the geomorphic expression of the area ( Figure 4b) and is classified into four broad categories: Madhupur Clay residuum, Alluvium silt and clay, Alluvium silt, and Marshy peat and clay. From a geological and geomorphological perspective, the city mainly covers two broad geomorphological units. One is high lands, which correspond to Madhupur terrace, and the other is lowlands, which correspond to alluvial deposits, mainly comprising swamp or depression, back swamp, channel bar, flood plain and valley fill units [68,69]. The relatively flat topography also perfectly resembles the geomorphic features of the area. The western part is marked by an abrupt change in elevation, whereas relatively gentle slopes are maintained in the eastern part. The surface elevation of gullies, valleys and bars is somewhat higher (>9 m) in the southwestern and northeastern parts than in the floodplain, depressions and abandoned channels, but the undulating topography was mainly created or developed by the rivers Buriganga, Balu and Turag. The general slope of the land surface is from north to south, more specifically towards the southeast, in which ground surfaces merge with the Recent floodplain of the rivers Buriganga, Balu and Shitalakhay. The elevation differences are remarkably reflected by distinct landforms: a high elevation is observed in the Madhupur terraces region, normally ranging from 15 m to >20 m above mean sea level (amsl), and a low elevation (<5 m amsl) is observed in the floodplain, depressions and abandoned channels ( Figure 5a).

Geology
Dhaka city is distinguished by a 400-to 500-m-thick unconsolidated sequence of fluvio-deltaic sediments which are overlain by Madhupur Clay or floodplain materials, which are 5 to 25 m thick [1]. This underlying unconsolidated layer acts as the main aquifer and is known to be a part of the Dupi Tila Formation [70]. The Pleistocene alluvial terraces are very weathered, oxidized, more compacted than the Recent floodplains and are represented as paleo-floodplains of the previous Ganges Brahmaputra Meghna (GBM) system [54,71]. The eroded part of the Madhupur Formation and the subsided part of the floodplain and depressed areas are believed to have been filled by the Bashabo Formation. The Holocene Bashabo Formation mainly consists of sand, silt and clay and serves as a localized shallow aquifer on top of the Madhupur Formation, merging with the floodplain deposits in the southeast part of Dhaka city [53,72]. The Madhupur Tract around Dhaka city is characterized by a large number of faults [70,73] that significantly control the local rivers and possibly affect the aquifer and river interaction. The rivers and stream systems flowing along this fault system possibly cut Dupi Tila sands [1,70]. A series of six secondary echelon faults occur in the western margin of the Madhupur Tract. A coupled force is acting that increases deep slip as well as strike-slip movement along the main Madhupur Fault. For this reason, the uplifted Madhupur Tract is moving towards the north, whereas the downthrown Jamuna floodplain is sliding towards the south [74]. It is expected that the Girujan Clay Formation is underlying the Dupi Tila Formation in Dhaka city. The Girujan Clay Formation was neither exposed nor penetrated by drilling but is assumed based on its extensive occurrences in a different place in Bangladesh [41]. The stratigraphic units of the study area, from oldest to youngest, are the Girujan Clay (Pliocene), Dupi Tila Sands (Plio-Pleistocene), Madhupur Clay (Pleistocene), Bashabo Formations (Holocene) and Alluvium (Recent).

Geology
Dhaka city is distinguished by a 400-to 500-m-thick unconsolidated sequence of fluvio-deltaic sediments which are overlain by Madhupur Clay or floodplain materials, which are 5 to 25 m thick [1]. This underlying unconsolidated layer acts as the main aquifer and is known to be a part of the Dupi Tila Formation [70]. The Pleistocene alluvial terraces are very weathered, oxidized, more compacted than the Recent floodplains and are represented as paleo-floodplains of the previous Ganges Brahmaputra Meghna (GBM) system [54,71]. The eroded part of the Madhupur Formation and the subsided part of the floodplain and depressed areas are believed to have been filled by the Bashabo Formation. The Holocene Bashabo Formation mainly consists of sand, silt and clay and serves as a localized shallow aquifer on top of the Madhupur Formation, merging with the floodplain deposits in the southeast part of Dhaka city [53,72]. The Madhupur Tract around Dhaka city is characterized by a large number of faults [70,73] that significantly control the local rivers and possibly affect the aquifer and river interaction. The rivers and stream systems flowing along this fault system possibly cut Dupi Tila sands [1,70]. A series of six secondary echelon faults occur in the western margin of the Madhupur Tract. A coupled force is acting that increases deep slip as well as strike-slip movement along the main Madhupur Fault. For this reason, the uplifted Madhupur Tract is moving towards the north, whereas the downthrown Jamuna floodplain is sliding towards the south [74]. It is expected that the Girujan Clay Formation is underlying the Dupi Tila Formation in Dhaka city. The Girujan Clay Formation was neither exposed nor penetrated by drilling but is assumed based on its extensive occurrences in a different place in Bangladesh [41]. The stratigraphic units of the study area, from oldest to youngest, are the Girujan Clay (Pliocene), Dupi Tila Sands (Plio-Pleistocene), Madhupur Clay (Pleistocene), Bashabo Formations (Holocene) and Alluvium (Recent).

Hydrostratigraphy and Aquifer Delineation
The study area has been geologically and hydrogeologically investigated by many researchers [1,[41][42][43]45,70,72,[75][76][77] but the geology, stratigraphy and hydrostratigraphy are not well defined because of the limited depth investigated. Almost all the investigations have been limited to the upper part of the Dupi Tila Formation (<170 m). In this study, an attempt has been made to delineate the more detailed hydrostratigraphy and aquifer geometry at a greater depth (up to 463 m). The thickness of the hydrostratigraphic units varies laterally and vertically because of surface topography and changes in the depositional conditions [71], that is, the heterogeneity of the layers [71,78]. A hydrostratigraphic reconstruction was performed based on lithologs from the Bangladesh Water Development Board (BWDB), the Department of Public Health Engineering (DPHE) and the Dhaka Water Supply and Sewerage Authority (DWASA). Golden Software's Strater ® 5 was used for the interpretation of the lithology and for generating cross-sections. First, the location of the borehole logs was mapped using Strater ® 5. Suitable lines were chosen for the lithological cross sections. The borehole logs were then plotted along the cross-section line. The aquifers and aquitards were delineated manually, based on lithology, using the drawing option in Strater ® 5. Similarly, the PL in each borehole log along the cross-section was inserted and then connected manually to show the PL. The generalized hydrostratigraphy was established based on lithology and, chronologically, from the surface up to the explored depth (463 m), the units are: localized shallow aquifer, aquitard, Upper Dupi Tila Aquifer (UDA), aquitard, Middle Dupi Tila Aquifer (MDA), aquitard, Lower Dupi Tila Aquifer (LDA) and aquitard.
The top aquitard unit is composed of silty clay, sandy silt and clay. The depth to base varies from 5 to 64 m and the average thickness is 30 m (Figure 5b-d). The UDA is overlain by this top aquitard and underlain by another aquitard. This aquifer is semi-confined. The UDA consists of unconsolidated sands of varying size: very fine sand to fine sand, fine sand to medium sand, and medium sand to coarse sand occasionally with gravel ( Figure 5b-d). The depth to base of UDA is generally 90 to 170 m and, on average, the bottom depth is 143 m. The average thickness is 98 m. In most parts of the city, the transmissivity values of UDA range from 400 to 500 m 2 /d and the storativity values range from 1.5 × 10 −6 to 9.5 × 10 −5 but are mainly around 2 × 10 −5 [79]. There are laterally discontinuous clay or silty clay layers in the Dupi Tila Formation, which are also regarded as local aquitards within UDA, MDA and LDA. The second aquitard is present between UDA and MDA, and is composed of sandy silt, silty clay and clay, with an average thickness of 32 m. MDA is mainly composed of fine-to-medium sand and medium-to-coarse sand, occasionally with gravel. The depth to its base ranges from 190 to 295 m, with an average thickness of 86 m. In MDA, the average transmissivity, storage coefficient and hydraulic conductivity in Mirpur are reported to be 1756.14 m 2 /d, 0.0029 m/d and 22.16 m/d, respectively, whereas in Tejgaon, they are 1854.94 m 2 /d, 0.0019 m/d and 32.03 m/d, respectively [80]. A very thick aquitard underlies MDA and the depth to its base ranges from 378 to 384 m (in the few boreholes that crossed this layer). The LDA is predominantly composed of fine-to-medium sand and medium-to-coarse sand. The depth to its base varies from 420 to 426 m and the average thickness is 40 m (in the few available boreholes). Very few lithologs are available for the purpose of illustrating the LDA (Figure 5d).

Abstraction Scenario for an Expanding Dhaka City
The groundwater exploitation scenario is directly linked with the population of the area, urbanization and city expansion. In Dhaka city, groundwater extraction has been reported since the 16th century: Guru Nanak sunk a dug well due to the scarcity of water [82]. In the 1980s, dug wells were also used as a source of drinking water, and still today some dug wells are found in the alluvial cover on top of the Madhupur Clay. In 1874, Nabab Khaja Abdul Ghani first established a water treatment plant in the Chadnighat area near the bank of the river Buriganga to supply water to Dhaka city [45]. The DPHE introduced groundwater development for a public supply and drilled the first borehole in Dhaka in 1941 [45,57]. Then the use of groundwater increased gradually over time as a source of water production [45]. Since 1963, DWASA has been continuously expanding its coverage area with updated services [45,57]. Figure 6 shows the development of production wells, along with water demand and production and population with time. The independence of Bangladesh in 1971 prompted the advance and development of the groundwater supply [30]. The city had a population of less than two million and groundwater abstraction was less than 50 × 10 6 m 3 /year during that time, and only 47 DWASA deep tube wells (DTWs) were operating [26]. From 1970 to 1985, almost all the DTWs were installed in UDA in the vicinity of the rivers Buriganga and Turag, whereas later on, well construction spread to the center and the northern part of the city. In 2000, the population increased to almost 10 million. DWASA installed a total of 308 DTWs with 413 × 10 6 m 3 /year production capacity, but demand was 548 × 10 6 m 3 /year. In addition, more than 900 private deep wells (with a minimum of 0.057 m 3 /s discharge) were also Water 2021, 13, 1357 9 of 28 operating in the UDA [46]. Since 2004, DWASA has drilled wells to larger depths (>170 m) in MDA and LDA [80] because of the huge drop in the PL of UDA. In 2013, the supply capacity was raised to 883 × 10 6 m 3 /year through the use of 644 DTWs and four water treatment plants [26]. In 2018, the population increased to 20 million. To meet the water demand, both the production capacity (931 × 10 6 m 3 /year) and the number of DTWs (887) were increased [26].

Methodology: Groundwater Level and Piezometric Data Collection
Water table data and PL data were collected from BWDB and DWASA in order to prepare long-term hydrographs and piezometric maps of UDA and MDA. DWASA assigned the Institute of Water Modelling (IWM) to establish twin monitoring wells in order to observe the PL of both UDA and MDA [83]. DWASA installed 10 twin monitoring wells in 2005, but data recording was not continuous, and they installed four more twin wells in 2016 ( Figure 6). Twin wells are paired wells for monitoring PL in the same area but at different depths, that is, one well in UDA (<170 m depth) and another in MDA (>200 to <290 m).
The water table and PL data were collected from 13 observation wells for the period of 1969 to 2018 from BWDB. Measurements were carried out using manual groundwater measuring devices for 11 observation wells and automatic recording devices for two observation wells (Banani and Lalbagh). The starting dates of measurements by all the piezometers and their depths were not the same. Initially, most of these observation wells were dug wells measuring the water table; then these were subsequently replaced by piezometers. The weekly data from dug wells reflect the water table of that particular area. In some cases, new observation wells were installed in the same or different locations at greater depth, to replace the older wells. In the case of DWASA production wells, the static water level of production wells during well construction was collected. No data were taken from private wells. The PL elevation is calculated with respect to the Public Works Department (PWD) reference datum. The PWD datum is 0.46 m below the MSL datum [84].
In order to understand the trend of the depression cone in Dhaka city, piezometric maps of UDA and MDA were constructed based on data from observation wells, twin wells and DWASA production wells (Figure 7). To construct piezometric maps, weekly PL data were converted to yearly averages for each observation well from 1980 to 2018. The piezometric maps for UDA were generated for every 10 years for the period from 1980 to 2010, and for the years 2015 and 2018. Piezometric maps for MDA were also generated for the years 2005, 2010, 2015, 2016, 2017 and 2018. DWASA has drilled in MDA since 2004 due to the vast decline of the PL in UDA; for this reason, piezometric maps for this aquifer were prepared starting from the year 2005. The PL data were analyzed using Microsoft Excel. The base map for constructing piezometric maps was prepared using QGIS and the manually prepared PL elevation contours were digitized using QGIS.

Long-Term Hydrographs: Analysis of PL Fluctuation
The old part of the city was developed in the vicinity of the Buriganga and Turag rivers and, therefore, most of the observation wells are situated in the south-central and western parts of the city. However, twin wells are situated sparsely throughout the city and show the changes of PL in peripheral and central parts of the city. In Dhaka city, most of the observation wells show a continuously declining trend ( Figure 8).

Motijheel (South-Central Part of Dhaka City)
The Motijheel area became earmarked as the center of regional administration and business activities before 1971 [62,64]. The water table data (dug well, 6.22 m deep) in Motijheel showed seasonal cyclicity with no declining trend from 1979 to 1984. In 1985, the dug well was replaced by a drilled piezometer. Since then, this well showed a continuously decreasing trend; in the first two years (1985 to 1987) PL declined exponentially with an average of 4.8 m per year drawdown, and from 1987 to 1989, the lowering trend was linear (Figure 8a). In 1989, the observed PL was −16.03 m PWD. Groundwater was exploited at this depth before the construction of the piezometer. After 1985, most of the dug wells were replaced by tubewells, and production wells were also drilled. Water was supplied through pipelines from this part of the city. Thus, the first two years show an exponential decline near the pumping area. After 1987, the well construction was extended to the center and northern parts of the city. The rate of pumping after 1987 in the area of Motijheel decreased gradually and, at the same time, abstraction was disseminated over a larger part of the city, followed by the PL decline becoming linear in this part of the city.  (Figure 9d). From 2000 onwards, Banani became an expensive diplomatic area in Dhaka city. Thus, the water demand in this area has been very high since then, and a sharp decline has been observed in this area in both aquifers. In 2018, the PL was −71.24 and −61.5 m PWD for UDA and MDA, respectively. The Nasirabad (southeast), Gandaria (south west) and Goranchatbari (central west) twin wells showed linear declines since 2016 but both aquifers' hydrographs were almost coinciding and parallel to each other. In 2018 in Gandaria, the PL was −48.75 and −48.23 m PWD, whereas in Goranchatbari the PL was −40.12 and −40.32 m PWD for UDA and MDA, respectively (Figure 9e,f)). Most of the twin hydrographs suggest that the declining rate was higher in MDA than in UDA from 2016 to 2018 (Table 1, Figure 9), not only in the central region but also in the peripheral region (Figure 9g  The water table map of the shallow aquifer for 1980 was prepared based on six dug wells with depths from 6.22 to 12.85 m. The highest water table, at 6.0 m PWD, was found in Sabujbagh and the lowest water table, at −0.5 m PWD, was recorded in Mirpur ( Figure 10). The highest water table was observed in the highest topographic area toward the south-central region and lowest water table in the lowest topographic area in the direction of the west-central part of the city. The average fluctuations of the water levels of the surrounding rivers (Buriganga, Turag Tongi Khal, Balu and Shitalakhya) of Dhaka city were from 7.2 to 0.8 m during the wet and dry seasons, respectively, since 1970. The water table and surrounding river water levels are directly correlated to each other. Almost all the dug wells were abandoned or replaced by drilled observation wells after 1984, and this is why the water table map is only available for 1980. Only a few dug wells were monitored, and therefore no isolines were drawn on the map. The water table in shallow aquifers is replenished every year during the monsoon season through the vertical percolation of rainwater [86,87] or floodwater in Bangladesh. In Dhaka city, dug wells are no longer in use because of the scarcity of lifting technology, poor maintenance and inferior water quality [45]. dug wells were abandoned or replaced by drilled observation wells after 1984, and this i why the water table map is only available for 1980. Only a few dug wells were monitored and therefore no isolines were drawn on the map. The water table in shallow aquifers i replenished every year during the monsoon season through the vertical percolation o rainwater [86,87] or floodwater in Bangladesh. In Dhaka city, dug wells are no longer i use because of the scarcity of lifting technology, poor maintenance and inferior wate quality [45].

Spatial and Temporal Distribution of Piezometric Level for UDA
Although groundwater development started for public supply in 1941 in Dhak [45,57,88], due to the unavailability of piezometric data, it was not possible to make pie zometric maps before 1980. The existing extraction of groundwater formed a local cone o depression in 1980. The PL contour of 0 m follows the outline of the city in that period The contours of −1 and −2 m PWD encompass the local depression due to groundwate pumping. Though PL was very close to the ground surface in the 1980s, a depression o −3 m PWD was noted in the south-central part (Motijheel area) of the city (Figure 11a).
The shape, extent and magnitude of the depression cone changed significantly in th year 1990 (Figure 11b). The cone covered a larger area than in the past and the lowes depression was further extended to the southeastern part of the city, in the Sabujbagh area, where the lowest PL was observed at −17.57 m PWD. The PL elevation contour −1 m PWD covered most of the Sabujbagh and Khilgaon areas, and part of the Motijheel area

Spatial and Temporal Distribution of Piezometric Level for UDA
Although groundwater development started for public supply in 1941 in Dhaka [45,57,88], due to the unavailability of piezometric data, it was not possible to make piezometric maps before 1980. The existing extraction of groundwater formed a local cone of depression in 1980. The PL contour of 0 m follows the outline of the city in that period. The contours of −1 and −2 m PWD encompass the local depression due to groundwater pumping. Though PL was very close to the ground surface in the 1980s, a depression of −3 m PWD was noted in the south-central part (Motijheel area) of the city (Figure 11a).
The shape, extent and magnitude of the depression cone changed significantly in the year 1990 (Figure 11b). The cone covered a larger area than in the past and the lowest depression was further extended to the southeastern part of the city, in the Sabujbagh area, where the lowest PL was observed at −17.

Discussion
The analysis of long-term hydrographs and piezometric maps reveals that the pattern and magnitude of the PL decline from the 1980s onward are directly proportional to the population growth associated with the city's expansion ( Figure 2) and number of DTWs installed (Figure 6) in particular parts of the city in a certain period. For a better understanding of the trends in PL decline, synthetic graphical overviews of piezometric trends of all the UDA observation wells in Dhaka city were plotted ( Figure 13). In these plots, each observation well is plotted as an individual piezometer to observe the evolution of PL. Based on the synthetic graphical overview of the piezometric trends (Figures 13 and  14) and the long-term hydrographs, it is evident that PL is declining at an increasing rate due to abstraction for urban demand. In the 1980s, a local cone with the lowest decline of -3 m PWD was observed in the south-central part of the city, in the Motijheel area. The center of the cone shifted towards the southeastern part of the city in 1990, with the maximum depression at −17.57 m PWD in the Sabujbagh area. The declining trend laterally extended to the south-central part and the northwestern part of the city.
In 1990, most of the observation wells showed linear decreasing trends, but a very steep declining trend was observed in Dhanmondi (southwestern region) and Motijheel areas. After 1990, most of the piezometers with depths of 30 m ran dry in the south-central to central regions and also in the southwestern area (Motijheel, Sabujbagh and Mohammadpur) (Figure 8a,b,g). In 2000, the depression cone was marked by the lowest PL elevation at −49.5 m PWD in the Ramna and Motijheel areas. Due to population growth in the 1990s and 2000s, the new Dhaka north expanded to its present boundaries in Mohammadpur, Mirpur, Tejgaon, Banani, Baridhara and Uttara [62,64]. The massive abstraction by these DTWs caused a PL decline, centered on the area of pumping (Figure 11b), but laterally the area of influence was observed outside of the city before 2000.
After 2000, most of the observation wells with depths up to 50 m ran dry in Sabujbagh, Banani, Mirpur and Dhanmondi, except in the northeastern and southwestern areas (e.g., Sutrapur well). Due to the huge PL decline, in 2004 DWASA drilled wells to greater depths than they had previously [80] and installed more than 100 DTWs between 2000 and 2005 both in UDA and MDA [26]. In 2005, most of the observation wells, including twin monitoring wells, declined very sharply and linearly, except for the observation wells in the southwestern region (Lalbagh and Gandaria) ( Figure 13). In 2005, a dramatic decline was observed in the entire city due to pumping by the private and DWASA wells. The highest rate of drawdown in UDA (5.1 m/year) was observed in the Mirpur area, followed by a constantly lowering rate. At the same time, in MDA, three sparsely distributed depressions were observed throughout Dhaka city, with the lowest PL elevations of −36.82 m PWD in the southeastern part in Shympur and Sabujbagh areas, −34.62 m PWD in the northwestern region in the Mirpur and Pallabi areas, and −21.35 m PWD in the Dhanmondi area in the south-central to southwestern part of the city.
Since 2005, the PL has declined more sharply in MDA than in UDA in the central part of the city, as is clearly reflected in the synthetic graphical overviews of the piezometric trends of both aquifers (Figure 14), the twin hydrographs ( Figure 9) and the piezometric maps of UDA and MDA (Figures 11 and 12). In 2010, only a few observation wells were working in UDA and these mainly showed both linear and exponential declining trends. In 2010, the depression cone of UDA was marked by the lowest PL contour at −70 m PWD in the Khilgaon and Tejgaon areas (Figure 11d). In MDA, the lowest PL at −68 m PWD was observed in the central part, which indicated that a very sharp decline of PL occurred in MDA in 2010. A regional depression cone was formed in this aquifer in 2010 and the area of influence included the peripheral parts of the city (Figure 12b).  Since 2016, the cone with the lowest PL contour of −65 m PWD extended towards the northwestern part of the city. In 2017 to 2018, the depression contours were similar to each other, but the cone slightly extended from the southeast towards the south-central direction ( Figure 12 and Figure 16). Since 2016, the twin wells of the central part and peripheral part of the city showed a sharper drawdown in MDA than UDA due to intensive extraction from MDA (Figures 13 and 14). From 2017 to 2018, the highest rate of decline was observed in the southeast part of the city, in the Khilgaon area, both in UDA (4.0 m/year) and MDA (5.74 m/year) (Figure 9a). The relatively lower drawdown observed in the eastern region in UDA was possible due to lower urbanization and the construction of most of the DTWs in MDA.
The depression cones reflect the pattern of change of groundwater extractions and city expansion in different time periods (Figures 2 and 6). Massive and continuous extraction has modified the directions of groundwater flow by reversing hydraulic gradients towards the central and south-central to southeast parts of the city, where the lowest depression was noticed in both aquifers [89]. The 3D models of both aquifers show the radial flow towards the center of the depression cones in Dhaka city (Figures 15 and 16). These two models illustrate the lateral and vertical spread of depression cones and the peripheral piezometry in the UDA and MDA in Dhaka city.
The consequences of this over-pumping may also prompt environmental hazards like subsidence, negative impacts on water supply, drought and ecological imbalance [14], as well as permanently damaging the aquifers [90]. The vertical recharge of UDA and MDA is very slow, as they are semi-confined aquifers, and lateral flow from the surroundings provides most of the groundwater flow to the depression cone (Figures 15 and 16). For these reasons, the peripheral parts of the city, even the less urbanized parts, also experienced the lowering trend of PL in both aquifers.
The different water balance components of the groundwater reservoir in Dhaka city are indicated in Figure 17. Inflow into the semi-confined Dupi Tila aquifer system in Dhaka city is caused by lateral inflow (L) from the surrounding areas and vertical leakage from the overlying and underlying clay layers. The importance of inflow from the rivers (IR) is governed by the hydraulic resistance between the river and groundwater. Several studies [41,43,77] have postulated direct contact between the river and the aquifer. However, the drilling data we consulted all showed that the upper aquitard is continuous, and this should limit direct infiltration from river water into the semi-confined aquifer. An indirect connection through the aquitard is nevertheless suggested by the steepened gradient of the depression cone in the southwest, near the Buriganga River. Outflow from the Dupi Tila aquifer system in Dhaka city in today's situation occurs fully through pumping (Q). The water balance is then given by: Lateral inflow + Leakage + Inflow from Rivers − Pumping = Change in Storage (1) An estimation of the contribution of the storage decrease to the pumped discharges is given below. In 2020 around 900 million m 3 /year was pumped in the study area, which has an area of ca. 358 km 2 . This accounts for a discharge flux of 2515 mm/year in 2020. The intensive exploitation of the aquifer started around 40 years ago, and assuming for simplicity a linearly increasing trend in pumping rates, the total discharge flux over the 40 years is around 50,300 mm, or a water column of more than 50 m. During these 40 years, PLs have been lowered systematically. The most recent piezometric maps show that the average PL under the city is now around 60 m below PWD-deeper in the center of the depression cones, but higher near the rivers, in both the UDA and MDA aquifers. With an average thickness of resp. 98 m (UDA) and 86 m (MDA), a total of 184 m, and using an estimated value of 10 −5 m −1 for the specific elastic storage, the total storativity of the pumped aquifer system will be 0.00184. The drawdown of 60 m has thus released a total of ca. 110 mm of water from storage during the last 40 years. Storage has thus contributed 0.11/50 = 0.2% of the pumped discharges. The pumped water is thus mainly derived from lateral flow and vertical leakage (part of which is derived from river water). Quantification of the relative contribution of both components requires a groundwater model, which we intend to build in the future. lateral flow and vertical leakage (part of which is derived from river water). Quantification of the relative contribution of both components requires a groundwater model, which we intend to build in the future.  As only a small fraction (<1%) of the pumped groundwater comes from storage, the present piezometric distribution is in a near steady state flow regime and the declining piezometric trend over the last decades is mainly due to a systematic increase in pumping rates and not to depletion of the groundwater storage caused by a lack of inflow. The depth of the depression cones is thus correlated with the abstraction rates.
Rainwater harvesting and the use of treated surface water for different purposes could be good options to reduce groundwater exploitation. Artificial recharge could also be an option to maintain a lasting water supply and to prevent further declines in piezometric levels in the future. Groundwater flow modeling can be used to quantify the water balance components of the aquifer system of Dhaka city, and we recommended that this is done in the future. As only a small fraction (<1%) of the pumped groundwater comes from storage, the present piezometric distribution is in a near steady state flow regime and the declining piezometric trend over the last decades is mainly due to a systematic increase in pumping rates and not to depletion of the groundwater storage caused by a lack of inflow. The depth of the depression cones is thus correlated with the abstraction rates.
Rainwater harvesting and the use of treated surface water for different purposes could be good options to reduce groundwater exploitation. Artificial recharge could also be an option to maintain a lasting water supply and to prevent further declines in piezometric levels in the future. Groundwater flow modeling can be used to quantify the water balance components of the aquifer system of Dhaka city, and we recommended that this is done in the future.   As only a small fraction (<1%) of the pumped groundwater comes from storage, the present piezometric distribution is in a near steady state flow regime and the declining piezometric trend over the last decades is mainly due to a systematic increase in pumping rates and not to depletion of the groundwater storage caused by a lack of inflow. The depth of the depression cones is thus correlated with the abstraction rates.
Rainwater harvesting and the use of treated surface water for different purposes could be good options to reduce groundwater exploitation. Artificial recharge could also be an option to maintain a lasting water supply and to prevent further declines in piezometric levels in the future. Groundwater flow modeling can be used to quantify the water balance components of the aquifer system of Dhaka city, and we recommended that this is done in the future.

Conclusions
In Dhaka city, hydrogeological data assisted in the delineation of three subsequent semi-confined aquifers (based on lithology) from the surface to the explored depth (463 m): the Upper Dupi Tila Aquifer (UDA), Middle Dupi Tila Aquifer (MDA) and Lower Dupi Tila Aquifer (LDA). The Dupi Tila aquifer system in the study area is a clear example of an overexploited aquifer, shown by continuously declining piezometric levels (PLs). Unplanned urbanization and a vast population increase have altered the natural groundwater dynamics in this city. The documentation of hydraulic heads for both UDA and MDA provides important information for groundwater resource management. The longterm hydrographs and piezometric maps of UDA and MDA provide both temporal and spatial scenarios regarding PL in Dhaka city. The exploitation and supply scenario are closely correlated with population and city expansion for the period from 1963 to 2018. The supplied water consists of more than 78% groundwater, pumped from DWASA and private wells. A sharp decline in PL values due to large abstraction was observed after 1990 in UDA and after 2005 in MDA. Large amounts of abstraction led to severe lowering of PL in both aquifers in 2010. In 2018, the lowest PL values were −89 and −70 m PWD, observed in the south-central part of the city, in UDA and MDA, respectively. In the center of the depression cone in 2018, the average rate of decline was 4.0 m/year in UDA and 5.74 m/year in MDA; the average rate of decline in the entire city was 2.25 m/year in UDA and 2.8 m/year in MDA. If the PL decline continues at this rate in the future, this will deplete the aquifers. Huge and prolonged extraction have changed the flow directions of groundwater in both aquifers by reversing hydraulic gradients towards the central and southcentral to southeastern parts of the city, where the peaks of the depression cones were observed in both aquifers. The long-term development and management of groundwater resources warrants comprehensive and multi-objective approaches that must include the cost-effectiveness of overexploitation and sustainability. In order to addressing sustaina-

Conclusions
In Dhaka city, hydrogeological data assisted in the delineation of three subsequent semi-confined aquifers (based on lithology) from the surface to the explored depth (463 m): the Upper Dupi Tila Aquifer (UDA), Middle Dupi Tila Aquifer (MDA) and Lower Dupi Tila Aquifer (LDA). The Dupi Tila aquifer system in the study area is a clear example of an overexploited aquifer, shown by continuously declining piezometric levels (PLs). Unplanned urbanization and a vast population increase have altered the natural groundwater dynamics in this city. The documentation of hydraulic heads for both UDA and MDA provides important information for groundwater resource management. The long-term hydrographs and piezometric maps of UDA and MDA provide both temporal and spatial scenarios regarding PL in Dhaka city. The exploitation and supply scenario are closely correlated with population and city expansion for the period from 1963 to 2018. The supplied water consists of more than 78% groundwater, pumped from DWASA and private wells. A sharp decline in PL values due to large abstraction was observed after 1990 in UDA and after 2005 in MDA. Large amounts of abstraction led to severe lowering of PL in both aquifers in 2010. In 2018, the lowest PL values were −89 and −70 m PWD, observed in the south-central part of the city, in UDA and MDA, respectively. In the center of the depression cone in 2018, the average rate of decline was 4.0 m/year in UDA and 5.74 m/year in MDA; the average rate of decline in the entire city was 2.25 m/year in UDA and 2.8 m/year in MDA. If the PL decline continues at this rate in the future, this will deplete the aquifers. Huge and prolonged extraction have changed the flow directions of groundwater in both aquifers by reversing hydraulic gradients towards the central and south-central to southeastern parts of the city, where the peaks of the depression cones were observed in both aquifers. The long-term development and management of groundwater resources warrants comprehensive and multi-objective approaches that must include the cost-effectiveness of overexploitation and sustainability. In order to addressing sustainability, a long-term management plan that includes the management of unplanned urbanization, population migration or control of growth and decentralization of industrial and economic zones is urgently needed, otherwise the situation may go beyond our control.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.