Prediction of groundwater drawdown and quality at Kyzylzharma well field in the arid Syr Darya Basin, Kazakhstan

Abstract

The geologically complex Kyzylzharma groundwater field is located in the south-central part of Kazakhstan in the lower Syr Darya Basin. It supplies the 243,000 population of Kyzylorda City by drinking and agricultural water needs. Numerical modeling was used to predict the consequences of increased groundwater withdrawal for future water supply needs. The results displayed a mean squared error for the groundwater simulations of about 0.6 m and was thus acceptable. The validated transmissivity was between 2 × 10⁻⁷ and 2 × 10⁻⁵ m/d. These parameters showed that the planned groundwater withdrawal will result in a depression cone reaching 90–100 m below present levels in 2040. Maximum groundwater level drawdown may reach 130 m below present levels. This drawdown increases risks for brackish saltwater intrusion into the main groundwater aquifer. The results point at the depletion of groundwater resources in the Syr Daria Artesian Basin and the risks for groundwater quality, particularly increase in mineralization. The key outcome is the recognition that effective joint management of the Syr Daria Artesian Basin’s groundwater resources, involving both Kazakhstan and Uzbekistan, necessitates development and operation of a joint numerical aquifer model. Building this model is a crucial tool for assessing the sustainable groundwater resources in the region.

Keywords:

• Sustainable groundwater
• Syr Darya Artesian Basin
• transboundary groundwater
• Uzbekistan

REVIEWING EDITOR:

• Sanjay Kumar Shukla, Edith Cowan University, Australia

SUBJECTS:

• Agriculture & Environmental Sciences; Earth Sciences; Civil, Environmental and Geotechnical Engineering

1. Introduction

Groundwater is becoming increasingly scarcer especially in arid areas (Zektser & Everett, 2001; Chen et al., 2010; Micklin, 2016). Therefore, it is increasingly important to estimate sustainable availability, quality of groundwater and risk in the long-term (Mukherjee et al., 2020; Rosenberry et al., 2015; Al-Katheeri, 2008; Tang et al., 2004). Current methods for groundwater risk include overlay index methods (Torkashvand et al., 2021), statistical methods (Su et al., 2019, 2020; Nsabimana et al., 2022) and process simulations (Geng et al., 2023; Machiwal et al., 2018). For process simulations, various types of software are used to model groundwater dynamics based on the groundwater flow equations (e.g. Konikow & Reilly, 1999; Condon et al., 2021; Malekzadeh et al., 2019). Groundwater transport is generally governed by changes in head described by a system of differential equations (Borevsky, 1988). By solving this system of equations, future changes in groundwater levels due changes in groundwater exploitation can be estimated. For this purpose, it is necessary to know the governing parameters and initial conditions (Musin & Khramchenkov, 2019). Numerical modeling of groundwater is thus used for solving various hydrogeological, hydroengineering and geo-ecological problems (Bundschuh & César Suárez, 2019; Bear & Cheng, 2010; Zhou & Li, 2011; Langevin et al., 2017; Dvoretsky et al., 2009). Similarly, different modeling techniques can be used to predict the future groundwater quality (e.g. Eid et al., 2023).

Kyzylorda region in the arid Syr Darya Basin, is one of Kazakhstan’s poorest and most vulnerable areas (Dukhovny et al., 2007; Lobanova et al., 2021). The poverty rates have more than tripled in the latest decades. About 800,000 people inhabiting the area are mainly working in agriculture. Core agricultural activities are irrigated rice production, cattle and fishery. The region is responsible for about 90% of Kazakhstan’s rice production (World Bank, 2019). The mean water consumption of the agricultural sector in Kyzylorda area is about 3.5–4 km³ per year corresponding to about 1/3 of the Syr Darya flow. Improved efficiency of the irrigation and conjunctive surface-groundwater use are the regional authoritýs priority.

The Syr Darya is one of the two main rivers in the endorheic Aral Sea Basin (Micklin & Williams, 1996; Pan et al., 2023). The other river is the Amu Darya. Syr Darya Basin encompasses the most productive agricultural regions in Central Asia (Ostrovsky, 2007; Moerlins et al., 2008). Kyzylorda, in the lower Syr Darya basin, is the largest city and the administrative center of the region with the same name. Kyzylzharma groundwater well field is in the central part of the region in the built-up area of Kyzylorda City and the only source of drinking water supply (Figure 1). It is confined to the Upper Cretaceous aquifer complex with the main operational Upper Turonian aquifer. The Syr Darya Artesian Basin constitutes a complex multi-layer pervious and semi-pervious flow system connected to the Aral Sea that is shared with Uzbekistan (Panichkin et al., 2017; Liu, Wang, et al., 2020). Complex hydrogeology together with lack of transboundary data and lacking long-term surface and groundwater observations make it difficult to model the system with reliability (Rafikov & Gulnora, 2014).

Figure 1. The experimental Kyzylzharma groundwater well field study area at Kyzylorda, Kazakhstan, in the Syr Darya Artesian Groundwater Basin.

In view of the above, the objective of this study was to set up, calibrate and validate a groundwater quantity and quality model that can be used to predict future groundwater levels and chemical content in the Kyzylzharma groundwater well field using solutions for stationary and transient states. This is necessary for the development sustainable future groundwater withdrawal that can help water managers and decision makers to improve the dwindling freshwater supply in the region. The developed groundwater simulation model is a first major tool for the management and rehabilitation of the downstream Aral Sea and the recognition that effective joint management of the Syr Daria Artesian Basin’s groundwater resources, involving both Kazakhstan and Uzbekistan, necessitates the development and operation of a joint numerical aquifer model. Building this model is crucial for assessing sustainable future groundwater resources in the region.

2. Materials and methods

2.1. Site description

The climate of the study area is classified as BWk according to Köppen (1936) and Geiger (Beck et al., 2018), characterized as dry arid middle-latitude desert. The mean annual rainfall is about 140 mm and mean monthly temperature varies between about 28 °C in July and −6 °C in January. The actual evapotranspiration is 700–1100 mm per year (Liu, Huang, et al., 2020).

The study area is a low-lying alluvial plain that crosses the Syr Darya Basin from southwest to northeast (Figure 1). The altitude varies from 130 to 136 m above mean sea level in the southeast to 112–116 m in the northwest. The general slope of the alluvial plain is to the northwest. The northeastern part of the study area is a flat surface with many takyrs, small ridges of well-fixed sand up to 3–4 m high, and remnants of dry riverbeds up to 1 m deep. The appearance of the plain changes in the central part of the area with alternating lakes, marshlands, salt marshes and croplands crossed by canals, irrigation ditches, and dry beds of ancient rivers. Sand dunes up to 2 m high are developed locally. The southwestern part of the territory is like the northeastern part of the area in terms of relief but distinguished by much smaller development of takyrs and a larger number of sand dunes. The height of sand dunes here reaches 5–6 m.

The Syr Darya in the area is characterized by strong meandering and oxbow lakes. The natural banks of the river are steep but not more than 0.7 m above the water surface in spring and not more than 3–4 m in low water. Numerous dams have been built along the riverbanks to protect the area from flooding. The Syr Darya has a flow rate of up to several hundred cubic meters per second. Until 1972, the river flow rate varied from 240 to 1700 m³/s. After the construction of the Shardara Reservoir and regulation of the river flow in 1974, its rate does normally not exceed 700 m³/s, and in the summer period (near Kyzylorda) it is reduced to less than 100 m³/s. The width of the Syr Darya is 140–600 m with a depth of 2.1–5.0 m. The flow velocity of the river is 0.8–1.7 m/s. The river freezes in December and opens from the ice in February–March. River water is characterized by a generally high mineral content and presence of microbes indicating influence from agricultural production waste discharged from countries upstream (RGP, 2023).

2.2. Hydrogeological setting and numerical modeling

To complete the input database for the model setup, the following information and data sources materials were used:

• Reports on exploration and assessment of operational groundwater reserves of the Cretaceous aquifer complex for the area corresponding to the Kyzylzharma groundwater well field (Syr Darya Artesian Basin, 1992);

• Results of previously performed numerical modeling on the regional scale for the Eastern Aral Sea hydrogeological conditions (Veselov & Panichkin, 2004; Panichkin & Miroshnichenko, 2014; Panichkin et al., 2017);

• Groundwater observations at wells drilled in the Cretaceous aquifer complex (pumping tests, well flowrate, groundwater levels and groundwater quality for the studied period).

Representative hydrogeological cross-section of the Kyzylzharma groundwater well field and conceptual schematization for the model simulations are shown in Table 1. The cross-section of the study area was schematized as five pressure-permeable layers separated by three weakly permeable ones (Table 1). Each permeable layer in the model corresponds to one layer in the section, and each weakly permeable layer corresponds to three refined layers to exclude disconnection of the grid cells. The regional Syrdarya Artesian Groundwater Basin contains five main water-bearing geological layered deposits corresponding to the Neogene–Quaternary, Upper Eocene, Upper Turonian-Senonian, Upper Albian-Cenomanian, and Lower-Middle Albian-Jurassic (Panichkin et al., 2017). The water-bearing geological layers are divided by three regional deposits of impermeable or semi-impermeable material the Paleogene, Lower Turonian and Lower-Middle Albian. These formations are distributed throughout the basin. The Upper Turonian aquifer represents the main groundwater bearing aquifer for the Kyzylzharma groundwater well field.

Table 1. Representative hydrogeological cross-section of the Kyzylzharma groundwater well field and conceptual schematization for the model simulations.

Item no	Name of aquifer and complex	Schematization of the model
		Type of layer	Layer number in the model*
1	Maastrichtian aquifer		2
2	Coniac-Campanian argillaceous weakly permeable complex		3–5
3	Coniac-Campanian aquifer		6
4	Upper Turonian aquifer		7
5	Lower Turonian argillaceous weakly permeable complex		8–10
6	Upper Albian-Cenomanian aquifer complex		11
7	Lower-Middle Albian argillaceous weakly permeable complex		12–14
8	Lower-Middle Albian aquifer complex		15

* Layer No. 1 of the model corresponds to impermeable Paleogene clays.

Hydrogeological conditions in the plane were schematized in the form of a rectangular modeling area measuring 61.3 × 50 km. Operational wells were specified in the model structure with time-varying flow rate boundary conditions. All existing and simulated wells operating the Upper Turonian aquifer were in the central part of the modeling area to minimize the influence of boundary conditions. The simulated area in the plane was approximated by an orthogonal irregular grid, equal to M × N = 143 × 122 columns and rows, respectively. Grid mesh size varied from 250 m in the central part of the study area to 750 m at the periphery of the model domain (Figure 2).

Figure 2. Grid design for simulations of the studied Kyzylzharma groundwater well field.

The GMS (Groundwater Modeling System) software was used to simulate the groundwater flow and water quality at the Kyzylzharma well field (e.g. Panichkin et al., 2017; Niswonger et al., 2011; Harbaugh, 2005). The software was developed by Brigham Young University, USA (GMS, 2021) but builds on the MODFLOW groundwater simulation code. The GMS is designed to model groundwater flow and transport of dissolved components in saturated and unsaturated zones including advanced tools for conceptual modeling, automatic parameter identification, construction of flow lines, stochastic modeling, geostatistics, visualization of modeling results in plane and section, construction of 3D stratigraphy models and generation of 2D and 3D networks.

Model calibration included calibration for steady and non-steady state conditions and post-processor integration to minimize convergence errors between calculated and observed values. The initial data preparation included the creation and digitalization of maps of the model top layers in absolute elevation, thickness of the separating layers, preliminary groundwater levels and groundwater salinity. For this purpose, the authors used the ArcGIS geographic information system. Using the abovementioned maps, we created layers for the model and converted them into special data formats, which were imported to the GMS modeling system. To perform the conceptual modeling, we used the 2D Grid, 3D Grid, 2D Scatter Point, Map, CAD Data and Image Data modules of the GMS software.

The purpose of the numerical model calibration was to define values of permeable and weakly permeable layers’ transmissivity. Steady state calibration was performed for observations in the beginning of 1974 to reveal hydrogeological conditions that existed for the conditionally undisturbed groundwater situation. The calculated groundwater levels were compared with the observed heads in wells during that period. Initial values of transmissivity for separating layers were obtained from pumping tests. The purpose of the non-steady calibration was to clarify the specific yield of water-bearing rocks and replicate the groundwater abstraction period from 1974 to 2010 from all boreholes divided into stress-periods (1 year) with an operational rate specified.

The results of the non-steady state calibration validated the convergence of the observed heads in 2010 collected from monitoring wells. The specific yield parameters were taken from the transient state calibration results. Calculation iterations were repeated until calculated and observed well heads matched satisfactorily. Based on this modeling experience, it was assumed that the permissible error in the convergence of observed and calculated results should not exceed 10–15% (Doherty & Hunt, 2010).

Three-dimensional steady-state transport of groundwater flow in heterogeneous and anisotropic strata is described by Bear and Cheng (2010): (1) $$\frac{\partial }{\partial x}\left(K_{\mathit{\text{xx}}}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{\mathit{\text{yy}}}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{\mathit{\text{zz}}}\frac{\partial h}{\partial z}\right)-W=\mathrm{0.\hspace{1em}}(x,y,z)\in \Omega \mathrm{ }$$(1) with initial conditions: $$\begin{array}{cc}h(x,y,z)=H_0(x,y,z)& \mathit{\text{on}}\hspace{0.17em}\Omega \end{array}$$ and boundary conditions: $$\begin{array}{cc}h(x,y,z)=H(x,y,z)& \mathit{\text{on}}\hspace{0.17em}{B}_1\\ K_n\frac{\partial h(x,y,z)}{\partial n}=Q(x,y,z)& \mathit{\text{on}}\hspace{0.17em}{B}_2\end{array}$$ where K_xx, K_yy and K_zz are hydrogeological conductivity along the direction of coordinate axes x, y and z (Lt⁻¹); h is head (L); W is volumetric flow per unit volume, representing water inflow and outflow (t⁻¹); S_s is specific storage of the porous media (L⁻¹). In general, the S_s, K_xx, K_yy, K_zz and W are function of spatial coordinates (S_s = S_s(x,y,z), K_xx(x,y,z), etc.). Ω is modeling domain, H₀(x,y,z) is head distribution for initial conditions; B₁ is constant head boundary, H(x,y,z,t) is head along B_1, B₂, the boundary with a given water flow rate, K_n is hydraulic conductivity along the normal to the boundary, B₂, Q(x,y,z) is groundwater withdrawal per unit volume; λ is a coefficient characterizing the relationship with the external groundwater transport, and h_B is groundwater level in the external environment.

Three-dimensional non-steady transport of groundwater flow in a heterogeneous and anisotropic medium is described by Bear and Cheng (2010): (2) $$\frac{\partial }{\partial x}\left(K_{\mathit{\text{xx}}}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{\mathit{\text{yy}}}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{\mathit{\text{zz}}}\frac{\partial h}{\partial z}\right)-W=S_s\frac{\partial h}{\partial t}$$(2) with initial conditions: $$\begin{array}{cc}h(x,y,z,t)=H_0(x,y,z)& \hspace{0.17em}\Omega ,\hspace{0.17em}t=0\end{array}$$ and boundary conditions: $$\begin{array}{cc}h(x,y,z,t)=H(x,y,z,t)& {\Gamma }_1,\hspace{0.17em}t\ge 0\\ K_n\frac{\partial h(x,y,z,t)}{\partial n}=Q(x,y,z,t)& {\Gamma }_2,\hspace{0.17em}t\ge 0\end{array}$$ where t is time.

The groundwater flow transport of dissolved chemical components is generally described by: (3) $$\frac{\partial (\theta C^k)}{\partial t}=\frac{\partial }{\partial x_i}\left(\theta D_{\mathit{\text{ij}}}\frac{\partial C^k}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(\theta {\nu }_i{C}^k\right)+q_s{C}_s^k$$(3) where C^k is concentration of the dissolved kth component (ML⁻³); θ is porosity of the transport medium (), x_i is coordinate corresponding to the ith Cartesian coordinate axis (L); D_ij is hydrodynamic dispersion coefficient tensor (L²T⁻¹); ʋ_i is transport velocity (LT⁻¹); q_s is volumetric flow per unit volume of aquifer, representing the source (positive) or sink (negative) (T⁻¹); and C_s^k is concentration of the dissolved k^th component in the source or effluent (ML⁻³).

2.3. Sampling and analyses

Water in the river, canals and irrigation ditches is usually low saline, but can reach 1.5–1.8 g/L (Schettler et al., 2013). Sulfate content can be up to 748 mg/L and the total hardness up to 14 mg/L, which makes it unsuitable for household drinking purposes without treatment. The groundwater is generally brackish in the Maastrichtian, Cognac-Campanian and Upper Albian-Cenomanian aquifer complexes. For this reason, it was important to combine the groundwater level drawdown simulations with groundwater quality simulations (Akramkhanov et al., 2011; Kulmatov, 2014). Groundwater is often brackish in Kazakhstan and previous studies have shown that groundwater is affected by salinity in the study area (Absametov et al., 2023; Zektser, 2006: Smolyar et al., 2012; Sagin et al., 2017). Based on these previous studies, the groundwater salinity of the Maastrichtian aquifer was set to a constant 4 g/L in the model setup. Initial values of groundwater salinity for Cognac-Campanian, Upper Turonian and Upper Albian-Cenomanian aquifer complexes to solve the concentration simulations were set based on borehole observations. In total 128 groundwater wells were sampled in 2020. Samples were taken from depths ranging from 425 to 575 m within the Upper-Turonian aquifer and analyzed regarding different chemical constituents. The chemical analyses generally followed the standards of APHA (1995) and DeZuane (1997). The groundwater samples were first filtered and then stored in acidified polyethylene containers. The pH and total dissolved solids (TDS) were quantified through a potentiometric technique, employing a fluid analyzer (Transmitter M-400 by Mettler Toledo, Columbus, OH). For the determination of Ca²⁺, Mg²⁺, HCO₃⁻ and CO₃²⁻, the titrimetric approach was utilized. The F⁻, Cl⁻, SO₄²⁻, NO₃⁻, NO₂⁻, NH₄⁺, I⁻, Br⁻ and K⁺ were analyzed using capillary electrophoresis (Capel-105M, Lumex Instruments, Mission B.C., Canada). The Na⁺, B⁺, Fe²⁺, Si⁴⁺, Li⁺, Sr²⁺, Cs⁺ and Rb⁺ were assessed by inductively coupled plasma spectrometry (ICPE-9800, Shimadzu Corporation, Kyoto, Japan). Groundwater chemical concentrations were compared to concentrations in samples from the Syr Darya. Graphical representations were drawn using Piper (1944) plots and Durov diagrams (Durov, 1948) using AquaChem Water Quality Analysis Software version 10.0 (AquaChem, 2021).

3. Results and discussion

A Piper plot was constructed to categorize the hydrochemical facies of groundwater (Figure 3). The cationic triangle revealed that 50.1% of the collected samples fell into the Cl-SO₄-Na category, while the remaining samples belonged to a non-dominant category. In the anionic triangle, Cl⁻ dominated in 93.3% of the samples, whereas the samples from the Syr Darya River exhibited a prevalence of Na⁺. Within the diamond shape Piper diagram, the groundwater samples from the Upper Turonian aquifer were classified as Cl-SO₄-Na water type with permanent hardness attributed to reverse ion exchange. The salinity in these samples stemmed from the notably high concentrations of chloride and calcium ions, particularly in samples distant from the Syr Darya River, which serves as the source of freshwater in the region.

Figure 3. Graphical representation of the hydro chemical facies and factors controlling the water chemistry: a) Durov diagram; b) Piper plot.

In addition, researchers in the field of hydrogeology have employed the Durov diagram, which connects pH, TDS and major ions, as a visualization tool (Durov, 1948; Figure 3(a)). The Durov diagram aids in elucidating three primary processes: mixing/dissolution, ion exchange and reverse ion exchange (Figure 3(a)). All samples, with a TDS exceeding 1000 mg/L, were positioned within the reverse ion exchange zone.

Anthropogenic pollutants such as the release of surplus irrigation water from agricultural fields and waste disposal (Srinivasamoorthy et al., 2008; Jacks et al., 1999), as well as the deposition of atmospheric chloride (Cl⁻), can lead to increased levels of chloride in the environment (Biswas et al., 2012). When the ratio of sodium (Na⁺) to chloride (Cl⁻) exceeds one, the weathering of silicate minerals may emerge as a notable contributing factor (Meybeck, 1987).

The Syr Darya water was characterized by a generally high mineral content and presence of microbes indicating influence from wastewater (Table 2). The observed concentrations were compared to the WHO guidelines for drinking water. However, except for a high total microbial count (TBN), the river water did not display any noteworthy pollutants.

Table 2. Syr Darya water quality.

No	Chemical indicator, unit	Concentration	WHO (2017) permissible drinking water limit
1	Temperature at the time of sampling,^oC	3 to −27	–
2	pH	7.2–8.7	6.5–8.5
3	Odor at 20^oC and 60^oC, points	1	≤2
4	Turbidity, mg/L	7.6–52.8	5–10
5	Color of water, degree	20.8–152.5	5 TCU
6	Ammonium NH4, mg/L	0.08–1.5	–
7	Chloride Cl^-, mg/L	80–100	250
8	Nitrite NO2, mg/L	0.04–0.16	3
9	Nitrate NO3 mg/L	10.8–13.0	50
10	Sodium Na^+, mg/L	150–170	200
11	Chloride Cl, mg/L	73.9–104.5	200–300
12	Iron Fe, mg/L	0.16–1.75	0.3
13	Total hardness, mg/L	8.5–14.0	10–500
14	Sulfate SO4, mg/L	382–748	250–500
15	General mineralization (dry residue), mg/L	810–1840	–
16	Aluminum Al, mg/L	– 0.23	0.9
17	Alkalinity, mol-eq/L	3.5– 4.0	–
18	Oxidation, mg О2/L	1.36–3.52	–
19	Fluoride F, mg/L	0.24–0.67	1.5
20	Manganese Mn, mg/L	– 0.1	–
21	Copper Cu, mg/L	– 0.048	2
22	Polyphosphate, mg/L	0.14–0.61	1
23	Free carbon dioxide CO2, mg/L	– 32.2	–
24	BOD5 (biochemical oxygen demand), mg О2/L	0.39–1.95	5
25	Dissolved oxygen, mg O2/L	2.75–3.98	4–6
26	Suspended solids, mg/L	4.0–108	30
27	Total dissolved solids TDS, mg/L	90–1300	–
28	Total α-radioactivity, Bq/L	0.09	0.5
29	Total β-radioactivity Bq/L	0.1	1
30	Coliform index	68–702	10 coliform organisms per 100 mL
31	TBN (total microbial count)	45	100/mL to 500/mL

Schematic maps of the transmissivity of separating layers (Coniac-Campanian, Lower Turonian and Middle Lower Albian argillaceous sediments), were created according to the results of the non-steady state calibration and shown in Figure 4. Minimum hydraulic conductivity for the Cognac-Campanian argillaceous displays a weak permeability equal to 2·10⁻⁷ - 5·10⁻⁶ m/d. For the Lower Turonian complex values of hydraulic conductivity varied from 4·10⁻⁶ to 1·10⁻⁵ m/d. For the layer of Lower-Middle Albian complex hydraulic conductivity was in the range of 5·10⁻⁶ to 2·10⁻⁵ m/d. For the modeled area there was a tendency for hydraulic conductivity to increase from north-west to south-east. The graph of modeled and measured heads (Figures 5 and 6) shows that the accuracy of the steady state calibration agreed in general with observations.

Figure 4. Transmissivity of the different weakly permeable geologic layers.

Figure 5. Comparison between calculated and observed heads for steady-state conditions (time step = 1 day, 1974).

Figure 6. Calculated groundwater levels of the Upper Turonian aquifer from the steady state calibration (1974).

Calculated long-term groundwater balances for steady state conditions in 1974 are shown in Table 3. This period was chosen as it represents almost natural conditions with little human influence.

Table 3. Long-term groundwater balances based on steady state calibration for 01/01/1974.

Aquifer and geologic complex	Recharge, m^3/day				Discharge, m^3/day
	Inflow through side boundaries	Inflow from layers above	Inflow from layers below	Total inflow to the aquifer	Outflow through the side boundaries	Outflow from layers above	Outflow from layers below	Total discharge of the aquifer (complex)
Maastrichtian	808.68	0	0	808.68	−697.48	0	−111.16	−808.64
Coniac-Campanian	1785.39	111.15	164.63	2061.17	−1921.65	0.00	−139.49	−2061.14
Upper Turonian	2629.48	139.49	404.80	3173.77	−2987.45	−164.63	−21.64	−3173.72
Upper Albian-Cenomanian	3394.32	21.63	99.95	3515.90	−3075.90	−404.80	−35.06	−3515.77
Middle Lower Albian	769.45	35.08	0	804.54	−704.54	−99.97	0	−804.51
Total for the simulated area	9387.32				−9387.03

The calculated groundwater levels were compared with observations from the monitoring wells. The values of the specific yield were adjusted in this calibration step. Initial values were obtained from the pumping tests. Iterations were repeated until a satisfactory match between the calculated and measured groundwater levels in groundwater wells were achieved. Based on the modeling, it was assumed that the acceptable error in solving the problem regarding head values should not exceed 10–15%.

The calibrated specific yield coefficient for the rock formations were as follows: 7∙10⁻⁶ L/m for the Maastrichtian aquifer horizon, 1∙10⁻⁵ L/m for the Cognac-Campanian aquifer horizon, 1∙10⁻⁵ L/m for the Upper Turonian aquifer horizon, 1.1∙10⁻⁶ L/m for the Upper Albian-Senonian aquifer complex, and 4∙10⁻⁶ L/m for the Lower-Middle Albian aquifer complex.

The aim of solving the inverse non-steady state conditions was to calibrate the specific yield values. The model reproduced the changes in hydrogeological conditions from 1974 (undisturbed period) to 1 October 2010. The time step of the inverse non-stationary solution was 1 year. Water extraction from production wells was carried out from the Cognac-Campanian, Upper Turonian and Upper Albian-Senonian horizons during this period (Figure 7).

Figure 7. Calculated non-steady state groundwater levels for the Upper Turonian aquifer in 2010.

The calibrated specific yield obtained from the model are consistent with experimental data and the results of previous studies (Trushel, 1993; Panichkin et al., 2011). As seen from Figures 8 and 9, simulated groundwater levels coincided reasonably well with representative observations. Thus, it can be concluded that the accuracy of solving the inverse non-stationary problem met the requirements for the model of the Kyzylzharma well field. The model adequately reflected the hydrogeological conditions that existed during the undisturbed period (1974), as well as the changes that occurred because of groundwater exploitation from 1974 to 2010. The model can thus be used for predictive purposes.

Figure 8. Observed and simulated groundwater level for well no. 351.

Figure 9. Observed and simulated groundwater levels for well No. 2561.

The calibrated model was used to predict the groundwater situation near to 2040, if 113 existing and 67 projected wells were used for exploiting the Upper-Turonian aquifer. Total productivity of all wells would be in the order of 100,000 m³/day. This situation is shown in Figures 10 and 11.

Figure 10. Prediction of groundwater levels in about 2040 in the Upper Turonian aquifer.

Figure 11. Drawdown of groundwater level in about 2040 compared to 2010 in the Upper Turonian aquifer.

According to Figures 10 and 11, we can conclude that at the end of the forecast period near to 2040, the increased groundwater withdrawal in the Upper Turonian aquifer will form a depression cone, covering the entire study area. Compared to 2010, the levels may decrease by 90–100 m in the center of the depression cone. Along the boundaries of the simulated area, levels may decrease by 4–10 m. The maximum groundwater table drawdown in the production wells may reach 130 m. This is, however, still less than the permissible maximum drawdown corresponding to 150 m set by the regional authorities.

Groundwater quality changes in about 2040 were predicted using the described simulation approach. Simulations of water quality were performed for the main Upper Turonian aquifer in the Kyzylzharma well field using the MT3D module of the GMS modeling system. By future groundwater table drawdown, there are risks for saline water intrusion from the brackish groundwater in the Maastrichtian, Cognac-Campanian, and Upper Albian-Cenomanian aquifer complexes. A predicted groundwater table drawdown of about 130 m in 2040 may result in the highest groundwater salinity equal to 1.33 g/L. This brackish water may reach the wells with the highest flow rate and the maximum drawdown (Figure 12).

Figure 12. Predicted mineralization of groundwater (total dissolved solids) in about 2040.

Thus, it is likely that overexploitation of groundwater from the Upper Turonian aquifer of the Kyzylzharma well field will result in groundwater quality deterioration (see also e.g. MacDonald et al., 2016). The regional authorities have set a maximum limit to the acceptable total mineralization rate of the groundwater equal to 1.5 g/L in about 2040. Similarly, the sulfate concentration in the groundwater is likely to increase due to the groundwater drawdown. The simulations indicated that the maximum sulfate concentration in the groundwater may reach 0.40 g/L (Figure 13). The permissible maximum sulfate concentration has been set to 0.50 g/L by the authorities.

Figure 13. Predicted increase in sulfate concentration in groundwater in around 2040.

As seen from the above simulation results, future water demand can be satisfied by groundwater exploitation of the Upper Turonian aquifer complex of Kyzylzharma well field. This will, however, result in a substantial drawdown of the groundwater table in the central parts of Kyzylorda City and as well affect the border of the urban settlements. Consequently, this may induce risks for subsidence of the ground surface. Also, the current plan for groundwater exploitation is likely to result in groundwater quality deterioration. Thus, it is necessary to implement a set of measures to continuously monitor water withdrawal and salinity in the well field. The monitoring results can be continuously included in the present numerical groundwater model for further updating of prediction results and evaluating alternative solutions and more efficient groundwater management of the aquifer.

Central Asia is a region known for its water scarcity (Chynybaeva, 2023; Wang et al., 2022; Satke, 2021), primarily stemming from the management of surface water resources. This challenge became more pronounced following the emergence of independent states in 1991 (Wikipedia, 2023). Particularly, the downstream nations soon recognized the risk for substantial water losses due to the absence of agreements regarding shared rivers. This realization prompted them to initiate negotiations. The successful conclusion of multiple agreements concerning transboundary surface water resources in Central Asia is a significant achievement in the context of attaining SDG 6.5.2 (UN, 2023).

Water resources play a key role in the economies of countries of the arid region in Kazakhstan. The Central Asia Climate Change Conference in partnership with the World Bank posed the question, ‘What is the extent of the water crisis in Central Asia and what are the possible solutions for the communities living in those areas’? (CACCC, 2023). To answer this question, it is necessary to perform predictions. This article proposes one of the possible approaches for this purpose by conducting numerical modeling and finding prognostic scenarios for the well fields in an area characterized by water scarcity and frequent periods of drought. It is, however, not the only solution. Increasing population combined with increasing water scarcity due to climate change requires a stronger emphasis on demand management. In the case of Kyzylorda groundwater resources, increased exploitation must be balanced by economic return and risks for groundwater quality deterioration and the overall sustainability of the region’s natural resources. The Kyzylzharma well field is a part of the extensive Syr Darya Transboundary Artesian Basin (Liu, Wang, et al., 2020; Kokimova, 2019). Joint management and protection of transboundary groundwater resources are currently impossible without a groundwater agreement between Kazakhstan and Uzbekistan. This is the main transboundary problem on groundwater in this region. Only the Kazakhstani side has agreed to conduct a survey of the groundwater situation. In line with this, our results revealed two primary issues: the depletion of groundwater resources in the Syr Daria Artesian Basin and the potential decline in groundwater quality, particularly an increase in mineralization. Thus, the key outcome of this study is the recognition that effective joint management of the Syr Daria Artesian Basin’s groundwater resources, involving both Kazakhstan and Uzbekistan, necessitates the development and operation of a joint numerical aquifer model. Building this model is a crucial tool for assessing the available groundwater resources in the region. To facilitate this, it is imperative to:

• Transit the forthcoming numerical model into a continuously operating model in both neighboring states, Kazakhstan and Uzbekistan.

• Foster international cooperation for joint, optimized management of the aquifer’s groundwater resources based on agreed scenarios and the exchange of hydrogeological monitoring data.

• Implement ongoing monitoring of groundwater, encompassing piezometric levels, flow rates, and groundwater quality in all operational wells, regardless of ownership or purpose. Regularly assess the technical and environmental condition of water intake wells.

• Enhance national legislation regarding mandatory groundwater monitoring.

• Regulate water intake well flow rates in strict accordance with approved and agreed upon exploitable resource values for both countries.

• Modernize the groundwater monitoring system by installing state-of-the-art equipment for recording well flow rates and groundwater levels. Implement data quality control measures in line with international standards. Develop groundwater quality monitoring programs covering the entire aquifer.

• Promote international cooperation on groundwater quality issues between neighboring countries.

The groundwater is the only source of drinking water supply for Kyzylorda City and surrounding settlements. However, it also serves the thirsty agriculture and the main rice production areas in Kazakhstan. As such, the groundwater needs to be prioritized and protected for drinking water needs. In view of this, there are major reasons to manage the demand side for agricultural water use and look at ways to improve the efficiency of agricultural water use. As rice requires a lot of water, it is also necessary to try to resort to rice varieties and crops that are more water efficient.

5. Conclusions

• The results of this study have improved the knowledge on the complex geological formation of the lower Syr Darya Artesian Groundwater Basin. It is especially important to determine the degree of transmissivity and transport of groundwater between the different separating layers (Coniac-Campanian, Lower Turonian and Middle Lower Albian argillaceous sediments) constituting the main aquifer. This knowledge can be used to better understand how to use the areás groundwater in a sustainable manner.

• The mean squared error for the steady-state groundwater table simulations of the Maastrichtian aquifer was equal to 0.80 m, for the Coniac-Campanian 0.63 m, for the Upper Turonian 0.41 m, for the Upper Albian-Cenomanian 0.37 m, and for the Middle Lower Albian 0.61 m. It can, thus, be stated that the error in the results was tolerable.

• The calibrated transmissivity of the Maastrichtian, Coniac-Campanian, Upper Turonian aquifers, Upper Albian-Cenomanian, and Middle Lower Albian aquifer complexes agreed well with observations and previous studies (Veselov & Panichkin, 2004; Panichkin & Miroshnichenko, 2014). Minimum transmissivity (2·10⁻⁷ – 5·10⁻⁶ m/day) was found for the layer of Coniac-Campanian argillaceous sediments. The transmissivity values were higher for the Lower Turonian argillaceous layer varying from 4·10⁻⁶ to 1·10⁻⁵ m/day. For the Middle Lower Albian argillaceous sediments layer, the transmissivity ranged from 5·10⁻⁶ to 2·10⁻⁵ m/day. The results of the modeling revealed a tendency for the transmissivity of the argillaceous sediments to increase from northwest to southeast. Simulated transmissivity values agreed well with previous studies (Trushel, 1993; Panichkin et al., 2011).

• The model values of the specific yield varied from 1.1·10⁻⁶ for the Upper Albian-Cenomanian aquifer complex to 1·10⁻⁵ for the Coniac-Campanian and Upper Turonian aquifers. For the Middle Lower Albian aquifer complex, the specific yield was equal to 4·10⁻⁶. In the non-steady state calibration, groundwater levels agreed well with observations during 1974 to 2010 and previous studies (Pan et al., 2023; Shapiro et al., 1992).

• The results show that the planned groundwater withdrawal will result in a depression cone reaching 90-100 m below present levels around 2040. Along the boundaries of the depression cone, a groundwater decrease of 4–10 m may be expected. Maximum groundwater drawdown in the production well may reach 130 m below present levels. The groundwater table drawdown increases risks for brackish saltwater and sulphate rich water intrusion into the main groundwater aquifer.

• The water availability in the Kazakhstani part of the aquifer is approximately ten times lower than in the Uzbek part. Given the substantial pressure on resources in Kazakhstan, particularly in regions where it is the sole source of fresh water, and considering rapid population growth, it is advisable to implement a series of national-level demand-management measures in Kazakhstan.

• Overall, the continued development of the model relies on data availability. Addressing this challenge necessitates the establishment of cooperation processes and trust-building between Kazakhstan and Uzbekistan. Data and information exchange will play a vital role in the secure management of the transboundary groundwater aquifer.

• The developed model of the Kyzyljarminskoye groundwater field serves as a practical example and a viable avenue for further groundwater resource management in the arid Syr Darya region and the transboundary Syr Darya Artesian Groundwater Basin.

Author contributions

Conceptualization, Absametov M. and Murtazin Ye.; methodology, Ospanov K.; software and validation, Murtazin Ye.; formal analysis, Absametov M.; resources, Rakhimova V., and Berndtsson R.; data curation, Ospanov K.; writing original draft preparation, Rakhimova V.; writing - review and editing, Berndtsson R.; visualization, Ospanov K.; supervision, Berndtsson R.; project administration, Murtazin Ye.; funding acquisition, Absametov M. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

We appreciate the support of the Ministry of Science and Higher Education of the Republic of Kazakhstan.

Disclosure statement

The authors declare no conflict of interest.

Data availability statement

The data used in this study are available upon request from the corresponding authors.

Additional information

Funding

This research was funded by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. BR10965134 ‘Assessment of fresh groundwater resources as the main source and long-term reserve of sustainable drinking water supply for the population of the Republic of Kazakhstan’).