Advancing from Conceptual to Numerical Modeling: Reducing Uncertainty in the Data-Scarce Dili Intergranular Aquifer System, Timor-Leste

1. Introduction

Groundwater is essential for the economic and social well-being of urban populations in developing countries [1,2]. The growing dependence on groundwater, coupled with diminishing resources, has prompted the exploration of effective scientific methods for groundwater investigation and management [3], which are needed to understand the changing and complex nature of groundwater and its impact on ecosystems and communities [4]. To learn more about local groundwater systems, to predict their future states and to manage them sustainably, it is important to understand their current status [3,5,6,7,8]. Such an endeavor requires, as a first step, a solid conceptual understanding of an aquifer system, that is, a detailed picture that describes the hydrogeological processes occurring within the system. It may include geological, geophysical, hydrological, hydrogeochemical, and other details according to the specific context of the aquifer and the study objectives [9]. Comprehensive knowledge of local lithology is crucial to understanding the movement of water and its quality within an aquifer. Three-dimensional (3D) lithological modelling is typically used as a baseline in groundwater studies to assess the hydrogeological heterogeneity of the subsurface and define hydrostratigraphic units (HSUs) [10,11]. A hydrostratigraphic unit can be defined as a part of a body of rock/soil that forms a distinct hydrologic unit with respect to the flow of groundwater [12]. It was later redefined by [13] as "a body of rock/soil distinguished and characterized by its porosity and permeability." Such information is key for a solid hydrogeological conceptualization [10], especially in regions characterized by significant hydrogeological heterogeneity [10,14]. Usually, a conceptual model uses cross sections, fence diagrams, maps, and tables to show how HSUs are spread out, how groundwater flows, what processes are controlling recharge, and what the hydrogeological parameters estimates are [10].

An effective assessment of groundwater resources can greatly benefit from reliable and consistent groundwater simulations [7,14]. Realistic hydrogeological conceptual models are essential for designing numerical flow models [10,11,4,15]. Numerical groundwater modeling can be applied to many groundwater investigations, such as to better represent the hydrogeological system and refine the conceptual model, to understand the flow system, to study the recharge dynamics and rates, to assess the groundwater resources' potential, to test different groundwater management approaches, and to investigate future impacts on the hydrological system from changes to stress factors such as recharge and abstraction [16]. Numerical groundwater modeling is a way to simplify the ultra-complex reality of the hydrogeologic system of an area.

Dili, the capital of Timor-Leste, is experiencing increasing demand for freshwater due to both population and economic growth. The city is located on alluvial plains of loose sedimentary deposits, overlaying the Dili Sedimentary Basin, which contains highly productive aquifers that store enormous volumes of water and have facilitated the settlement of a large number of people [17]. The Dili population totally relies on the coastal Dili intergranular aquifer system (DIAS) for water supply, which, as such, is the object of this study. Many have studied groundwater in Dili City [17,18,19,20,21,22,23,24,25,26,27,28]. The studies conducted have contributed to this field. However, DIAS still significantly lacks hydrogeological knowledge, including robust recharge estimates and dynamics and a solid conceptual model encompassing a 3D lithological model to describe layered sediments and hydrogeological units.

Rapid urban growth and development in Dili have raised concerns about impacts on groundwater resources, and therefore a deeper understanding of how urbanization and increased groundwater extraction affect local water supplies is needed. Due to population growth and more economic activity, the use of groundwater in Dili is expected to rise significantly in the next few decades [22]. There is currently no proper groundwater management in the Dili basin [24], which will undoubtedly lead to sustainability problems related to water availability and quality, including seawater intrusion. The aquifers in the Dili Basin have good yields with 30 production wells calculated to have produced 10.6 million cubic meters (mcm) from January to October 2019, with most extraction in the western zone [17]. This estimate does not include thousands of existing private wells [17]. Further, collecting accurate data on water abstraction is challenging [8]. On the other hand, dropping water tables in some production wells have been reported [19]. Therefore, increased long-term withdrawal, without considering the aquifer behavior in response to recharge and pumping, will put DIAS at serious risk of degradation and compromise the ability to meet water needs. Despite some hydrogeological studies, including a simple modeling attempt, many questions remain about DIAS’s conceptual model. This, in turn, means water managers and decision-makers have little information or guidance to work with in relation to sustainable pumping, drilling, and contamination risks. Section 2, Description of Study Area, provides a more detailed discussion of existing studies and remaining hydrological questions.

The only groundwater numerical model of DIAS to date, to the best of our knowledge, estimated the safe yield for Dili's groundwater to be around 0.28 m³/s [22]. The authors considered a homogeneous aquifer (i.e., no different stratigraphic units) and homogeneous recharge across the basin. The lack of detailed data on sedimentary layers, recharge dynamics, and clarity regarding DIAS’ conceptual hydrogeological model prevented the authors from developing a more complex model. Only 14 pumping wells were used, which further limited their ability to model the system realistically. Therefore, further research is needed to improve the conceptual and numerical modeling of DIAS. Understanding the timing and magnitude of changes in aquifer storage is critical to assessing groundwater availability in the DIAS and its sustainable use. As such, the aim of this study is to gain a greater understanding of the Dili basin's hydrogeology and predict future groundwater availability. This work is done by creating a numerical groundwater flow model that allows us to simplify the complexity of the local geological and hydrogeological contexts. Hydrogeologists, planners and regulators need a flexible, advanced tool to identify key areas and rates of groundwater recharge and discharge and predict groundwater dynamics—this is usually done through numerical groundwater flow models [29]. MODFLOW is the most common groundwater modeling code used globally [30], and as such was chosen for this modeling exercise. More specifically, the objectives of the study are (i) to understand which variables might impact the water balance more (e.g. recharge, evapotranspiration, hydraulic properties, abstraction rates and locations, and hydrostratigraphy); (ii) to predict changes to the water balance according to different scenarios of recharge and abstraction changes; and (iii) to propose future work, arising from the knowledge and uncertainties here gained, that allows for further inquiry into the hydrogeology of the system and more-informed management decisions. This paper does not include the strategies for sustainable groundwater management concepts, nor do we propose policies and regulations to protect groundwater resources; we will discuss these in a separate paper. No previous studies explain in detail the lithological conceptual model and numerical model of the groundwater flow system in the whole study area. Therefore, this study fills the major gaps by creating a thorough conceptual model of the DIAS and predicting the groundwater flow patterns according to present and future abstraction and recharge dynamics.

2. Description of study area

2.1. Geographical context

Dili has a population of 280 000, considering the administrative posts of Dom Aleixo, Nain Feto, and Cristo Rei [31]. The rapid growth has led to more groundwater needed to sustain daily life and economic activities. This research focused on the alluvial aquifers of the Dili Basin, also known as DIAS, which is bounded by the sea to the north and limited by mountainous areas to the south (Figure 1). The study area is mostly flat, covering around 26 km², with elevations between 1 and 63 meters above sea level. It is highly urbanized, mostly covered by residential, office, and commercial buildings. Topography inclines from south to north with a slope gradient of 0.5-6%, which can be > 6% around the foot of the mountains in the southern area [25]. The aquifers are recharged by rainfall and river infiltration [19,21,25] and most likely also mountain-block recharge from the southern hills. The region’s hot tropical climate is heavily influenced by the West by West Pacific Monsoon and its mountainous topography. The wet season takes place between December and May and the dry season between June and November [32]. The mean annual precipitation is 994 to 1223 mm/year [25]. The mean temperature is around 27°C in coastal areas and around 22°C in the highlands [19].

2.2. Characteristic of Dili intergranular aquifer system

Geologically, the DIAS corresponds to the Dili Sedimentary Basin, which is overlapped by the downstream (flat) part of five watersheds: Comoro, Maloa, Lahane, Taibesi, and Becora. The DIAS can be divided into three zones according to different river basins (Figure 1) and recharge dynamics (as demonstrated later). For those familiar with Dili, Zone 1 represents the areas of Comoro, Madohi, Bebontuk, and Manleuana; Zone 2 represents Fatuhada, Kampung Alor, Bairo Pite, Farol, Colmera, Kaikoli, Vila Verde, Maskarinhas, Balide, Lahane, Audian, and Bidau; and Zone 3 consists of Kuluhun, Taibesi, Becusse, Camea, Becora, and Bedois.

Previous hydrogeological studies conducted across the area are summarized in Table 1. The DIAS stands on top of the Permian Aileu Complex, formed mainly by phyllitic schist, metapelitic, metasandstone and metabasic rocks [33,34]. The Quaternary sediments, such as sand, gravel and clay, are deposited above the metamorphic rocks. DIAS is a multiple aquifer system with various proportion layers of sand and gravel divided by clay [17,19,20,21]. The deposits comprise typical alluvial sediments, including gravel, fine to coarse sand, silt, and clay, arranged in lenses that are graded laterally and vertically [17,19]. Lithological borehole logs around the central to eastern areas (zones 2 and 3 show different sedimentary strata with different hydraulic properties from those in the western area (zone 1). These represent the general hydrogeological context that prevails in the central to eastern area, where the presence of other fluvio-deltaic complexes forms important aquifers. The depositional pattern and aquifer thickness within the sediment basin are linked to the current rivers, which comprise both floodplain and terrace alluvial deposits. In contrast, the coastal and central regions consist of marine sediments at specific depths and fluvial sediments exclusively in the upper layer, as documented by various drilling logs [17]. Due to limited data, the thickness of the basin is unknown, but some drilling reports indicate it may be more than 150 meters thick in the central area and gradually thin toward the south. The water level is shallower along the northern coast and deeper around the southern area, making it more vulnerable to seawater intrusion and bacteriological and chemical contamination in the northern zone [23]. Since the city lacks an urban plan and groundwater policy, groundwater quality is threatened by uncontrolled domestic wastewater and solid waste, as indicated by the presence of bacteria and nitrates [23]. Some wells have also been found to contain elevated levels of calcium, magnesium, and heavy metals [25].

The Comoro watershed has by far the largest catchment area, enhancing the infiltration potential due to larger volumes flowing through in this western part of DIAS. The results of pumping tests confirm this area's high permeability. The groundwater pumping rate through production wells is >10 l/s, while the pumping around the central to eastern area is <10 l/s, as observed in the field. The relatively homogeneous coarse sediment materials in the Western area reflected a good aquifer potential from groundwater exploitation [17,20,21,23,25], but less of the confining bed may also allow for groundwater contamination and seawater intrusion risk. Groundwater abstraction reduces storage and baseflow, which in turn reduces streamflow and impacts low environmental flows during dry periods [35], which is a risk around the Comoro River. In addition, several boreholes also report that groundwater quality (electrical conductivity value could be >1500 μS/cm) becomes poor when drilling reaches depths > 55 m around the North to Central Zone.

Groundwater recharge results from direct infiltration of precipitation and river infiltration resulting from runoff from the high mountains in the southern area ,17,20,2125. Additionally, it is likely that mountain-block recharge occurs in the interface of the Southern mountains with the Dili sedimentary basin, especially in the Comoro catchment because of its larger area. The connection between groundwater and rivers is complex because of the multi-layered flow nature of DIAS, the dry-wet season patterns system of Dili. Typically, the Comoro River flows continuously during the wet season and remains dry during the dry season. It is possible, however, to observe a small wetland with abundant water at the mouth of the Comoro River during the dry season, just above sea level, which is most likely an expression of the groundwater table, representing fresh water sitting on top of the seawater. The data indicates that the water table close to the coast remains reasonably constant, which is because of the constant hydraulic head at the coastline. All other rivers typically run because of rain events and are dry otherwise.

Inadequate regulation of groundwater extraction leads to the uncontrolled operation of many unreported wells. Many private and domestic wells operate without meters, causing an underestimation of actual withdrawals. In 2021, Comoro had 4,373 wells as identified [36]. In the central and eastern zones there were fewer wells, probably because of geological conditions that are less favorable for drilling and/or water quality, although there is no comprehensive well register for those areas. There is no updated information on this matter.

3. Materials and Methods

The methodology for this study includes, firstly, a combination of several data sets to create a lithological model, which is followed by the hydrogeological conceptual model. Secondly, a steady-state groundwater flow model is set up and calibrated, which includes estimations of groundwater recharge, sensitivity analysis of hydraulic parameters and other water-balance variables, and definition of boundary conditions. Analyses of model results assisted in refining the conceptual model. Thirdly, future possible scenarios of abstraction and recharge are tested through a transient groundwater numerical model, which helps understand how DIAS will respond to changes in hydrological stresses over the next 30 years.

3.1. Conceptual model development

3.1.1. Data Source

Hydrogeologists use surface and subsurface data to reconstruct these systems to enhance their understanding of hydrogeological systems and improve quantitative modelling [37]. The borehole database created for this work includes information such as elevations, details of sediment layers, and water levels, which was collected from local drilling companies (e.g., H2O Pty. Ltd. 2007–2018; Mira Mar Unip. Ltd. 2014–2023; Geotechnic Pty. Ltd. 2019–2022); spatial geological data (including shapefiles) from the Instituto do Petroleo e Geologia (IPG) database (2014); and the Government of Timor-Leste (2020). Additional groundwater depth data was obtained from: (i) IGTL's database, which was then combined with field data collected for this study (2018–2023); (ii) the Government of Timor-Leste (2020); and (iii) four piezometer monitoring data from 2020–2021, obtained from the National Directorate for Water and Sanitation Regulations (ANAS/DNRAS, 2021). These data were combined with a 5-meter resolution LIDAR digital elevation model (DEM) obtained from the IGTL database and used to determine the DIAS water table elevation.

3.1.2. Lithostratigraphic model of DIAS

In this study, individual borehole coordinates, elevations, and hydrostratigraphic unit intervals derived from borehole lithological logs were added to the Leapfrog Geo software (Figure 2). This is a commonly used, user-friendly application for three-dimensional modeling of subsurface geology and hydrostratigraphy, surface, and subsurface data while interpolating borehole data to construct a three-dimensional spatial model [38]. The conceptual model for the DIAS is based on the surface elevation DEM. Lithological data from 157 borehole logs, including information gathered from field investigations, were used to conceptualize the sedimentary deposits into six strata and to subdivide the hydrostratigraphy into a two-unit hydrostratigraphy (Figure 7-b) with respect to the study area boundary. Delineation of these units divides the sedimentary layer into zones of different hydraulic properties and geometries. This division allows for a more comprehensive mapping and understanding of the groundwater flow system [39].

3.1.3. Hydrogeological Conceptual Model of DIAS

The conceptualization is a necessary simplification of the natural system, as it is not feasible to incorporate all of its complexities into a computer model [40]. Lithostratigraphic models help make more realistic hydrogeological conceptual models that summarize the hydrogeological knowledge of the modeled site, giving a framework for designing numerical models [10]. Usually, a conceptual model uses cross sections, fence diagrams, and tables to show HSUs’ architecture, maps of groundwater heads and flow directions, analysis of recharge/discharge dynamics and rates, and estimates of hydrogeological parameters [10].

The conceptual modeling of DIAS includes lithostratigraphic modeling, understanding surface-groundwater interactions, recharge dynamics, and the analysis of groundwater levels and flow directions. This procedure, in turn, allowed the selection of appropriate boundary conditions for the numerical model. Crucially, the numerical modeling calibration process allowed for refinement of the conceptual model. This was done by analyzing the performance of the model against measured water levels according to different conceptual understandings of recharge rates and processes. This approach is in close alignment with modeling best practices suggested by [41]: "There should be an ongoing process of refinement and feedback between conceptualization, model design, and model calibration such that revisions and refinements to the conceptual model can be made over time."

3.2. Steady-state numerical groundwater modeling using visual MODFLOW

Visual MODFLOW Flex was employed to simulate groundwater flow [42]. MODFLOW simulates groundwater dynamics within the aquifer by accounting for confined and unconfined conditions, using the finite difference method and a block-centered approach. It incorporates multiple recharge and discharge sources, such as recharge, groundwater extraction by wells, runoff from precipitation, flow to/from riverbeds and sea, and spring flow [43]. It can simulate multiple or single aquifer layers, accommodate various boundary conditions and anticipate future aquifer conditions using different pumping and climate scenarios [44]. Many studies have successfully applied combinations of GIS and MODFLOW to develop conceptual models and groundwater flow models in their study areas [10,16,30,28,5]. MODFLOW can be a flexible and powerful tool for simulating groundwater flow in a range of subsurface environments, from simple to complex [45].

3.2.1. Model setup and assumptions

The hydrostratigraphic and conceptual models were used as the basis for the steady-state groundwater flow model. The study area was discretized by dividing it into 125 rows and 250 columns. Thus, the total number of cells over the area is 31,250, with each cell measuring 45 m x 45 m in each layer. The dominant topographic slope of the area determined the main regional flow direction, which guided the orientation of the grid.

The following assumptions were made for the modeling:

• The river partially penetrates the aquifer system and has vertical banks;
• The river is not separated from the aquifer by any confining material;
• The influx into the system is primarily through recharge due to rainfall;
• Naturally, rivers interact with groundwater through river leakage. The phenomenon occurs in two forms: inflow (positive leakage), where water flows from the river to the aquifer, and outflow (negative leakage), where water flows from the aquifer to the river;
• The flow in the aquifer system is in steady state;
• Each lithological unit is homogeneous and isotropic;
• The groundwater system is assumed to be a porous medium, which has homogeneous and isotropic properties in the horizontal direction (Kx = Ky);
• The horizontal hydraulic conductivity value is assumed to be 1 order of magnitude greater than the vertical hydraulic conductivity value (Kx=Ky=10Kz);
• The southern hilly part and the basement, in the form of metamorphic rock, are assumed to be an impermeable layer and are therefore made inactive cells in the model;
• The model was calibrated for the period 2008–2023, using drilling and pumping data (2008–2022), recharge and evapotranspiration data (2021–2023), and river conductance and boundary condition data. Other hydrogeological parameters, including vertical hydraulic conductivity, specific yield, and specific storage, were assigned based on reference values from the literature [46,47,48], due to the absence of local field measurements.

Visual MODFLOW Flex imports data in ASCII or Surfer (.grd) formats [49]. In this study, surface elevation, layer elevation, and groundwater level data were formatted as Surfer (.grd) files and used as model inputs. [49]. The shape file layers or input parameters necessary for simulating groundwater flow in Visual MODFLOW Flex include (1) active and inactive cells, (2) elevations of the top and bottom of each layer, (3) horizontal and vertical hydraulic conductivity, (4) recharge rate, and (5) initial hydraulic heads.

Figure 3. Flowchart showing the development of the lithological model

Figure 3. Flowchart showing the development of the lithological model

2.2.2. Aquifer hydraulic parameters

Aquifer parameters, such as hydraulic conductivity (K), specific storage (S_s), and specific yield (S_y), govern groundwater flow [50]. Based on the lithological description and limited aquifer test data, the domain was divided into six units, each with a different K. Alluvial sediments typically exhibit significant heterogeneity, with K values varying across various orders of magnitude [51]. Since there are limited aquifer properties data for the area from pumping tests, we use the existent local data in combination with reference for K [47] (Table 2). As for the other hydraulic properties, namely S_s and S_y, we used MODFLOW default values in the absence of any accurate local data. The initial K values were modified within appropriate ranges during steady-state calibration to obtain best-fit models. Initially, the model used a single K value for both clay layers. However, after calibration, it became necessary to adjust these values separately to better reflect field conditions observed in some wells. Field conditions indicate that the upper clay layer consists of sandy clay, while the lower layer is made up of pure clay. Due to their differing properties, separate K values were assigned to each layer. The initial and calibrated parameters in the study area are presented in Table 2.

Figure 4. Spatial distribution of the lithostratigraphic units, each corresponding to calibrated horizontal hydraulic conductivities for the groundwater flow model.

Figure 4. Spatial distribution of the lithostratigraphic units, each corresponding to calibrated horizontal hydraulic conductivities for the groundwater flow model.

3.2.3. Groundwater Recharge

Groundwater in storage has a changing value that fluctuates in response to seasonal and long-term changes in recharge and discharge rates to and from the aquifer. Recharge is usually defined as the water from precipitation and surface-water bodies that infiltrates deep enough below the ground surface to reach the water table, where it joins the groundwater zone [52,53]. Recharge cannot be directly measured [54,55], and it is difficult to be estimated with a high degree of certainty [56]. In the absence of data that allow us for more comprehensive recharge calculations, the empirical recharge estimate methods (Equations 1-4) were also used to estimate how groundwater recharge is distributed. Since DIAS has different spatial hydrological conditions (Figure 1), it is relevant to input the recharge rate of each zone. The quantification of recharge in the study area was done using the simple water balance equation, where recharge (R) is equal to precipitation (P) minus real evapotranspiration (ETr) minus run off (R_o). ETr was calculated by applying Turc model [57].

Surface runoff was calculated by applying the empirical runoff equation from the Department of Agriculture of India (1990), which has been shown to be suitable for the island of Java (Institute Teknologi Bandung, 2001) and has been used in various Indonesian studies [58,59] for Jogjakarta and Sumbawa, Indonesia. DIAS is consistent with the Indonesian geographical context; thus, this method is considered adequate to be applied for our case study.

DIAS is recharged by rainfall and river infiltration [16,18,22], and most likely also mountain-block recharge from the southern hills. River infiltration and mountain block are processes particularly relevant for the Comoro River, which has a large, mountainous 207.29 km² basin. Typically, mountain blocks recharge constantly and slowly, bringing groundwater into the aquifers at the interface between the mountain and the edge of the sedimentary basin. Additionally, recharge at the mountain front is also likely to play a role in the overall recharge process, which takes place through the infiltration of water from streams located at the mountain front. It is likely that groundwater is directly subject to evapotranspiration when it is very shallow, which happens close to the coastline and riverbeds (especially the Comoro River). The DIAS area was divided into three main recharge zones: west, central, and east, based on the distribution of monitoring stations and the river catchments of each area (Figure 1). As such, the recharge calculations as per Equations 1–4 were performed for each of those areas, which only represent direct rainfall recharge.

As there is no data available to calculate these additional sources of recharge (i.e. rivel infiltration, and mountain block and front processes), we opted for calculating recharge for the upper Comoro catchment using the same methods as above and adding that recharge to zone one (as zone one is the end of the Comoro catchment). We assume the same recharge rate in mm/year over zone one, which in turn means a reduction in total volume recharge of 88% because of the small area compared to the total upstream area. This might be a gross simplification and limitation for recharge calculations, but in the absence of any data that can allow for more accurate estimations, it is a reasonable approach. Future studies should refine recharge calculations based on target fieldwork and new datasets.

For the DIAS direct recharge calculations, three years of rainfall data (2021–2023) from five stations (Table 6) were obtained from DNRAS, the National Directorate of Meteorology and Geophysics (DNMG) and IGTL. Temperature data were sourced from DNMG records (focusing on the flat area). However, due to limited climate data for the same period in the upstream Comoro basin, rainfall records from 2012 to 2019 were used [17]. We acknowledge these data sets are very limited, but the data availability in Timor-Leste is a major limiting factor for hydrogeological work, and our analysis is the best that can be done with available information in a timely manner to model the aquifer.

3.2.4. Boundary conditions and initial conditions

To describe how water enters or leaves the simulated aquifer system, it is important to consider the boundary and initial conditions, the elevation of the top and bottom layers, the recharge rate, K, the abstraction of groundwater from wells, rivers and bodies of seawater. These factors were used to run the flow model. The next step was to construct the model, primarily by converting the conceptual model into input files for the numerical model. Finally, calibration was carried out to demonstrate that the model could produce field-measured heads and flows, which were then used as calibration values.

Figure 5. (a) Map showing grid pattern: active cells, inactive cells, and boundary conditions, (b) Map of watershed boundary of each zone.

Figure 5. (a) Map showing grid pattern: active cells, inactive cells, and boundary conditions, (b) Map of watershed boundary of each zone.

The following boundary assumptions were made for the modeling:

No-flow boundary: The model base is assigned to a zero-flow boundary, as the bedrock underlying the principal aquifer is poorly permeable. The no-flow boundary was assigned in the hilly region as an inactive cell located in the southern part of the study area, which contains metamorphic rocks; therefore, the western, eastern, and southern boundaries are not included in the model;

• Constant head boundary (CHB): The constant head boundary in this model is used to represent boundary conditions that have a fixed head and do not change with time. In this model, the constant head boundary is used to represent the coastline, so the constant head boundary is used to simulate groundwater discharge to and from the sea;
• River: The river boundary conditions represent the Comoro, Maloa, Lahene, Taibesi/Kuluhun, and Becora/Bidau rivers, which generally flow from south to north. All rivers are dry during the dry season and flow only during the rainy season. In the dry season, the Comoro River retains water only in its upstream section, while the downstream section is dry, indicating it is a losing stream. This suggests that surface water from the Comoro River infiltrates into the aquifer, contributing to groundwater recharge. The general parameters of the river boundary conditions are given in Table 3.

In the river package of MODFLOW, river boundary conditions are used to represent natural boundaries like rivers in the Dili area. Groundwater level data from monitoring wells helped establish the initial head and served as calibration data for the steady-state model.

Regarding the governing equation for analyzing the river boundary condition, we follow the MODFLOW manual classic river package or the MODFLOW-6 gwf-riv package. The river package requires river stages, riverbed elevation, river conductance, and aquifer head. We considered river stages, bed elevation, and conductance in three segments for each of the four selected rivers. In the absence of field data, the initial conductance of the riverbed was set up as default from MODFLOW and was adjusted during the calibration period.

3.2.5. Modelling

MODFLOW is a modular three-dimensional finite difference groundwater modelling code [60]. In it, the developed conceptual model is converted into a mathematical model that simulates groundwater flow indirectly by means of a governing equation thought to represent the physical processes that occur in the system, together with equations that describe heads or flows along the boundaries of the model [61]. The equation representing the transient flow of a compressible fluid in a non-homogeneous anisotropic aquifer is derived by integrating Darcy's law with the continuity equation. Groundwater modelling was applied using the equation for three-dimensional groundwater flow [48,60].

3.2.6. Steady state calibrations

In the absence of pre-development data, or any other meaningful dataset, the steady-state groundwater flow model has been constructed by using the water level of September to November 2021 and calibrated to represent the natural flow system, where it represents a natural balance regime of inputs and outputs from the aquifer system. The calibration involved matching the observed heads in the aquifer with the hydraulic heads simulated by Visual MODFLOW. Three years of measured water level data from monitoring wells combined with initial water level data from domestic wells within the model area were used for calibration, with hydraulic conductivities and local recharge rates as calibration variables.

To date, just four monitoring wells have been installed on the north coast and in the western part of Dili in recent years [17,22], but they do not cover the central and eastern areas. We selected DM1, DM3 and Kuluhun-A and then combined it with water level data from 30 households for steady-state calibration, with the hydraulic values used for the model heads, with the observed heads representing the natural groundwater piezometric surface.

A series of model runs were carried out to match the observed and simulated water level contours by varying the K values and boundary conditions. The calibration parameters were modified until a satisfactory match between simulated and observed hydraulic heads, considering a root mean square error (RMS) < 10%, and baseflow values was achieved to better represent the actual conditions.

3.2.7. Sensitivity analysis

Sensitivity analysis is a critical component of modeling studies due to the uncertainties associated with estimating aquifer parameters, stresses, data availability and boundary conditions. The aim is to assess the impact of different model parameters and hydrological stresses on the aquifer system and to identify the most sensitive parameters that require focused investigation in subsequent research [50,62]. In the present study, the sensitivity of the parameters was assessed by varying the values of four variables: recharge, evapotranspiration, river conductance and hydraulic conductivity.

3.2.8. Predictive transient modelling

Using the steady-state calibrated groundwater flow model and applying the steady-state heads calculated by the model as initial heads for the transient model runs, we considered current abstraction and recharge from DIAS as an input to a transient model to test the current conditions to the best of our ability. Following this, we considered different scenarios of abstraction accounting for population/demand growth and of groundwater recharge accounting for possible changes to land use and climatic conditions (Figure 6). The intention was to test how the DIAS hydrogeological system will respond to plausible future scenarios of abstraction and recharge. The predictive transient model was performed for 30 years, representing the period 2024-2054, with changing conditions every 10 years to account for socio-environmental dynamics.

3.2.9. Estimation of current groundwater abstraction

Groundwater abstraction in the study is through dug wells, tubewells, and boreholes. There are three main types of abstractions, namely domestic, public supply and private for commercial enterprises (bottling water companies, hotels, etc.). Water use per capita varies significantly around the world, typically ranging from 110 to 160 liters per capita per day for middle-income households [63]. In determining the pumping rate in the MODFLOW Flex modelling in this case, the approach used refers to the daily water consumption habits of urban communities in Dili City, which we consider to be between 100 and 120 liters per person per day.

The average Dili City household consists of five to six family members [30]. Accordingly, each household's estimated daily water demand falls between 700 and 1,300 liters, or 0.7 and 1.3 m³/day. From this range, a mean value of 1.0 m³/day was chosen to serve as the modeling's representative pumping rate for domestic abstraction. This figure gives a realistic approach to the groundwater pumping load in the flow simulation and represents typical water consumption conditions in Dili City's residential areas.

Production wells operate continuously for 24 hours a day. There is no complete pumping dataset for public supply nor commercial wells. Therefore, we estimated pumping for these types of wells based on the limited available recorded abstraction data and on maximum rates recorded in pumping tests. The amount of groundwater used by industry and pumped from different public wells varies greatly. Small industries use an average of 144 m³/day, while large industries and production wells for public water supply use an average of 576 m³/day. Thus, we used this range for input as groundwater abstraction for commercial enterprises and public supply wells on MODFLOW through the well package.

3.2.10. Suggested future scenarios of abstraction and recharge

We expect significant population growth in Dili, along with economic development and associated lifestyle changes. This will most likely increase groundwater demand and pumping, although it is not realistically possible to predict by how much. Thus, we created extreme but plausible scenarios, in which it was assumed that domestic pumping would increase by 32% and commercial/public supply pumping by 45% every ten years. This projection assumes an increase in domestic, commercial, and public supply pumps but also a higher rate of coverage from the public supply system, which is still limited around Dili. We increased the pumping rate in each MODFLOW well instead of adding more wells, as we don't know where new wells would be.

Climate and land use change are likely to alter groundwater recharge dynamics and rates at DIAS, although it is not possible to predict, at least with current knowledge, how that will occur. Land use changes will likely decrease infiltration due to deforestation and impermeabilization, which might lead to an increase in runoff, diminishing diffused recharge but potentially increasing localized recharge (e.g., through riverbeds). The effects of climate change in Dili are unclear, as is its impact on recharge. Thus, we considered it useful to consider three possible recharge scenarios in relation to current conditions for each 10 years of future simulations: decrease by -7.5%, increase by +7.5% and no change. Although +7.5% is an arbitrary figure, it signifies a substantial change in recharge that can generate plausible future scenarios, which help hydrogeologists, decision makers, and managers understand the aquifer's future hydrogeological conditions.

In total, nine scenarios were created with combinations of pumping and recharge rates, summarized in Table 4. Scenarios 1, 2 and 3 were run from 2024 to 2034, considering current pumping conditions and different recharge percentages compared to estimations of the current conditions (0, +7.5, and -7.5%, respectively). Scenarios 4, 5, and 6 replicate scenarios 1, 2, and 3, respectively, while introducing a modeling period from 2034 to 2044 that accounts for a 32% increase in domestic pumping and a 45% increase in commercial/public pumping, along with the same changes to recharge conditions as previously noted (0%, +7.5%, and -7.5%). Scenarios 7, 8, and 9 repeat scenarios 4, 5, and 6 (respectively) and add a modeling period from 2044 to 2054, considering a further 32% increase in domestic and a 45% increase in commercial and public pumping in relation to 2034–2044 and the same changes to recharge conditions as previously (0, +7.5%, and -7.5%, respectively).

4. Results

4.1. Conceptual model of DIAS

The simplified lithological model has shown a great correlation with the borehole lithological logs (Figure 7-a). The thickness of Quaternary sediments typically increases from south to north. Drilling data confirmed that the sedimentary basin is thicker (> 140 m) around the central to western area, while the eastern area is < 100 m. The sediments of the western area are dominated by poorly- to well-sorted gravelly sand that vary in both size and composition, while the central to eastern area of Dili is formed by poorly sorted gravelly sand. The clay layer (thickness <10 m) is underlain in some areas by a thick granular zone, but the clay distribution is not continuous.

Table 5 summarizes the characteristics of the six identified lithological layers. The first layer is composed of gravelly sand with various proportions of these materials and has a thickness ranging from 1 to 68 m. The second layer is composed of clay, with a thick range from 0 to 30 m. The lateral distribution of the clay layer is uneven; it is absent in places, and it is thicker around zone two (Caicoli). The third layer consists of sand with gravel and boulders with occasional clay; the estimated thickness of this layer is 2 m to 108 m, being thicker around the Comoro area (zone one). The fourth layer is composed of unevenly distributed muddy sand mixed with clay, with a thickness ranging from 0 m to 20 m and being thicker around Zone 2 (Caicoli, Bairopite, and Hudilaran). Layer five is composed of sand with a variable thickness ranging from 24 m to 64 m, which is thicker around the Comoro area. Beneath these sediments, there is the metamorphic bedrock (layer six).

Although confined conditions may exist in specific locations, particularly in zone 2, due to the lack of monitoring data, it is not possible to distinguish between wells under confined or unconfined conditions. Thus, at the basin scale, the aquifer system can be considered to range from unconfined to semi-confined. Therefore, we modeled DIAS as an unconfined system. While this is a limitation that may be an overgeneralization, it is the best approach available with the current data. Future modeling efforts could address this issue if new data becomes available.

The observed vertical and lateral lithological variability indicates a heterogeneous sedimentary system that influences aquifer hydraulic properties. The presence of discontinuous clay layers may locally affect groundwater flow and storage conditions.

Based on the subsurface lithological model (Figure 7-a), groundwater modeling is performed considering two HSU. HSU-One represents bedrock, consisting of fresh to weathered schist. It can locally store and transmit water through fractures and weathering, but it is primarily composed of fresh schist that acts as an aquifuge or is impermeable—and it is modeled as such. HSU-Two represents Quaternary deposits mainly composed of three different confined and semi-confined aquifers separated by clay layers. HSU-Two contains five layers that control groundwater flow, as described in Table 5.

4.1.1. Groundwater level and flow direction

There is no comprehensive groundwater monitoring program in Dili. There are only four monitoring wells, of which only two are currently monitored monthly. For inaccessible wells, we assumed the groundwater level recorded after drilling was the initial water level and used it to construct the flow direction map. Otherwise, we took water level measurements, at different times of the year according to logistics. Although it is expected that static water levels may vary with the seasons, static water level measurements were considered representative of groundwater levels in the aquifer and acted as a surrogate for ambient groundwater conditions. Figure 7-b shows groundwater elevations ranging from 1 to 28 meter above sea level (masl). The general flow direction is from south to north, while it is southeast to northwest in the eastern region. Moreover, the wells located close to the coastline have a depth to the water table of ~0.5-2 m, which may allow for evaporation to occur directly from the groundwater.

4.1.2. Recharge

Rainfall amounts varied from 922 to 1,229 mm/year. Real Evapotranspiration values ranged from 568 to 622 mm/year. Recharge rates are 115 mm/year in the eastern region (zone 3), 78 mm/year in the central region (zone 2), and 65 mm/year in the western region (zone 1). For the upstream of Comoro, the recharge calculated was 247 mm/year. As explained in the methodology, for considering the mountain block recharge and river infiltration in addition to the direct rainfall diffused recharge, the recharge from the upstream of Comoro was added to zone one at the same rate (resulting in a total recharge volume reduction of 88%, due to Zone 1's smaller area compared to the total water volume available from the upstream area for recharge). These changes resulted in a total recharge rate of 312 mm/year in the western region (Zone 1). This is due to the region's size and its role as the main catchment area that carries water from the Comoro River system.

4.2. Steady state flow model

4.2.1. Calibration

Figure 8 shows a comparison of observed water levels and simulated hydraulic heads in the DIAS using a 1:1 correlation line. Table 7 shows that root mean square (RMS) and residual mean for almost all the sites during the calibration validation are low and within acceptable limits. Potential discrepancies between simulated and observed levels can be explained because of water withdrawal effects (evident in a few areas of intensive pumping), poor data availability, and errors in measurements (head) and estimations (baseflow, recharge). Initially, it was assumed that all cells in the original active zones of the model were uniformly wet. However, some of the model's original active cells, such as those in the eastern region, experienced drying out during the simulation and calibration processes. Dry cells represent areas where the water table is below the cell's bottom elevation (see Figure 9 and Figure 10). This may occur because of the boundary conditions chosen, because there is no flow, and/or because of the hydrodynamics in this thin and elevated part of the basin.

When using the MODFLOW default value for river conductance, a very high RMS was obtained. However, when the river conductance value was manually adjusted, the analysis revealed that a higher value produced a lower RMS. In addition to the RMS results, the river conductance value affects the type of river included in the model. The rivers in the area are of the losing stream type. However, the model results that use the default conductance value do not accurately represent the losing stream type. After adjusting the river conductance value to match the losing stream type observed in real conditions, the model results better represented these hydraulic dynamics, indicating that the final conductance values were more realistic.

We tried calibrating the model using two recharge scenarios. Firstly, we considered only direct diffuse rainfall recharge to the aquifer. This procedure resulted in an error message stating that the recharge was too low and the model was unable to run. Therefore, we considered the additional recharge coming from the upper catchment of the Comoro River (Table 6), and the calibration results were satisfactory with RMS <10% (Table 7). The result is an indication that recharge originating from the upper catchment is significant for DIAS water balance. This is an example of how the numerical model can be useful for a better understanding of the conceptual model.

4.2.2. Modeled groundwater balance

Recharge to the aquifer system is from precipitation, at 103 mcm, and river leakage, at 169 mcm. The total inflow to the aquifer is 272 mcm and the total outflow is 273 mcm (see Table 8). This results in a quasi-neutral water balance of ∆S = -1 mcm, which is consistent with what is expected in steady state – this ∆S represents 0.4% of the total inputs/outputs, which can be considered negligible seeing the data limitations and assumptions in the model. Table 8 shows significant outflow to the sea and no inflow, as well as significant infiltration from the rivers and outflow from the aquifer to the rivers depending on K and available water resources. There is also significant loss of water through evapotranspiration directly from the water table, which is not surprising given the very shallow water tables near the coast. The model shows that rivers mainly supply groundwater, which is consistent with the above remarks regarding losing conditions and field observations of river dynamics.

4.2.3. Groundwater level variation

The numerical model is an unique tool for understanding the spatial distribution of groundwater heads over time and according to different aquifer stresses (e.g., pumping, recharge) and boundaries. According to the steady-state model (Figure 9), the groundwater regime has the same general flow directions, from southwest to north, and is consistent with measured levels (Figure 7-b), with groundwater levels ranging from 0 to 66 m (Figure 9).

4.3. Transient model and scenario testing

The results in Table 9 and Figure 10a demonstrate the significant impact of pumping on the water balance of the system. Introducing pumping at the estimated current rates and maintaining recharge unaltered resulted in ΔS = -27 mcm by 2034. Increasing pumping created stronger groundwater deficits of –30 mcm by 2044 and –34 mcm by 2054 (Figure 11), indicating greater stress on the DIAS, especially in zones two and three where the water supply is deteriorating (see Figure 10). Infiltration from rivers is the main source of water in the system. Interestingly, outflow to the sea is the main water loss by 2034 and 2044, while abstraction becomes the main outflow by 2054, which is an indication of potential seawater intrusion as abstraction increases. As shown in Figure 10, several areas along the northern coast, including Madohi, K. Alor, Fatuhada, Farol, and Formosa, are at risk of saltwater intrusion due to the declining water table.

Despite the increase in pumping, the results from Table 10 and Figure 10 b show that the increase in recharge has a strong effect over time, contributing to a less negative ΔS: -25 mcm by 2034, -24 mcm by 2044 and -21 mcm by 2054 (Figure 11). The results indicate that recharge is a key aspect influencing the water balance and maintaining the sustainability of the system. Thus, significant future efforts should be invested to obtain more accurate estimations, and consequently reduce uncertainty around aquifer dynamics, and in measures to increase recharge. Infiltration from the rivers is the main source of water to the system. In line with S1, S4 and S7, these simulations show that outflow to the sea is the main water loss by 2034 and 2044, while abstraction becomes the main outflow by 2054, which is an indication of potential seawater intrusion as abstraction increases even if recharge also increases.

The results from Table 11 and Figure 10 c show that the decrease in recharge has a strong effect over time, contributing to an increasingly negative ΔS: -28 mcm by 2034, -37 mcm by 2044 and -47 mcm by 2054 (Figure 11). The results show that if recharge decreases overtime, the sustainability of the system is threatened. This is a further indication that significant future efforts should be invested to obtain more accurate estimations, and consequently reduce uncertainty around aquifer dynamics, and in measures to increase recharge. Infiltration from the rivers is the main source of water to the system. In line with S1, S3, S4, S7 and S9, these simulations show that outflow to the sea is the main water loss by 2034 and 2044, while abstraction becomes the main outflow by 2054, which is an indication of potential seawater intrusion as abstraction increases and recharge decreases. S8 represents the worst-case scenario for the water balance according to the assumptions considered, combining the highest abstraction with lowest recharge rates.

Figure 10 shows the water balance for all scenarios, according to the different components. It illustrates that river leakage is a strong component of the water balance, as the main source of water input into the system. This is compatible with field observations that large volumes of water flow through the rivers during the wet season and infiltrate quickly. Thus, it is crucial that future studies focus on estimating the potential for river infiltration more accurately and on surface-groundwater interactions more generally. Reducing river infiltration uncertainty will reduce water balance and conceptual uncertainty greatly and will allow for (better) aquifer management. As for outflows, flow to the sea and abstractions are the main components, depending on which scenario we consider. Thus, it is critical to have better estimates of groundwater pumping in Dili to constrain the water balance properly. A future numerical groundwater model considering density differences between fresh groundwater and seawater will also be important to better estimate discharge to the sea, and the risk of seawater intrusion – although that is of limited relevance until meaningful monitoring datasets exist.

Figure 9 shows the distribution of groundwater levels for the steady-state model, while Figure 12 shows equivalent maps for the nine different scenarios. Table 12 summarizes the variation in groundwater levels for the different scenarios according to the aquifer zones. The comparisons below are in relation to transient groundwater levels.

In scenario 1 (current pumping and recharge) he groundwater flow pattern remains similar to steady state (Figure 12a). Groundwater head ranges from 0 to 66 masl around zone 1, from 0 to 27 masl around zone 2, and from 0 to 33 masl around zone 3.

Scenario 2 (current pumping, +7.5% recharge) shows increased groundwater heads (Figure 12b) of 1-2 masl around zone 1, of 2-3 masl around zone 2, and of 2-3 m around zone 3.

Scenario 3 (current pumping, -7.5% recharge) shows a decline in groundwater levels (Figure 12c). Groundwater head drops by approximately 1–3 masl (the red zone narrows) around zone 1. There is a more noticeable drop in head, along with dense contours and steep gradients, around zone 2, where the head drops by approximately 2–3 masl. The green-yellow zone shifts northward around zone 3 (Table 12). The noticeable differences include the 12 m contour line around zone 1, which shifted ~100 m to the south due to water extraction from boreholes INAP 1 and 2, Acad Policia, G4, Top Fresh, and Comoro A; the 34 m contour line around the Comoro River also moved south by ~40 m; the 8-meter contour line moved southward around zone 2; in zone 3, the 22 m contour line shifted eastward, closely followed by the 20 m contour line near borehole Becora and Kuluhun; and in the same zone, contour lines between 30 and 32 m exhibited a similar pattern.

In scenario 4 (Figure 12d) (increased pumping, current recharge), groundwater heads show little differences.

Scenario 5 (Figure 12e) (increased pumping, +7.5% recharge) shows groundwater heads increase 1-2 masl around zone 1, the head increases 2-3 masl with sparse contour around zone 2 and the head increases 3-4 masl around zone 3.

Scenario 6 (Figure 12f) (increased pumping, -7.5% recharge) shows accelerated decline in groundwater levels by 1-3 masl around zone 1, by 2-4 masl around zone 2 and by 2-3 masl around zone 3 (Table 12). A cone of depression is observed around zone 1 due to abstractions, and the 12 m contour line shifted southward ~ 200 m. In zone 2, the 10 m contour line shifted southward ~ 50 m and the 16 m contour line around zone 3 is also shifted southward ~ 50 m.

For scenario 7 (Figure 12g) (further increased pumping, current recharge), the system remains stable but groundwater levels are trending downward, where cones of depression appear around zones 1 and 3.

Scenario 8 (Figure 12h) (further increased pumping, +7.5% recharge) shows head increases by 1–2 masl around zone 1 and by 2–3 masl around zones 2 and 3. In this scenario is noticeable taht the 4 m contour line shifted northward around Zone 2, as did the 16 m contour line around Zone 3.

Scenario 9 (Figure 12i) (further increased pumping, -7.5% recharge) shows a drop in the groundwater head of 1–3 masl around zone 1, 2–4 masl around zone 2 and 2–3 masl around zone 3. This resulted in the 4 m contour line being shifted southward in zone 2, while the 16 m contour line was also shifted southward, forming a cone of depression around zone 3. There is also the formation of a cone of depression between contour lines 8 and 12 m in zone 1, which includes the INAP 1 and 2 wells, the Acad Policia well, the G4 well, the Top Fresh well, and the Comoro A and B wells.

Comparisons between current and changed recharge conditions over a 10-year period (2024–2034) show that reducing recharge by 7.5% results in a significant decrease in groundwater levels, particularly in zones 1 and 2. Conversely, a 7.5% increase in recharge shows the opposite effect. The results consistently identify the western Comoro area (zone 1) as the primary recharge zone, while zones 2 and 3 showed less recharge. The results of the 20-year models (scenarios 4, 5 and 6) show the same trends: a decrease in recharge of 7.5% leads to a significant drop in groundwater levels in zones 1 and 2, while an increase in recharge of 7.5% stabilizes the system. Zone 1 continues to act as the dominant recharge zone, while zone 3 is the most sensitive to recharge fluctuations.

Reduced recharge accelerates depletion, most prominently in zone 2, over a 30-year period (2024–2054), making conservation strategies even more critical. Conversely, increased recharge increases aquifer reserves and maintains the hydraulic gradient. Across all tested timeframes, zone 1 consistently functions as the main recharge area, while zone 2 becomes the primary discharge zone and experiences the most pronounced fluctuations. Zone 3 also benefits significantly from increased recharge, as evidenced by its clear up-and-down responses.

5. Discussion

5.1. Comparison with previous modelling of DIAS

Despite a previous attempt to model DIAS, our work is highly relevant and novel because our modelling objectives, methodologies, datasets, complexity, boundary conditions, lithological heterogeneity, various recharge mechanisms, future scenarios, and results differ significantly [24], as summarized in Table 13. A key distinction between the studies is the objectives, which in our case were to refine the conceptual model using the numerical model and understand the hydrological variables that are more influential for the water balance, while theirs were to estimate safe/sustainable yield. All model choices and datasets have a strong influence on the results. An important example of a fundamental difference is the sea boundary, which we modeled as a fixed head, and modeled it as free drainage [24]. It is reasonable to assume that the average sea level fixes the head at sea, a common practice in modeling. An important consequence of this distinction is the water inflowing from the sea to the aquifer in our model, which in a simple way (because we did not consider density differences) represents potential seawater intrusion, an important threat for the sustainability of this coastal aquifer. This scenario was not identified because of the different boundary condition choices [24].

We recognize we benefited from having more robust and complete datasets, and thus we were able to create a more complex and realistic numerical model. The model takes into account lithological heterogeneity, various recharge mechanisms, boundary conditions, and future scenarios involving different recharge rates and pumping activities. Nonetheless, we acknowledge the importance of the first modelling attempt, according to the data available to them [24], which allowed us to improve upon their work (Table 13).

In what relates to recharge, our study is the first to distinguish three different recharge zones for the DIAS and to mention mountain-block recharge as a relevant mechanism (Table 6). We calculated recharge using empirical recharge estimation methods outlined in equations 1-4. The current study indicates a recharge rate of 312, 78, and 115 mm/year for zones 1, 2, and 3.

A previous study reported a recharge rate of 107 mcm/year based on the hydrological water balance method [16]. Another study reported recharge rates of 118 mm/year, despite not specifying the method [24]. Our study demonstrates that recharge and infiltration from rivers are essential elements of the water balance and should be integral considerations in future hydrogeological research and management strategies for DIAS, a point overlooked in prior studies.

5.2. Implications of the model for the sustainable use of groundwater

The numerical model here presented is, to the best of our knowledge, the most advanced tool for sustainable groundwater management in Dili, and anywhere in Timor-Leste. The main implications for sustainable groundwater management arising from this model are discussed as follows.

The western zone of the aquifer is more productive, and has higher groundwater potential, than the central and eastern zones due to its recharge dynamics. Consequently, if water extraction increases across DIAS, the water levels are expected to fall in the central and eastern zones over the next 30 years (Figure 9), which is depicted by the model's ability to forecast spatio-temporal changes in groundwater levels and flow directions. The model also shows that the north coast could be contaminated by seawater if extraction increases (Figure 9). As shown before, there is very little data hydrogeological data for the central and western zones, including a record of domestic wells, and as such acquiring more data for those areas is essential for aquifer protection.

The model evaluates potential impacts of climate and land-use changes, and water use, through different groundwater recharge and abstraction scenarios. When recharge drops, while extractions stay the same or increase, DIAS extraction rates exceed recharge capacity. Different water usage scenarios involving domestic and non-domestic users show the negative impacts of over-abstraction, providing useful evidence for setting groundwater extraction limits and establishing zoning policies. Three scenarios (S2, S5, and S8) specifically showed that the north and central regions can be drastically affected, with maximum drawdowns of 6 to over 10 m.

The model shows recharge, directly from rainfall, through river infiltration and/or from mount-block processes, is a key component of the water balance. Thus, significant effort needs to be put into characterizing recharge processes and calculating recharge rates better. Furthermore, creating conditions that allow for greater recharge will significantly impact the aquifer’s long-term sustainability. Such an effort could involve reducing impermeable areas in Dili, increasing infiltration throughout the river catchments (e.g., through reforestation and riverbank and soil protection), reducing river flow velocity and increasing infiltration in riverbeds. The measures to increase recharge would have the multipurpose of helping to reduce and manage flood risk, which is a chief water-related problem in Dili, and a threat to property and human lives. As majority of DIAS’s water supply relay on groundwater, a comprehensive groundwater management strategy is therefore required, considering environmental, technical, social and political factors. So far, the policymakers may have drafted policy, but it has not been translated into reality. For instance, the amount of groundwater abstraction in each zone remains undefined. Moreover, Dili is characterized by poor urban plans and a lack of wastewater management that can cause the shallow groundwater to be contaminated [21,22]. This study can assist policymakers in estimating sustainable extraction levels to prevent over-extraction, especially in the central and northern parts of the DIAS. Policymakers are also suggested to draft and enforce drilling policies to control and organize groundwater drilling activities, and to collect drilling data and create databases for increasing the hydrogeological knowledge and to keep understanding the evolving nature of the state of the aquifer system. Acquiring climatic data and crating databases shared across institutions is also fundamental to understand the aquifer system dynamics, which are presently nonexistent. It is necessary to improve groundwater levels and quality (e.g. salinity) monitoring to effectively control groundwater abstractions and their impacts on the aquifer. The groundwater monitoring system is important for (i) assessing the state of the aquifer in terms of how it will respond to pumping and recharge; (ii) evaluating the risk of seawater intrusion, contamination, and water availability for domestic purposes; (iii) evaluating the recharge and discharge dynamics from the river to the aquifer; and, (iv) monitor areas with larger aquifer withdrawals, as indicated by cones of depression. We propose the groundwater monitoring system shown in Figure 13, which is based on scenario modeling for the worst-case scenario (scenario 9 – low recharge, high abstraction). Ultimately, the management strategy should promote groundwater conservation and efficient use, and prioritize storage enhancement and distribution capacity to sustainably meet the population and environmental needs.

5.3. Limitations and the uncertainties in the model

The modelling attempt here presented is affected by severe limitations, the chief being data availability: (i) there is no groundwater monitoring network in Dili (only two wells are being measured monthly); (ii) there are very few measurements of pumping data; (iii) there is no register of domestic wells; (iv) there is no public register of drilling reports or well logs; (v) there are no river flow measurements; (vi) there is no continuous and comprehensive climatic monitoring; (vii) there is no field-based data that allows for the estimation of recharge, and spatio-temporal variability, with typical methods, such as the chloride mass balance or water-table fluctuation, and for the estimation of mountain-block recharge; (viii) there is no data on seawater intrusion; and, (ix) there is very limited information to distinguish between unconfined and confined conditions.

These data-availability limitations impact the ability to model the system accurately. Thus, the model here present can be considered an initial model that explores the conceptual uncertainties of the hydrogeological system but also that provides a solid base for future, more complex modelling efforts. The main limitations of this model are: (i) recharge was calculated with a single method based on limited, short-term local climatic data, which allowed very little insight into spatio-temporal variability and mountain-block recharge; (ii) many assumptions have to be made about abstractions from the system, as very incomplete data exist on this; (iii) despite the thousands of wells existing in Dili, we only had access to 157 to model the lithostratigraphy of the basin, which limited our ability to do so comprehensively at the basin scale; (iv) very little data exist on the aquifer hydraulic properties and their spatial variability, and as such we used the few pumping tests available and literature values to characterize the aquifer; (v) as this is an initial modelling effort, and with little data, we did not consider density differences and seawater intrusion, although this model sets the basis for such types of modelling considerations in the future; and (vi) as there are no measurements of river flow, the model is very limited in the ability to model the interactions between surface and groundwater, especially river infiltration and aquifer discharge.

Despite this long list of important limitations, the model represents a significant increase in the hydrogeological understanding of the Dili basin, including conceptual and how it may respond to climatic and abstraction stresses in decades to come. We consider it a fundamental tool for building new knowledge upon, and for science-informed groundwater management.5.4. Future work

Based on the study's findings, further research is needed in the areas of monitoring, conceptual understanding, modelling and management. As for monitoring, we suggest future work to include: (i) a comprehensive monitoring network of groundwater levels and salinities, potentially including data loggers in strategic wells, including assessing the potential for see water intrusion, and water-table fluctuations for groundwater recharge estimation; (ii) comprehensive monitoring of climate variables, both in the Dili area and upstream in the Comoro catchment, which can include collection of rainfall samples for characterizing rainwater quality and calculating recharge through the chloride mass balance method; (iii) a monitoring network of surface-groundwater interactions and how it changes throughout the year according to climate variability; and, (iv) monitoring of abstraction in a selection of public, commercial and domestic wells.

The monitoring efforts mentioned above can assist in a better conceptual understanding of DIAS. The main conceptual questions are: (i) what is the spatio-temporal variability of recharge?; (ii) What are the dynamics and rates of mountain-block recharge?; (iii) What are the dynamics of surface-groundwater interactions, and how do they influence recharge and discharge?; (iv) What areas are more at risk of seawater intrusion?; (v) In which zones do we find confined conditions, and how do hydraulic properties vary spatially with lithology, including in the riverbeds?; and, (vi) How do groundwater levels change throughout the year according to climate and pumping variability?

Comprehensive monitoring and enhanced conceptual understanding will allow for more robust modeling. Aspects to be developed in future modelling efforts are (i) density-driven seawater intrusion modelling; (ii) surface water flow modelling; (iii) fully coupled surface-groundwater modelling; (iv) automatic modelling calibration; (v) recharge modelling; (vi) a groundwater model considering unconfined and confined conditions; (vii) a model testing more complex scenarios of future recharge and demand and impacts on aquifer depletion, sweater intrusion and interactions with surface water; (viii) a model testing the effectiveness of different groundwater management approaches; and, (ix) contamination modelling.

6. Conclusions

Hydrogeological conceptual and groundwater flow modeling are essential tools for assessing groundwater availability and quality, understanding the processes controlling groundwater behavior, and ensuring the protection of groundwater resources. This approach was here applied to the coastal urban area of DIAS, which lacks groundwater management, has very little data and monitoring, and faces increasing unregulated groundwater drilling and abstraction. The current study, through lithological and numerical modeling, increases the conceptual understanding of DIAS’ hydrological system, highlights the knowledge gaps and data needs, illustrates the main processes controlling groundwater behavior and the water balance, and predicts the groundwater response to future scenarios of increased abstraction and variable recharge. Importantly, it provides recommendations for future work, including monitoring, conceptual understanding, modeling, and management approaches. The present study is an essential tool to assist policymakers in the sustainable development of water resource management and mitigating environmental risks in the city of Dili. It can serve as guidance in how to study similar aquifer systems in the world where little data and conceptual understanding exist.

The lithological model reveals that unconsolidated sediments such as gravelly sand, clay, and muddy sand, deposited above metamorphic bedrock, form the DIAS. The aquifer system consists of three unconfined to semi-confined aquifers, intercalated by two semi-confining aquitard layers. The bedrock is considered an aquifuge with low hydraulic conductivity. DIAS can be divided into three zones according to river catchments (five rivers flow into Dili’s sedimentary basin). Zone one has the thickest part of the sedimentary basin. It is characterized by poorly to well-sorted sediments, with thinner clay layers and thicker aquifers. Zones two and three are characterized by poorly-sorted sediments, with thicker clay layers and thinner aquifers. Thus, zone one has the highest aquifer potential in Dili. This portion is lithologically related to the alluvial plains of the Comoro River, the biggest in the region. The rivers flow from south to north, and so does groundwater.

Rainfall-distributed recharge rates were, based on climatic data, estimated to be variable across the aquifer: 65 mm/year for zone one, 78 mm/year for zone two, and 115 mm/year for zone three. Infiltration from the rivers and mountain-block recharge are additional recharge processes in the system, although there is virtually no data to calculate these. Recharge for the upper catchment of the Comoro River was calculated as 247 mm/year, and we added a portion of this recharge as an input to the numerical model to try to account for possible mountain-block recharge.

The steady-state model, through its calibration, showed good agreement between observed and modeled groundwater levels from 33 wells. It shows that rainfall diffuse recharge is not enough as input to the system, and additional recharge was needed for the model to run and produce better calibration results (<10% RMS). Transient modeling, considering current estimations of pumping and recharge, shows a negative water balance, with decreasing water levels and an increase in the risk of seawater intrusion. Future potential scenarios of increased pumping and variable recharge show that the system is at risk of a negative water balance, which is especially sensitive to recharge.

The results allow us to suggest future work for monitoring, conceptual understanding and modelling. We recommend a monitoring network of groundwater levels and salinities, abstraction, river discharge, surface-groundwater interactions and climatic variables. The results will allow for better conceptual understanding, especially in terms of the recharge calculations, areas of risk of seawater intrusion and aquifer depletion, and surface-groundwater interactions. Future modeling efforts can include seawater intrusion, contamination, more complex surface-groundwater interactions, a risk of subsidence, and automatic calibration. Without systematic monitoring and regulatory frameworks, groundwater management in Dili will remain highly uncertain and unsustainable.

Finally, this work allows us to suggest the following for sustainable water management: (i) thorough monitoring groundwater, climatic and surface water networks; (ii) regulation policies for pumping and drilling; (iii) measures to increase recharge and minimize flooding (e.g., reforestation of the river catchments, development of infiltration ponds throughout the catchments, prevention of soil impermeabilization and erosion, and structures throughout the rivers to gradually promote infiltration); and (iv) investing in thorough groundwater studies. This is crucial to address knowledge gaps and improve the management framework. The study promotes evidence-based groundwater governance in urban areas where data is scarce, such as Dili City. It also supports the achievement of several Sustainable Development Goals (SDGs), including SDG 6 (Clean Water and Sanitation), SDG 9 (Industry, Innovation and Infrastructure), and SDG 11 (Sustainable Cities and Communities) [64].