Hydrogeology and geochemistry of a tectonically controlled, deep-seated and semi-fossil aquifer in the Zambezi Region (Namibia)

Recent exploration has revealed that deep-seated and large groundwater reservoirs in Africa’s intracontinental basins can be regarded as an additional strategic resource for supply of drinking water. The origin, genesis and recharge of these groundwater reservoirs, however, are still poorly understood. A multidisciplinary approach involving remote sensing, geophysical surveys and hydraulic investigations, as well as hydrochemical and isotope studies, was pursued to gain better insight into the genesis and the potential of a recently discovered lower Kalahari aquifer (LKA) located in the Zambezi Region (Namibia). The study shows that regional tectonic activity associated with the propagation of the Okavango Rift Zone had a tremendous impact on the drainage evolution and hydrogeological setting of the region. Furthermore, there is geomorphological evidence that the LKA—prior to tectonic subsidence and burial—was part of a paleochannel of the upper Zambezi River. Hydraulic continuity could be confirmed by geochemical evolution down the flow path. Cation exchange combined with dissolution of calcite progressively produces alkalinity and sodium and consumes calcium in the north–south direction. Comparison of stable isotope content of the LKA with modern rainfall indicates that the recharge occurred under cooler climate conditions. Analysis of 14C concentrations and 36Cl/Cl ratios show that the age of the groundwater exceeds 100 ka and is hence older than presumed. It is concluded that the assessment of the sedimentology, tectonic structures and geochemistry are key factors for understanding both the paleoclimatic and modern recharge processes of deep-seated aquifer systems.


Introduction
Intracontinental basins often comprise significant and partially unexplored fossil or semi-fossil aquifers. Due to their frequently considerable size, large amounts of freshwater are stored in these deep-seated sedimentary aquifers. These aquifers provide a valuable additional source of drinking water, which may be used for industrial and agricultural purposes. Hence, deep-seated aquifers are of increasing relevance to water supply in arid and semi-arid regions of the world. The groundwater, however, is often of considerable age, and the hydrogeology, as well as the exploration potential and recharge processes of such deep-lying groundwater bodies, is poorly understood. Consequently, any ill-considered groundwater extraction can rapidly lead to permanent depletion of aquifer reserves, i.e. groundwater mining (Margat et al. 2006).
In northern Africa, scientific attention was drawn toward aquifers in extensive basins bordering the Sahara Desert to the north, east and south, including the North-Western Sahara Aquifer System (Sahara and Sahel Observatory (OSS) 2008b; Siegfried 2004), the Nubian Aquifer (Sultan et al. 2007(Sultan et al. , 2011 and groundwater of the Taoudenit (Fontes et al. 1991), and the Illumeden (Sahara and Sahel Observatory (OSS) 2008a) and Lake Chad basins (Bouchez 2015). In southern Africa, the Kalahari Basin was assumed to host several deep-lying sedimentary aquifers. This basin stretches over some 2200 km in a north-south direction from the SE tip of the Democratic Republic of Congo to South Africa, covering a total area of approximately 1.6 million km 2 (Haddon and McCarthy 2005). In the Northern Kalahari, here defined as the area to the north of the Proterozoic Ghanzi-Chobe orogenic belt (Fig. 1a), sedimentary depocenters with thickness exceeding 400 m are developed in the border area between Angola and Namibia and in the Okavango Rift Zone (ORZ) (Haddon 2005;Miller et al. 2010). In these areas, two deep-seated freshwater aquifers were discovered in recent years: the Ohangwena II aquifer (Fig. 1c), which is developed within a paleo-fan of the Cubango River (Lindenmaier et al. 2014;Wallner et al. 2017), and the Lower Kalahari Aquifer (LKA) in the Zambezi Region (ZR) of Namibia, which is the subject of this study. Recent geophysical investigations provide evidence that the Okavango Delta in Botswana is also underlain by a sedimentary freshwater aquifer (Podgorski et al. 2013(Podgorski et al. , 2015. Due to its sedimentological and tectonic position and the prevailing low hydraulic gradients in these basins, the groundwater is considered semi-fossil. This means that the infiltration of water occurred a considerable time ago under wetter climatic conditions (pluvials), but that the aquifer is probably still connected to the hydraulic cycle and replenished to some extent.
The LKA is situated in the eastern parts of the Zambezi Region between 17°30′S and 18°30′S and between 23°15′E and 24°30′E, and is bordered to the east and south by the Kwando and Linyanti rivers, respectively (Fig. 1). Based on exploratory boreholes drilled between 2003 and 2007, the LKA is found at depths from 130 to > 250 m as part of the sedimentary succession of the Kalahari Supergroup (Margane and Beukes 2007). It was recognized during the early stages of investigation that the development of the aquifer system in this area was affected by tectonic subsidence and uplift (Margane et al. 2005). The nature of past and recent recharge, however, remained unclear. The emphasis of the actual study is on the assessment of sedimentology and tectonic evolutions, as they are considered crucial for understanding both the paleoclimatic and the recent recharge processes. An interdisciplinary, stepwise approach was developed involving (i) remote sensing data for the detection of characteristic landforms and lineaments, (ii) geophysical surveys to define the aquifer structure, and (iii) hydraulic, (iv) hydrochemical and (v) isotope studies to assess and characterize hydraulic interactions and flow processes. The study includes a comprehensive review of the literature on the paleoclimatological and geological conditions in the Kalahari Basin.
In the early 1990s, an extensive drilling program for rural water supply was realized (Margane et al. 2005), and subsequent groundwater investigations were carried out at the beginning of the millennium under a bilateral technical cooperation project between Namibia and Germany. All these activities led to a fair amount of reliable hydrogeological information available for the ZR. However, compared to the adjacent Okavango Basin, very few international studies exist in this area.

Hydrology and climate
The catchment of the Kwando-Linyanti-Chobe system, whose headwaters lie in the Angolan highlands, covers an area of over 156,000 km 2 and is situated between 13°S and 19°S and between 19°E and 14°E (Fig. 1a). In the case of rare high floods, the Kwando River is connected to the adjacent endorheic drainage system of the Okavango River via the Magweqgana Spillway (Fig. 2). The source of the Kwando River is located at an elevation of about 1350 m above sea level (a.s.l.). The lowest parts of the catchment are formed by the Linyanti-Chobe floodplains, which lie at about 930-940 m a.s.l., and the adjacent endorheic Mababe (917-930 m a.s.l.) and the Makgadikgadi (900-930 m a.s.l.) depressions. Compared to neighboring rivers, the Kwando at the gauging station in Kongola is characterized by remarkably small average discharge of 31.9 m 3 /s, corresponding to 9 mm/a, and relatively low seasonal variation. This is indicative of high evapo-transpiration losses from widespread meanders, oxbow lakes, swamps and floodplains, and of losses during infiltration to the underlying shallow groundwater. This assumption is consistent with measurements of the average annual temperature and class A pan evaporation at the Katima Mulilo (17°30′S, 24°16′E, 946 m a.s.l.) and Sesheke stations (17°28′S, 24°18′E, 944 m a.s.l.) located nearest the study area (Fig. 2), which vary between 21.3 and 21.9°C and between 2215 and 2506 mm/a, respectively (Margane et al. 2005;Yachiyo Engineering Co. Ltd. 1995).
The climate of the ZR is referred to as a hot arid steppe climate (BSh) by the Köppen-Geiger system (Peel et al. 2007) and is characterized by cool and dry winters, with minimum temperatures during June and July, and a hot rainy season during the summer months from mid-November to March. The rainfall regime is referred to by Gasse et al. (2008) as uni-modal austral summer rainfall. Due to the convective nature of rainfall events, thunderstorms are common, and the rainfall regime shows strong variations from year to year as well as spatial variability. The ZR marks a transition zone between semi-humid and semi-arid conditions, as regional Fig. 1 a Map showing the study area, the position of the Okavango Rift Zone and the presumed maximum extent of paleolakes in the Makgadikgadi Basin. LKA = Lower Kalahari Aquifer. b Surface topography (vertical exaggeration ≈ 500) and conceptual geological cross section (vertical axis strongly exaggerated and not to scale). c Inlay map shows the approximate extent of the Kalahari Basin (after Haddon 2005) and the Kwando-Chobe Catchment. OKZA = Okavango-Kalahari-Zimbabwe Axis; GTA = Griqualand-Transvaal Axis; OhR, OmR and ZR = Ohangwena, Omaheke and Zambezi regions in Namibia, respectively precipitation in southwestern Africa is controlled by the seasonal shift of the Intertropical Convergence Zone (ITCZ) and the development of the Congo Air Boundary (CAB). The former marks the surface air boundary between the northeast and (dry) southeast trade winds, the latter between low-level tropical easterlies from the Indian Ocean and westerly unstable and moist air flows from the Congo Basin (e.g. Gasse et al. 2008;Nicholson 2011;Reason et al. 2006). Rainfall in the Northern Kalahari is hence associated with a complex interplay of airstreams originating from the warm western Indian Ocean and the South Atlantic off the northern Angolan and Congolese coastline. In contrast, the Atlantic coastline south of approximately 16°S is influenced by the Benguela Upwelling system of cold water, which stabilizes air and prevents convection throughout the year. The cold sea surface temperatures are therefore the reason for the pronounced aridity along the Namibian and southern Angolan coastlines.
The highest mean annual rainfall (MAR) in the region is observed in the Angolan highlands, with maximum values for the upper Okavango and Kwando river catchments of about 1200-1400 mm/a (Pombo et al. 2015;Steudel et al. 2013), whilst rainfall in the lowest areas of the basin-the Makgadikgadi Plains-is only 300-400 mm/a (Burrough et al. 2009;Robertson 2004). The MAR at the meteorological stations at Katima Mulilo and Sesheke is 650-700 mm, with 90% occurring during austral summer.

Geology and tectonics
Tectonically, the ZR forms part of the ORZ which has been described as an incipient or nascent continental rift basin (Modisi 2000;Modisi et al. 2000) that underwent early stages of continental extension. Continued seismicity in the area provides evidence for ongoing weak tectonic activity in the region (Hutchins et al. 1976;McCarthy 2013;Reeves 1972), which is linked to the propagation of the East African Rift System (EARS) (Bufford et al. 2012). The ORZ forms part of the southwestern branch of the EARS , and is defined by three characteristic Ben échelon^northeasttrending and probably syntectonic depocenters that correspond to about 100-km-long and 40-80-km-wide structural graben systems, namely the Ngami, Mababe and the Linyanti-Chobe depressions (Kinabo et al. 2007). In total, this complex graben and horst structure is about 450 km long and  (Fig. 1a). In the northeastern direction, the ORZ borders the Machile Basin, which also represents a graben structure and once formed the valley of the proto-Kafue River (Moore et al. 2013). Wanke (2005) determined from Landsat TM images that the Eiseb Graben in the Omaheke Region ( Fig. 1c) in eastern Namibia represents a southeastward extension of the ORZ.
The depressions are filled with predominantly fluviatile-lacustrine deposits of the Kalahari Supergroup. The Kalahari Supergroup is an Upper Cretaceous to recent sedimentary succession consisting of terrestrial, siliciclastic, unconsolidated to semi-consolidated sediments overlying basement rocks. The sediments are composed mainly of conglomerates, gravels, clays, sandstones and unconsolidated sands. The base of the Kalahari corresponds to the Cretaceous African (Erosion) Surface, which was described by Miller (2008) as a Bdeeply dissected pre-Kalahari terrain^. The interior basin developed in response to Bdownwarp^and epeirogenic uplift along several axes including the Okavango-Kalahari-Zimbabwe Axis (OKZA) during the Upper Cretaceous/Early Tertiary (de Wit 2007;Haddon and McCarthy 2005). Epeirogenic uplift alternated with extended episodes of tectonic stability and erosion. Hence, deposition was intermittent, with periods of non-deposition corresponding to pedogenetic phases and phases of pedogenetic cementation (Miller 2008). It is generally assumed that the unconsolidated sand in the Kalahari Basin was formed either by fluvial or aeolian transport of material or by weathering of in situ sandstones, and that both processes occurred simultaneously (Gärtner et al. 2013;Haddon 2005).
The Kalahari Group is commonly underlain by basalts of the Karoo Supergroup which are allocated to the Batoka, Stormberg or Rundu formations in Zambia, Botswana/Zimbabwe and Namibia, respectively. The Karoo basalts erupted mainly between 184 and 178 Ma (Jourdan et al. 2005(Jourdan et al. , 2007Jones et al. 2001 and ref. therein). Money (1972) and Ridgway and Money (1981) describe the Batoka Basalt in western Zambia as an often olivine-rich tholeiite. In the Barotse Basin of Zambia, the thickness of the basalts increases in the western direction, from 30 m at the eastern margin to a maximum observed thickness of 390 m in the south along the Zambezi River (Unrug 1987). In Botswana, the basalts are up to 1000 m thick (Key and Ayres 2000). Surface outcrops of basalt are present over a large area in the Zambezi valley upstream and below the Victoria Falls. Smaller outcrops form the rapids in the upper course of the Zambezi River, e.g. Ngonye Falls and Mambova Rapids (Fig. 2). Interpretation of aeromagnetic data shows that basalt suboutcrops are widespread throughout northeastern Namibia (Fig. 16.4a in Miller 2008) as well as in northern Botswana (Key and Mothibi 1999).

Hydrogeology
A continuous aquifer, referred to as the Upper Kalahari Aquifer (UKA), which is characterized by Quaternary unconsolidated sand and silt of the Kalahari Supergroup, covers the entire area of the ORZ. The mainly alluvial deposits can be divided into recent and older alluvium (Thomas and Shaw 1991). The recent alluvium is seasonally inundated by the Okavango, Kwando, Linyanti, Chobe and Zambezi floodwaters and represents the riverine floodplains of the Okavango Delta, the Kwando River, the Linyanti-Chobe swamps and the Eastern Zambezi Floodplains (Fig. 2). Older alluvium covering the remaining areas, in contrast, is characterized by virtual absence of modern surface discharge. Areas to the north and northwest of the ORZ, i.e. west of the Kwando River, are largely covered by ancient linear dune fields (O'Connor and Thomas 1999;Thomas et al. 2000), being part of the large Northern Kalahari Dune Field (Stokes et al. 1998). In the ZR, elevated areas to the north of the Kongola-Katima Mulilo Highway (Fig. 2) are covered by alluvial sand or reworked aeolian sand. The vegetation is dominated by Kalahari, miombo and teak woodlands, which are deep-rooting trees (Fanshawe 2010). The older alluvium to the south is entirely different and has been described as abandoned inland fans that are covered by loamy soils and dominated by mopane woodland (Baars 1996;Mendelsohn and Roberts 1997). Mopane trees are essentially shallow rooters and are known to be more salt-tolerant than most other savanna tree species. They are therefore indicative of higher clay content and poorly drained and often alkaline soils.
The LKA was discovered by means of exploratory drilling carried out between 2004 and 2007. Five boreholes were drilled into the LKA. The top of the LKA was identified at depths between 125 and 150 m (Margane et al. 2005;Margane and Beukes 2007). In the two northern locations (boreholes no. 041002 and 041004), the Karoo basalts were found to underlie the Kalahari at depths between 193 and 222 m. For the two boreholes no. 041006 and 200254 (Fig. 3) drilled further south, the base of the Kalahari was not reached at the final drilling depth of 189 m and 250 m, respectively. The sediments of the LKA are made of fine-to medium-grained sandstones with interbedded layers of clay and calcareous or silicified sandstone. The LKA appears to be hydraulically separated from the upper aquifer by a bluish-green layer of clay or clayey silt with sand, which in large parts of the area forms a continuous, low-permeable layer. Between 2007 and 2015, 16 additional boreholes with depths varying from 130 to 250 m were established to provide water for small settlements or cattle stations. Unfortunately, detailed lithological and hydraulic records for these boreholes are not available. In contrast to the UKA, the groundwater of the LKA is mostly fresh.

Paleo-drainage conditions
As the Kalahari Basin undergoes tectonic deformation, there is a strong structural control on its drainage network. It is nowadays generally assumed that the Limpopo River was the main river draining the central interior areas of southern Africa, and that the Okavango, Kwando, upper Zambezi, Kafue and Luangwa rivers formed part of this enormous paleo-drainage basin (Moore 2004;Moore and Larkin 2001). Flow of these northern tributaries into the Limpopo River was severed by downwarp of the Kalahari Basin and epeirogenic uplift along the OKZA and the Griqualand-Transvaal Axis, located further south. The OKZA thus broadly forms the divide between the endorheic parts of the northern Kalahari Basin and the Limpopo and Orange rivers. The timing of the uplift and formation of the endorheic interior drainage basin is uncertain, but the Late Cretaceous (c. 65 Ma) has been suggested as an upper constraint. A late Paleogene age (c. 25 Ma) may, however, be much more realistic. Several distinct phases of rejuvenation of uplift along these axes occurred during the Late Oligocene (c. 24-20 Ma), Late Miocene (c. 10 Ma) and Pliocene (Haddon 2005;Walford et al. 2005).
Abrupt changes in river courses are numerous in the area and are often interpreted as capture elbows formed by tectonic movements, changes in base levels and subsequent headward erosion. Tectonic movement also led to the formation of paleolakes of varying size covering the Makgadikgadi Basin and adjacent areas such as today's Okavango Delta and the ZR. Nugent (1990) recognized that the Upper Zambezi and the Mid-Zambezi system, which is located downstream of Victoria Falls, developed as separate rivers until comparatively recent times. Cotterill (2006) and subsequently Moore et al. (2013)   was diverted into the Makgadikgadi Basin. Basalt bars exposed at the Katombora and Mambova rapids are believed to represent faulted blocks of the Chobe and Linyanti horst, respectively, and mark the points behind which a large paleolake was impounded. As a result of the drainage evolution of the Okavango-Kwando-Zambezi river system, a characteristic distribution of older and recent alluvial and lacustrine sediments has developed, stretching in the NE-SW direction across the ORZ (Moore and Larkin 2001; Thomas and Shaw 1991) and often confined by faults and dune fields (McCarthy 2013). Drainage evolution models developed by Cotterill (2006) and Moore et al. (2013) suggest that the age of large paleolakes in the area exceeds 500 ka. The largest lake was named Paleolake Deception (PLD) after a small eastward-flowing ephemeral tributary to the Makgadikgadi Basin and was discovered by McFarlane and Eckardt (2008). It presumably filled the basin to an elevation of 995 m a.s.l. and possibly covered a total area of 173,000 km 2 , about 2.5 times the size of Lake Victoria (Fig. 1a). Paleolake Makgadikgadi (PLM) represents the 945-m a.s.l. stage; its extent is clearly visible in aerial photos and satellite images by the Gidikwe Ridge (Fig. 1a), which is interpreted as an offshore barrier bar along the western shoreline of the paleolake. PLM that may have reached a total area of 65,000 km 2 has been extensively described in the literature over the last 50 years (Burrough et al. 2009;Grove 1969;Heine 1987;Meier 2008;Ringrose et al. 2005;Shaw 1988; Thomas and Shaw 1991). Interpretation of helicopter time-domain electromagnetic (HTEM) data provides evidence that PLM extended northeastward into the region presently occupied by the Okavango Delta (Podgorski et al. 2013). A smaller paleolake is associated with an elevation of 930 m a.s.l. and referred to as Paleolake Thamalakane (PLT). Shaw and Thomas (1988) proved the existence of a paleolake in the ZR, which they named Paleolake Caprivi (PLC). The PLC probably corresponds to the PLT stage. This lake was originally believed to have existed until fairly recent times, but the results of the radiocarbon dating of shoreline Mollusca which suggested at least two lake stages at 15 and 2-3 ka BP are now considered incorrect (e.g. Moore 2004) due to the fact that carbon samples are often contaminated by modern or old (Bdead^) inorganic carbon.
The sudden change in the flow direction of the Kwando River along the border between Namibia and Botswana is obviously related to tectonic uplift along the Linyanti Fault (Fig. 2). Thus, the change in flow from a southern to eastern direction must be understood as a tectonically driven river diversion rather than capture by headward erosion (Moore 2004). Sand spits at the eastern shoreline of the Sowa (formerly Sua) Pan (eastern Makgadikgadi Pan, Fig. 2) were interpreted as relict fluvial channels of the former Kwando River before its flow was diverted (Moore and Larkin 2001). This indicates that the flow of the Kwando was directed to the south into the Makgadikgadi depression prior to tectonism related to the propagation of the EARS.

Paleoclimate
It is generally recognized that, despite some notable interruptions, the western, northern (Sahara-Sahel) and eastern regions in Africa experienced overall warmer conditions and aboveaverage rainfall during the Early to Mid-Holocene, which corresponds to a warm interglacial. The interval is known as the African Humid Period (AHP), lasting-depending on the author-from 14.8-11.0 ka to 5.5-4.0 ka calibrated years before the present (cal. BP) (e.g. deMenocal et al. 2000;Gasse 2000;Willis et al. 2013 and ref. therein). The warmer and wetter conditions have been linked to a northward migration of the ITCZ due to an increased, partially orbitally driven interhemispheric insolation contrast and an intensification of the African monsoon. Cooler tropical temperatures and more arid conditions persisted throughout northern and most of East Africa during the Last Glacial Maximum (LGM). Various terrestrial records indicate broad arid phases between 23 and 17 ka, and lake levels in Africa were almost uniformly low during the LGM (Feakins and deMenocal 2010 and ref. therein).
The climatic conditions during the LGM and Holocene in the southern rainfall zone, and in particular the Kalahari environment, are more complex and are still a subject of controversy. This may be attributed to the extreme climate variability of the Kalahari Basin due to its position between a tropical summer rainfall and temperate winter rainfall regime, the scarcity of terrigenous proxies in arid climates and difficulties inherent in the various available dating technologies. There is generally good agreement that colder temperatures prevailed in southern Africa during the LGM, and that the winter rainfall zone expanded northward and experienced wetter conditions (Chase and Meadows 2007). Noble gas temperature records from the Stampriet Artesian Basin south of Windhoek (24°S, 19.5°E) and from a southeastern Kalahari aquifer near Letlhakeng, Botswana (24.1°S,25°E) suggest that temperatures during the LGM were lower than today by 5-6°C (Kulongoski et al. 2004;Stute and Talma 1997). According to Stute and Talma (1997), higher than expected δ 18 O values in groundwater provide evidence that the Stampriet area, which is today under the influence of rainfall originating largely from the Indian Ocean, may have received more Atlantic moisture during the LGM due to a northward shift of the winter rainfall zone. The dominant influence of the migration of tropical and subtropical circulation features on climate variability in southern Africa has been widely emphasized (Stone 2014). A conceptual climate model suggests that during interglacials, warming of the Northern Hemisphere (NH)-as the hemisphere with more land-is accentuated, which results in an overall northward shift of the ITCZ. For the African continent, this leads to increased rainfall over subtropical northern Africa (as during the AHP) and to a decrease in precipitation in the summer rainfall region (Chiang and Friedman 2012;Railsback 2017;Schneider et al. 2014). Conversely, due to the glaciation over the NH, and since the Southern Hemisphere (SH) has the larger fraction of ocean, cooling over the NH during glacial periods is stronger, which pushes the ITCZ southward into the SH. This results in increased aridity of the Sahara-Sahel zone (as during the LGM) and potentially higher rainfall in southern African regions, whereas intermediate conditions presumably prevailed during cooler interglacial periods. This conceptual model of climate change thus underpins the observation that climate variation patterns in northern Africa are very different from those in the SH (e.g. Gasse 2000;McCarthy and Rubidge 2005). The best evidence of this conceptual model predicting insolationdriven increased rainfall over the southern interior during glacial periods is provided by recent studies on sediment cores off the mouth of the Zambezi River (Schefuss et al. 2011;van der Lubbe et al. 2014van der Lubbe et al. , 2016, in contrast to comprehensive but partially contradictory results from thermoluminescence (TL) or optically stimulated luminescence (OSL) dating techniques performed on dune sands (Ivester et al. 2010;Munyikwa et al. 2000;O'Connor and Thomas 1999;Stokes et al. 1997a, b;Thomas et al. 2000Thomas et al. , 2003. The linear dune fields are widely considered to represent the oldest morphological landforms in the northern Kalahari Basin (McCarthy 2013). Chronological reconstruction of geomorphological events suggest that the dunes are older than 2 Ma (Moore et al. 2013). This is particular obvious in the area southwest of the Okavango Delta where the Gumare Fault ( Fig. 2) truncates the east-west-trending dunes. To the east of the Gumare Fault, the dune field has been completely destroyed by flooding, which is a clear indication that the formation of the dunes predates faulting as well as flooding by the development of the mega-paleolakes in the Makgadikgadi subbasin (McFarlane and Eckardt 2007). The determined younger ages of dune sands may be partially explained by dune degradation, bioturbation and overturning, as highlighted by McFarlane et al. (2005); it is plausible that older sands were temporarily exposed to sunlight and subsequently buried again.
Controversy also exists about the development and age of large paleolakes in the Makgadikgadi Basin which were associated with wet phases. A number of studies from the Lake Ngami, Mababe and Makgadikgadi basins that applied OSL and radiocarbon dating methods provide evidence that at least seven paleolake full phases may have occurred in the Makgadikgadi Basin during the late Pleistocene and Holocene (Burrough and Thomas 2008;Burrough et al. 2007Burrough et al. , 2009Gamrod 2009;Ringrose et al. 2005;Shaw et al. 1997). Moore et al. (2013), however, using water budget estimates, showed that the paleolake in the Makgadikgadi Basin could only reach the level of Gidikwe ridge (i.e. the PLM stage) if the Okavango, Upper Zambezi and Kwando rivers all terminated in the Makgadikgadi Basin. Models suggesting that paleolake levels were controlled by major tectonic events rather than climate variations are therefore more convincing and plausible. Moore et al. (2013) proposed that large paleolakes only existed >500 ka years ago in the basin and that they stepwise desiccated after major tributaries of the Upper Zambezi paleo-river system were cut off by tectonic subsidence and lifting and resulting processes such as headward erosion and river capturing.

Remote sensing
A digital elevation model was applied to identify characteristic landforms including dune fields, river fans and terraces, and faults and stream lineaments in the project area. The model uses high-resolution topographic global elevation data generated from NASA's Shuttle Radar Topography Mission (SRTM). The raster data were downloaded from the U.S. Geological Survey's EarthExplorer website (USGS 2015) on 22 January 2015, and have a resolution of 1 arcsecond, or about 30 m × 30 m. In addition, a true-color image dated 8 August 2017, delivered by the European Space Agency's Sentinel-2 satellite mission, was also obtained from the USGS EarthExplorer website (USGS 2017). The RGB (red, green, blue) composite image was created from bands 4, 3, and 2 and has a resolution of 10 m × 10 m.

Geophysical survey
During the years 2003/2004, about 57 transient electromagnetic (TEM) soundings were collected in the ZR. In 2004/2005, the survey was supplemented by an additional 17 soundings in order to extend survey lines toward the north and west (Fielitz et al. 2004). Poseidon Geophysics (Pty) Ltd., Gaborone, collected these soundings as part of the Namibian-German Technical cooperation Project BInvestigation of Groundwater Resources and Airborne Geophysical Investigation of Selected Mineral Targets in Namibia^.
The TEM soundings were collected using a GDP-32 receiver and a TEM/3 receiver coil (single-channel magnetic antenna) with an effective coil area of 10,000 m 2 . In-loop configuration was applied. The transmitter coil was 200 m × 200 m, single-turn, and a GGT-10 transmitter was used. All equipment was supplied by Zonge International (Zonge 2018). In order to achieve investigation depths of some 200 m, the measurements were performed applying 20 A for 1 Hz, 15 A for 4 Hz, and 5 A for 16 Hz-the lower the frequency, the longer the transient and the deeper the investigation depth. The measurements were stacked and averaged, 64 cycles for 1 Hz, 256 cycles for 4 Hz, and 512 cycles for 16 Hz. Standard processing software by Zonge (SHRED, TEMAVG) was used for stacking, averaging and visualization of the data. For data inversion, the TEMIX Plus software program (Interpex Limited, Interpex 2016) was used, and a 1-D inversion was applied.
The soundings were positioned along 16 profiles, partly crossing each other, and with irregular spacing between the soundings. The location of the soundings was aimed at covering the supposed area of the deep freshwater aquifer including its borders. The data quality was very satisfactory. In most cases, it was possible to collect transients up to 20 ms and longer, with low noise. At some locations, the transients collected during the first and the supplementary campaign differed. The detected resistivity shift could not be attributed to technical causes and was considered a reflection of the geological heterogeneity close to the northern graben structure, where most of the supplementary measurements were located.

Investigations of groundwater hydraulics
Groundwater contours for the UKA are based on water table data obtained from the database of the Namibian Department of Water Affairs and Forestry (DWAF). The water table measurements were converted into piezometric heads using elevations obtained from SRTM imagery (1 arcsec or approx. 30 m × 30 m, see Chapter 2.2). The dataset comprised about 750 values for water tables that were taken between 1970 and 2015 during different seasons, i.e. rainy or dry season. Variogram analysis and universal kriging with linear drift was applied for regionalization of water level data and mapping of groundwater contours. In order to construct groundwater contours for the LKA, 16 water table measurements were taken during the end of the dry season from 10 to 23 November 2015. Another nine older water table measurements were added from the database. Contours had to be constructed manually, as kriging or another interpolation method proved unsuccessful due to the limited number and quality of data.
In addition, results from pumping test analyses were statistically examined. Previous pumping test results were reported by Margane and Bäumle (2004) and Margane and Beukes (2007). Additional test data from the DWAF archives were analyzed. Finally, the initial, pumped and residual water levels were measured manually and continuously recorded with pressure probes during the sampling campaign conducted in November 2015. The results allowed for rudimentary test pumping analysis. In total, 42 single-well pumping tests from 35 boreholes were able to be analyzed: 20 for the UKA and 22 for the LKA. Pumping duration ranged from 1 to 72 h, and pumping rates also varied widely, from 0.2 to 58 l/s.

Sampling and analytical procedures
Groundwater samples from 22 boreholes were obtained during the field campaign carried out in November 2015. The bulk of groundwater samples were taken along a north-south-directed transect at around 23°45′E (Fig. 3). The boreholes are mostly used for providing drinking water or water for stock in cattle camps, and are installed with a submersible pump or in lesser cases with a hand pump. Five of the sampled boreholes are used as monitoring wells by DWAF and are thus permanently equipped with data loggers. Except for the monitoring boreholes, completion reports containing a detailed lithological description were not always available. It was hence difficult to determine with certainty the position of the confining unit separating the LKA from the UKA. It was deduced from this that boreholes with depths exceeding 130 m could potentially tap the LKA. Fifteen of the 21 sampled boreholes fulfilled this depth criteria (maximum depth 250 m). It should be noted that for boreholes 037223, 200501, 200503, 200504 and 200620, the position of the screens is unknown, and that for boreholes 202147, 202148 and possibly 201624, the upper parts of the screens are located within the UKA (Appendix Table 2). During this and two previous campaigns, a number of additional samples were taken from boreholes of the UKA and the Kwando River, as well as from a PVC rainfall collector that was set up at Camp Kwando Lodge (23°19′E, 18°03′S, 963 m a.s.l.), about 25 km south of Kongola Bridge, for purposes of comparison.
Physicochemical parameters were measured using commercial portable meters (WTW/Xylem Group) as part of through-flow cells to avoid any contact of groundwater with atmospheric oxygen and carbon dioxide. In situ measurements included temperature (T), electrical conductivity (EC), pH, dissolved oxygen (DO) and oxidation-reduction potential (Eh)(measured against an Ag/AgCl electrode and corrected with respect to standard hydrogen electrode). The pump rate was measured volumetrically using a 20-l container and a stopwatch. Alkalinity and acidity were determined onsite by potentiometric titration on 100-ml unfiltered samples using a field burette. The normality of the titrant solutions (HCl, NaOH) was 0.047 N. The titration was considered important, as in previous campaigns it was observed that pH measured in situ was usually higher (more alkaline) than in samples in the laboratory. The lowering of the pH was explained by dissolution of atmospheric CO 2 in the sample water. A fixed-endpoint method (to pH 4.3), and USGS Inflection Point Titration (IPT) and GRAN function plot methods were applied for analysis (Rounds 2012). The data were entered in a customized spreadsheet and preliminarily analyzed on a laptop in the field. The samples were taken once all in situ parameters except Eh had completely stabilized, which took between 60 and 120 min. Eh values dropped continuously from the start of pumping, but toward the end of pumping, the fluctuations gradually diminished. The submersible pumps yielded pump rates in the range of 0.8 to 1.4 l/s, which resulted in a total abstracted volume of water of between 2 and 8.6 m 3 . Much smaller amounts, however, were extracted from boreholes with hand pumps. The samples were continuously kept at low temperature using cooling boxes and frozen cooling elements.
Inorganic components and the stable hydrogen and oxygen isotopes were analyzed by the BGR laboratory in Hannover, Germany. Carbon isotopes and chlorine-36 were analyzed by the external accredited laboratory, Hydroisotop GmbH, in Schweitenkirchen, Germany. The applied analysis methods are summarized in Table 1. The isotopic composition of hydrogen (δ 2 H) and oxygen (δ 18 O) was determined using a Picarro cavity ring-down spectrometer (CRDS) equipped with a CTC autosampler. The measurements were normalized to the Vienna Standard Mean Ocean Water/Standard Light Atmospheric Precipitation (VSMOW/SLAP) scale, where values of 0 and −428 ‰ for δ 2 H and 0 and −55.5 ‰ for δ 18 O were assigned to VSMOW and SLAP, respectively. The external reproducibility (standard deviation of a control standard during all runs) was better than ±1.0 ‰ for δ 2 H and ±0.1 ‰ for δ 18 O. The results for δ 18 O and δ 2 H are expressed in the conventional notation as relative deviations in per mil (‰ or 10 −3 ) from the VSMOW according to the following relationship (Coplen 2011): where R sample is the isotopic ratio of the sample ( 2 H/ 1 H or 18 O/ 16 O) and R VSMOW is the isotopic ratio of the standard. The deuterium excess D ex is calculated from the isotopic composition of δ 2 H and δ 18 O as follows: The D ex is about 10 ‰ for rain that plots near the Global Meteoric Water Line (GMWL), but is much lower for surface water and groundwater that have lost some volume to evaporation. The D ex thus serves as a valuable tool for tracing the origin of water.
Carbon-13 ( 13 C) content was determined using isotope ratio mass spectrometry (IRMS). The 13 C abundance is reported as the deviation in per mil (‰ or 10 −3 ) of the ratio 13 C/ 12 C in the sample relative to the 13 C/ 12 C ratio of the Vienna Pee Dee Belemnite (VPDB) standard, i.e. as δ 13 C ‰ VPDB. Measurement reproducibility of duplicates was better than ±0.3 ‰. Carbon-14 ( 14 C) was analyzed using accelerator mass spectrometry (AMS), and 14 C activity is given in percent modern carbon (%mC). 14 C content expressed in %mC represents the deviation from 14 C activity in the undisturbed atmospheric equilibrium that equals 100 %mC or 0.226 Bq/g of carbon. Chlorine-36 ( 36 Cl) is measured as 36 Cl/Cl ratios and is also determined by AMS. Chloride needs to be extracted from acidified solutions by adding excess AgNO 3 , resulting in the precipitation of AgCl. Samples with chloride content of less than about 20 mg/L are considered critical, as the mass of precipitated chloride may be insufficient for determining the radioactive isotope 36 Cl. Since sulfur concentrations are generally high in the study area, and because the sulfur-36 isotope interferes with 36 Cl during AMS analysis, sulfur has to be removed from the samples in an intermediate step. This is achieved by precipitation as BaSO 4 (with Ba(NO 3 ) 2 ) and subsequent removal by centrifugation. The 36 Cl concentrations are expressed as 10 −15 atoms of 36 Cl per total atoms of Cl.

Groundwater age dating
Radioactive tracers can be applied to determine the apparent age of groundwater. Radioactive tracers used for the longer timescales are 14 C and 36 Cl. Both are cosmogenic nuclides, naturally produced in the upper atmosphere from cosmic radiation. For radioactive tracers with known or constant input function, an apparent age can be derived from concentrations using the equation of radioactive decay. Radiocarbon has a half-life of 5730 years, which allows the determination of timescales between some thousand and approximately 35,000 years. 36 Cl has a half-life of 301,000 years and is suitable for determining ages on the order of 10 5 to 10 6 years (e.g. Kalin 1999;Phillips 1999;Mook 2000). Radioactive decay is described by the exponential function: where λ t is the decay constant of the tracer (λ 14 C ¼ 1:210•10 −4 a −1 ; λ 36 Cl ¼ 2:303•10 −6 a −1 ) C is the measured tracer concentration of the sample C re is the (initial) tracer concentration of recharge t is the time corresponding to the apparent groundwater age The apparent groundwater age of a sample is calculated as Phillips (2013) described an easy mass balance equation for the absolute 36 Cl concentration in an aquifer. The model incorporates radioactive decay of the recharged 36 Cl, underground production in the aquifer toward a secular equilibrium ratio of 36 Cl /Cl and an additional source of chloride, such as diffusion from an adjacent aquitard, which may have another 36 Cl /Cl ratio in secular equilibrium. The Phillips model expressed for the tracer concentration A 36 Cl as given in Suckow et al. (2016) reads as follows: If R se2 = R se1 , the equation simplifies and reads: and the apparent groundwater age is calculated as (Bentley et al. 1986;Phillips et al. 1986) Furthermore, if chloride is added from dissolution of halite, then R se = 0, and consequently (Phillips 2013): In this case, the equation for apparent groundwater age reads: where Cl and Cl re represent the measured chloride concentration and the chloride concentration of recharge in mg/L, respectively.

Remote sensing
The high-resolution DEM and true-color image proved suitable for detecting landforms and tectonic regimes in the study area. The Linyanti and Chobe faults are easily detectable in the DEM. The Linyanti fault line located south of the Linyanti River and the Magweqgana Stream is expressed in vertical displacements on the surface of about 4-8 m. Paleolakes likely existed several times in the down-faulted areas north of the Linyanti and Chobe faults. Lake Liambezi (Fig. 2), which episodically fills up after extreme flooding, can be regarded as a remnant of these paleolakes. The Sibbinda Fault to the north marks the boundary of a zone of gradually rising elevation and is considered a blind fault. The elevation south of the fault line is about 970 m above mean sea level (a.s.l.), whereas the plateau north of the fault reaches values well above 1000 m a.s.l. (Fig. 2).
The Kwando and Zambezi rivers have deposited fans into the ORZ similar to but smaller than the Okavango Delta to the southwest (Fig. 4). These fans were also recently discovered from digital elevation data by Burke and Wilkinson (2016). The Kwando Fan and parts of the Zambezi Fan are essentially inactive, because its formative rivers have incised into the relict fan deposits. The fans therefore form the upper part of the older alluvium. Figure 4 also illustrates the strong tectonic control on the drainage system in the study area. Apart from the obvious elbow-shaped river diversion near the Linyanti Fault, the course of the Kwando River is also deflected south of the Matabele-Mulonga Plains and north and south of Kongola. The Matabele-Mulonga Plains, in fact, are believed to represent an abandoned arm of the Kwando River, which hence initially joined the Zambezi in the area of the Barotse-Bulozi Floodplains (Moore et al. 2013). The near-parallel course of the Kwando and Upper Zambezi rivers, as well as the course of the eastern tributaries to the Zambezi such as the Lumbe and Njoko rivers, is also most likely an expression of the regional tectonic setting. The Ngonye Falls near Sioma provide evidence that the section of the Fig. 4 Prominent landforms in the border area between Angola, Namibia, Botswana and Zambia as seen from a Sentinel-2 RGB composite image mosaic taken on 8 August 2017; abbreviations of faults as in Fig. 2 Zambezi River downstream of the falls is a relatively young valley. The rapids are a clear sign of rejuvenation of an otherwise very mature longitudinal river profile. The Zambezi obviously once followed a route to the west of its modern course crossing the Siloana Floodplain, as indicated by the arrows in Fig. 4. Dune fields in the Northern Kalahari generally developed to the west of large river valleys which supplied the sands for dune formation. To the west of the Kwando and Okavango rivers, linear dune fields are well preserved. The fact that the Siloana Floodplainin contrast-is characterized by largely degraded dunes and areas without remnant dunes is another indication that the area has undergone extensive fluvial reworking.

Geophysics
The supplementary TEM measurements revealed the boundaries of the alluvial aquifer system toward the east and north (Fig. 5). The basalt appears at sounding 4_16 at profile 9, at sounding 1_200 at profile 1 and at sounding 3_50 at profile 3. All soundings north/west and north along these profiles show the indication of high-resistive basalt (Fig. 5c). The TEM inversion results show a decrease in resistivity toward the alluvial basin. Along profile 9, the basalt dips toward the southeast, from about 20 m depth at sounding 4_08 to about 100 m depth at sounding 4_16. The alluvial filling, with its relatively low resistivity of some 20 Ωm, is clearly seen from sounding 4_20 toward the southeast. From sounding 4_08 to sounding 6_5, the resistivity of the basalt layer drops from about 1000 Ωm at 4_08 to some 50 Ωm at 6_5.
At the southern end of the TEM profiles, no indication of higher resistivity at depth is detected. In Fig. 5c, the inversion result of the most southern sounding of the whole campaign (sounding 1_30) shows very low resistivity at depth (about 5 Ωm). The shape of the transient gives no hint of resistive layers at depth. However, owing to the low resistivity of the layers at this site (below about 120 m depth), the investigation depth of this sounding is limited. Although the figure displays the vertical profile to a depth of some 350 m, such investigation depth cannot be reached by the measurement layout as described above given the low resistivity of the alluvial layers. From theoretical calculation it can be concluded that a resistive layer deeper than about 280 m would remain unnoticed by the TEM sounding at this point. Figure 6 shows the piezometric surfaces of groundwater in the UKA and LKA. For the UKA, only minor manual corrections of automatically generated contour lines-mainly in areas with low measurement density-were required. Groundwater of the UKA flows in a north-south direction in the eastern parts of the ZR and in a west-east direction near the Kwando River floodplains (Fig. 6a). The latter is an indication that losing-stream conditions occur in the floodplains of the Kwando River.

Groundwater flow
The LKA appears to be confined throughout the entire area. Based on water level measurements taken in November 2015, the piezometric head in the LKA near borehole no. 041004 (Fig. 3) is over 14 m higher than in adjacent boreholes of the UKA. To the south, the measured head difference at the borehole group at site no. 200254 for the same period amounted to about 4 m. The piezometric head difference hence remains positive but gradually decreases toward the south. However, reliable groundwater contours could not be established until now. This is due to the fairly small number of boreholes, and hence limited data, and the prevailing low hydraulic gradients. In addition, most boreholes penetrate the aquifer only partially, and individual screen depths vary considerably. It is also probable that some boreholes are connected to the UKA. Despite these limitations, a general north to south flow direction could be established from hydraulic head measurements taken during November 2015 (Fig. 6a). The hydraulic gradient in this direction is estimated at 0.08 ‰.
The variable sedimentological environment, with phases of fluvial, fluvio-deltaic and lacustrine conditions, is not only responsible for the high lithological variability, but also explains the variations in hydraulic transmissivity. The box-plots shown in Fig. 6c represent the results from the available 42 single-well pumping tests. The LKA is characterized by higher transmissivity with a lower dispersion compared to the UKA. The median values of transmissivity amount to 59 m 2 /d for the LKA compared to 20 m 2 /d for the UKA.

Hydrochemistry
The general chemical groundwater composition exhibits strong heterogeneity, as shown in the percentile plots of Fig. 7. Due to the effects of evapoconcentration and dissolution of salts, groundwater often becomes brackish. Total dissolved solid (TDS) content of up to 600 mg/L can largely be explained by the dissolution of calcite and magnesite. Weathering of silicates also plays a role; silica content, in fact, is exceptionally high, sometimes exceeding 100 mg/L SiO 2 , due to the relatively high solubility of quartz at the prevailing groundwater temperatures of >25°C and as a result of very active CO 2 production in subtropical soils. TDS values of above 600 mg/L indicate a growing influence of evapoconcentration and/or dissolution of halite and gypsum, leading to the dominance of Na, Cl and SO 4 content in groundwater (Plata Bedmar 1999). Groundwater within or near the Kwando floodplains is generally fresh, with low TDS (median of 444 mg/L) and very low sulfate and chloride (not included in Fig. 7) content. Groundwater of this section is usually of the calcium-magnesium-bicarbonate type. Other areas show overall increased levels of sodium, chloride, sulfate and total mineral content. Boreholes in the older alluvium and in the Linyanti-Chobe swamps show extreme variations; brackish groundwater with high sulfate content is common for both areas. By contrast, the LKA shows a homogeneous composition with fair water quality (Fig. 7). Median values for TDS, sulfate and chloride are 959, 211 and 126 mg/L, respectively.
Chemical reactions such as dissolution, precipitation, ion exchange and redox reactions control and may alter the chemical composition of groundwater along the flow path. The presumed flow direction within the LKA is from the Namibian border with Zambia in the north to the border with Botswana in the south. The total flow length is approximately 80 km. The general evolution of groundwater chemistry leads from a Ca-Mg-HCO 3 to an Na-(HCO 3 + Cl + SO 4 ) type, as shown in the piper diagram (Fig. 8). Along the groundwater flow direction, a notable increase in TDS, chloride, sulfate and sodium content, as well as alkalinity, can be observed (Fig. 9). In contrast, calcium and magnesium content are overall very low, with values between 2 and 4 mg/L. Selected tabulated chemical analysis results obtained during this study are given in Appendix Table 2. Saturation indices calculated by the PHREEQC model (Parkhurst and Appelo 1999) show that the groundwater of the LKA is generally saturated with respect to calcite, magnesite and amorphous silica, and undersaturated with respect to all other major salts including gypsum, fluorspar and, obviously, halite. TDS values exceed 1500 mg/L in the southernmost areas near the Linyanti River. The correlation with flow distance is fairly strong, as shown by the Pearson product-moment correlation coefficient, r, of 0.88. Chloride concentrations increase from Fig. 6 a Piezometric surface and hydraulics of UKA and LKA in the Zambezi Region. Magenta and blue isolines denote hydraulic heads of the UKA and LKA, respectively. Broken lines show tentative contours in areas with insufficient data. b Variogram analysis prepared for construction of groundwater contours of UKA. Note that contours for LKA were drawn manually. c Box-whisker plots of transmissivity (T) from single-well pumping tests in the two aquifers below 100 mg/L to over 200 mg/L, with a median of 126 mg/L. The correlation is not as strong as for TDS (r = 0.62) due to a significant number of outliers. Note that borehole no. 037723, whose chemistry does not seem to fit into the regional picture, is located well outside the boundaries of the known extent of the LKA (Fig. 3). Sulfate content also increases steadily (r = 0.67) and reaches very high levels of above 300 mg/L. The SO 4 /Cl ratio (meq:meq) is on the order of 0.8-1.5, which is about 7.5-15 times that of the marine ratio of 0.104, and shows no clear trend in the N-S direction. Four rainfall samples collected at the beginning of the rainfall seasons in November 2014 and 2015, respectively, show high variations in the SO 4 /Cl ratio of between 0.5 and 5.0. Unfortunately, additional data, especially during the peak of the rainy season, is not available. Although deposition of ions in bulk precipitation (precipitation + dry fallout) in windy, semiarid environments might be higher than in humid and coastal areas, evapoconcentration of rainwater cannot explain the increase in chloride and sulfate along the flow path under confined conditions. Hence, the increase in chloride and sulfate is most likely predominantly caused by dissolution of halite and gypsum. The Br/Cl ratio can be indicative of halite dissolution. Whilst moderate evapoconcentration causes only minor changes in the Cl/ Br ratio, halite dissolution may increase the value of molar Cl/Br to several thousands (Braitsch and Günter 1962). Molar Cl/Br in marine water is about 655 ± 4, and that of continental rainwater between 50 and 650 (Alcalá and Custodio 2004;Fontes and Matray 1993). Observed values in the LKA show no trend in the N-S direction and vary between 950 and 1400. They are hence indicative of moderate halite dissolution. Higher and critical (>1.5 mg/L) fluoride concentrations occur in the south of the study area. Dissolution of fluorspar facilitated by low calcium content could be one explanation; however, this would not explain the increase in fluoride toward the south, as groundwater is undersaturated with respect to fluorite throughout the area. It is therefore more likely that desorption from clay minerals-facilitated by high pH (>8.5) and competition with hydroxyl ions-is the major mechanism, as described by Jacks (2016). Sodium experiences a strong and almost uniform increase (r = 0.81) along the flow path, with sodium content increasing from about 50 to >600 mg/L. The Na/Cl ratio (meq:meq) remains fairly stable, with only a minor increase toward the south. The values, typically in the range of 3-4, show that halite is not the only source of sodium in the system. The increase in total dissolved inorganic carbon (TDIC) and, therefore, total alkalinity is also strongly correlated with flow distance (r = 0.89). Maximum values exceed 100 mg C/L. Alkalinity clearly increases by a factor of 6 to 8 at a rate of about 1.54 mg C/ L per km flow distance. Furthermore, a striking correlation exists between the increase in alkalinity and in Bsodium-excess^, Na-Cl, defined as the molar difference between sodium and chloride in a sample. Na-Cl represents sodium that originates from sources other than halite dissolution, such as silicate weathering and ion exchange processes. The ratio of Na-Cl vs. (alkalinity + SO 4 ) within the LKA is almost consistently ±1, as shown in Fig. 9.
The processes that can best explain the geochemical evolution of groundwater in the LKA are cation exchange processes that consume alkaline earth metals and produce Na. The continuous removal of Ca by cation exchange along the flow path drives the progressive dissolution of calcite along the flow path and contributes alkalinity to the groundwater. In addition, sulfate is added by dissolution of gypsum. The geochemical evolution can be described by the chemical reactions: 2 Ca 2þ þ NaX→ 1 2 CaX þ Na þ ð11Þ and results in an Na-HCO 3 -SO 4 water type. Due to the observed low redox potential, and rotten smell at some sites, there is a possibility that alkalinity is also produced by reduction of sulfate following the equation However, considering the high SO 4 content in the LKA, this process is most likely of minor importance. The increase in alkalinity and the development of Na-HCO 3 water is a common phenomenon in deep-seated old groundwater, and was also reported in the Great Artesian Basin in Australia (Cresswell and Duesterberg 2012; Suckow et al. 2016), the Cretaceous sandstones of the Milk River Aquifer in Alberta, Canada (Phillips et al. 1995), and the Triassic-Cretaceous Guarani aquifer system in Brazil (Gastmans et al. 2010). Calcite is known to be abundant in the Kalahari sediments in the form of calcareous sandstone and calcrete. The prevalence of sodium as the main sorbate and the occurrence of gypsum are indicative of a more saline environment. In fact, the processes resemble that of freshening of a salty aquifer in coastal areas as described, for instance, by Appelo and Postma (2007) and van der Kemp et al. (2000). Freshening is also known to occur in aquifers within continental sedimentary basins such as the transboundary Guarani Aquifer system in South America (Houben et al. 2015;Sracek and Hirata 2002), which is composed of continental fluvial and aeolian sandstones. The Guarani Aquifer was deposited under arid conditions similar to the Northern Kalahari environment. Consequently, the system was originally brackish. It is assumed, however, that in the Guarani Aquifer, large amounts of salts may have already been flushed out after Cretaceous tectonic uplift. In contrast to the LKA, high sodium, chloride and sulfate concentrations in the center of the Guarani Basin possibly originate from mixing with water from the underlying, more saline unit or diffusion from incompletely flushed zones of lower permeability within the aquifer. What is most noteworthy about the downgradient chemical evolution observed in the LKA is that it confirms hydraulic continuity along a north-southdirected flow direction despite very low and uncertain hydraulic gradients, and provides evidence that the LKA is an Fig. 8 Piper diagram showing the general evolution of groundwater chemistry from a Ca-Mg-HCO 3 to an Na-(HCO 3 + Cl + SO 4 ) type; arrows show observed trends along the flow direction Fig. 9 Evolution of ion content and ratios in the LKA along the N-S transect. Unless otherwise stated, concentrations are represented by black hollow triangles, and ion ratios by red crosses; red lines denote trends in ion content active hydrogeological system that still forms part of the modern hydrological cycle.

Hydro-isotopes and apparent ages
Hydro-isotope analysis results are given in Appendix Table 3. Figure 10a shows the stable isotope composition of groundwater, surface water and rainfall samples. The bulk of the samples from the LKA show the lowest observed values of deuterium and oxygen-18 of between −56 and −62 ‰ and −8.7 and −8.1 ‰ VSMOW, respectively. These samples plot near the Local Meteoric Water Line (LMWL) determined by Beyer et al. (2016) for the Ohangwena Region in northern Namibia. Four exceptions occur, namely the samples from wells no. 201624, 202147, 202148 and 041002. Their position in the diagram is clearly below the LMWL. These four wells are located in the northern margins of the known extent of the LKA (Fig. 3). Of these wells, the LKA has been properly sealed off only in well no. 041002; the screen position at boreholes no. 201624, 202147 and 202148 is 120-132 m, 90-145 m and 85-129 m, respectively (Appendix Table 2), and suggests that groundwater pumped from these wells represents mixed water from the UKA and LKA, which can explain the uncharacteristic isotope signature. Together with wells of the UKA, these samples lie near an evaporation line (EL) and are isotopically enriched due to kinetic isotopic fractionation during non-equilibrium evaporation of water. This is in good agreement with stable isotope results from previous studies in the ZR, which are shown together with the new analyses in Fig. 10b. The EL determined from all available samples follows the relationship y = 5.47x − 12.4 (r = 0.99). Samples plotting along the EL are generally characterized by low D ex (Eq. (2)) of between 5 and −10 ‰, whereas most samples of the LKA have values between 7.7 and 10.3 ‰; the LKA values are closer to the D ex of the local and global meteoric water, which implies that the LKA is replenished by direct recharge and fast infiltration (e.g. during storm events) and that indirect recharge from rivers or open water areas is unlikely. The orange line in Fig. 10b represents a mixing line and indicates that mixing of different groundwater with various exposure to evaporation could also, to some extent, explain the stable isotope variations. Mixing of groundwater of different age and origin may be induced by sampling boreholes with long well-screen sections which are typical for semi-arid areas like the ZR. The number of surface water samples (Fig. 10a) is insufficient for a reliable interpretation. As is to be expected, the available data also plot on an evaporation line indicative of evaporative losses; the best fit to the surface water samples shown in Fig. 10a, however, appears to differ from the EL for groundwater.
The isotope content of modern precipitation differs significantly from the characteristic low isotope content obtained for the LKA. The weighted average isotope content at climate stations of the Global Network of Isotopes (GNIP) in southern Africa is added to Fig. 10a (GNIP 2017). The Online Isotopes in Precipitation Calculator (OPIC) predicts-based on interpolation of the isotopic composition of modern meteoric precipitation-values of deuterium and oxygen-18 of −28 ± 6 ‰ and −4.6 ± 0.8‰ VSMOW for the study area (Bowen 2017;Bowen and Revenaugh 2003). The much lower values for the majority of sites from the LKA (and a noteworthy number from the UKA) could be explained by the temperature effect on isotope fractionation, a change in seasonality (i.e. from summer-to winter-dominated) or a stronger influence Fig. 10 a Stable isotope binary plots for groundwater of the UKA and LKA, surface water and precipitation. Hollow symbols represent weighted average isotope content of precipitation at selected stations in southern Africa (GNIP 2017). Local Meteoric Water Line (LMWL) after Beyer et al. (2016). b Comparison of groundwater stable isotope composition of this study with previous results; *data from Margane et al. (2005); **data from Plata Bedmar (1999). Numbers in orange give mixing ratios between assumed end-members of continentality that could be produced by a change in the circulation pattern. As there is no evidence for a fundamental change in the rain-bringing climate systems in the region, the variation in average temperature is the most likely explanation. The lower isotope concentrations hence indicate that groundwater was formed under cooler than modern climatic conditions, such as during the Pleistocene, when average annual temperatures in southern Africa were about 5°C lower than today (Holmgren et al. 2003;Kulongoski et al. 2004). The amount effect, i.e. if the percentage of high-intensity rainfall events is increased, could also contribute to lowering the δ 18 O and δ 2 H values of rainfall.
The bulk of the LKA samples are characterized by remarkably high 13 C content of between −6.6 and −4.4 ‰ VPDB and very low 14 C activity below 6 %mC (Fig. 11a). Exceptions are again boreholes no. 201624, 202147 and 202148, which are assumed to partially tap the UKA, with δ 13 C < −11 ‰ VPDB and 14 C concentrations between 30 and 42 %mC. Borehole no. 037723, which is 164 m deep but located to the northeast of the proven extent of the LKA (Fig. 3), also shows both lower 13 C and higher 14 C concentrations. Borehole no. 041002 has 13 C content typical of the LKA but much higher 14 C activity. In an open system, the carbon-13 content of dissolved inorganic carbon in groundwater, δ 13 C DIC , depends on the 13 C content of soil CO 2 , δ 13 C CO 2 g ð Þ , which is generated by the degradation of biomass, as well as pH and temperature. The pH determines the relative abundance of carbonate species CO 2(aq) , H 2 CO 3, HCO − 3 and CO 2− 3 in the water, and temperature influences the fractionation enrichment factors (e.g. Clark 2015). Common trees and bushes (Acacia, Combretum, Sclerocarya) near the Tswaing Crater (25°25′S, 28°5′E) show organic δ 13 C org of −30.7 to −25.5 ‰, whereas grasses and macrophytes (e.g. Cyperaceae) show values of −13.2 to −11.4 ‰ VPDB (Kristen et al. 2010). The lower 13 C content of many grasses is due to the fact that they use C4 photosynthesis, whereas most savannah trees use C3 photosynthesis. C4 plants are able to minimize water loss in hot, dry environments due to improved photosynthetic rate and efficiency, and are therefore better adapted to drier climate and lower atmospheric CO 2 concentrations (e.g. Lara and Andreo 2011;Osborne and Sack 2012). Assuming that plants using C3 photosynthesis dominate, then δ 13 C CO 2 g ð Þ typically ranges from −23 to −20 ‰ VPDB (Clark 2015). At a temperature of 28°C, the expected δ 13 C DIC of groundwater in an open system is between −16.9 and −13.9 ‰ at pH 7, and −15.5 and −12.5 ‰ VPDB at pH 8. For a closed system, exchange with soil CO 2 is inhibited, and carbonate dissolution will reduce the initial 13 C content by half according to the stoichiometry of the carbonate dissolution. 13 C concentrations in calcretes, δ 13 C calcrete , worldwide range from −14 to 4 ‰, with a maximum of −6 to −4 ‰VPDB (Talma and Netterberg 1983). Using −24 to −21 ‰ for δ 13 C CO 2 aq ð Þ and −4 ‰ for solid freshwater carbonate, the resulting δ 13 C DIC of groundwater is −14 to −12.5 ‰ VPDB, and hence in a similar range as that of an open system. Samples from the UKA and the few outliers of the LKA samples have 13 C concentration similar to the values to be expected for open and closed systems. The majority of the samples from the LKA, however, have much higher values. One reason for this could be higher than presumed values of δ 13 C calcrete or a different CO 2 source. Unfortunately, sedimentary analyses of calcretes are not available in the study area. According to Talma and Netterberg (1983), soil δ 13 C CO 2 g ð Þ of −18 to −7 ‰ VPDB was measured in grassland-dominated areas of southern Africa. An abundance of C4 plants thus could be responsible for higher initial soil δ 13 C CO 2 g ð Þ . Assuming soil δ 13 C CO 2 g ð Þ to be about −12 to −10 ‰ during the time of recharge of groundwater, this results in δ 13 C DIC of groundwater of −8.5 to −7.5 ‰ for closed-system conditions. For the LKA, this could mean that grassland may have prevailed during recharge. Recharge under grass-covered areas is usually higher than under densely wooded grassland, even if mean rainfall amounts are reduced due to the absence of deep-rooting trees. This observation is particularly pronounced if storm events with high rainfall intensity prevail (e.g. Kim and Jackson 2012).
Another, more likely, explanation assumes that the groundwater of the LKA is very old and the δ 13 C DIC concentrations are dominated by the carbon isotope characteristics of rock. This is because the observed progressive dissolution of carbonate associated with an increase in alkalinity and DIC steadily dilutes the original δ 13 C DIC (and 14 C) content. Furthermore, due to the long residence time, solute-matrix exchange occurs and groundwater DIC will approach isotopic equilibrium with carbonate rocks (e.g. Han et al. 2014). Assuming pH > 8, temperature of 28°C and δ 13 C calcrete of −4 ‰, groundwater DIC at isotopic equilibrium with the rock matrix is enriched by about 1 ‰ compared to calcretes. The resulting value of −5 ‰ VPDB is in good agreement with measured δ 13 C DIC values. Similar effects have been observed in central parts of the Guarani Aquifer system, where dissolution of calcareous cement and isotopic exchange resulted in increased 13 C isotope ratios of about −6 to −8 ‰ (Gastmans et al. 2010).
Dissolution of carbonates is obvious from the observed geochemical evolution and increase in TDIC down the flow path. It is therefore tempting to assume that the very low observed 14 C concentrations are linked to the dissolution of dead carbon. 36 Cl/Cl ratios, however, are very low, suggesting that groundwater is much older than 50 ka and that 14 C must have virtually fully decayed due to its shorter half-life. This becomes obvious from the binary plot in Fig. 11b that shows the observed 14 C concentrations and 36 Cl/Cl ratios together with a theoretical isotope composition following simple exponential decay of the estimated respective initial concentrations. The plot can also be used to estimate the experimental detection limit for 14 C concentrations during this investigation, as suggested by Suckow et al. (2016). The purple horizontal line in Fig. 11b suggests that concentrations below about 5-6 %mC are affected by contamination of the samples by atmospheric CO 2 during sampling and/or analysis. This error range is plausible, as exposure of the samples to air could not be fully avoided, and CO 2 uptake is facilitated by the alkaline pH of the samples . Mixing between older groundwater, characterized by high 13 C content, low 36 Cl/Cl ratios and 14 C of zero, with younger water of the UKA could be responsible for the different isotope content of the group of outliers within the LKA. However, as the measured values do not follow a simple binary mixing line-as indicated by the orange lines in Fig. 11-a more complex scenario involving additional end-members of mixing is apparent.
The measured 36 Cl/Cl ratios R for samples of the LKA range from 7 × 10 −15 to 50 × 10 −15 (atoms:atoms), with a median of 20 × 10 −15 . The values for the three boreholes of the UKA and of the deep borehole 037723 located to the northeast of the LKA show significantly higher 36 Cl/Cl ratios (Fig. 11b). For three boreholes, namely 034611, 035060 and 200248, R could not be determined because the chloride concentrations of the samples were too low (i.e. <3 mg/L). The challenge in interpreting the 36 Cl data on the one hand is in estimating the initial 36 Cl/Cl ratios, R re , as well as the secular equilibrium 36 Cl/Cl ratios in the aquifer and of a potential additional source of Cl, R se . On the other hand, values of Cl vary considerably in the samples, from below 10 to 279 mg/L, which means that the origin of the Cl added after recharge needs to be clarified. The two highest measured values of R are 161 and 176 × 10 −15 , which set a lower constraint to R re . The uncorrected 14 C ages of these samples are 3651 and 9655 years, which are too low to have a significant impact on R re . In published studies applying the 36 Cl dating technique (still a rather modest number), estimates of R re are mostly between 100 × 10 −15 and 500 × 10 −15 . A value of 200 × 10 −15 thus appears suitable for this study. R se is estimated at 5 × 10 −15 ; this value is also similar to published data. Measured R values for the LKA are very low, and result in uncorrected apparent ages, i.e. assuming simple decay of the initial recharge concentration (Eq. (4)), of between 626 and 1433 ka. The addition of chloride could in principle be explained by variable evaporation rates in the vadose zone of the recharged water, or by addition of chloride along the flow path by dissolution of halite, mixing with brine/connate water or diffusion from adjacent aquitards. An evaporated recharge source yields a reasonable fit in a plot of R vs. the 36 Cl concentration, A, as shown in Fig. 11c. In this example, the average chloride concentration of LKA samples, amounting to about 130 mg/L, was chosen as input. The identified geochemical evolution down the flow path described above and the high molar Cl/ Br ratios between 950 and 1400 of groundwater from the LKA, however, suggest a different source. The average linear groundwater velocity determined by Darcy's law is 0.15 m/a applying estimates for hydraulic conductivity of 1 m/d, porosity of 0.2 and hydraulic gradient of 0.08 ‰. A value for the relative flux of chloride F of 5.9 × 10 9 atoms/L/a (Eq. (6)), corresponding to an addition of 3.5 × 10 −4 mg/L/a of Cl, leads to an increase from an initial chloride content of 20 mg/L to about 210 mg/L at a rate of 2.3 mg/L/km over the flow distance of 80 km (Fig. 12a). Two scenarios are investigated: one in which the additional Cl is added by halite dissolution (i.e. R se = ,0) and a second in which the chloride is added by mixing with connate water or diffusion from an aquitard (R se ≠ 0). In the first scenario, A 36 Cl is unaffected but R is Bdiluted^by the addition of dead chlorine (Eqs. (9) and (10)). In the second case, underground production, such as by diffusion from an adjacent aquitard, adds 36 Cl, and with increasing time, R in the aquifer will approach secular equilibrium R se (Eqs. (7) and (8)). In Fig 12, the piston-flow model (PFM) refers to simple radioactive decay of recharged 36 Cl in an advective transport model, ADH refers to the first scenario involving added dissolved halite, and ASC to the second scenario with an added source of 36 Cl by diffusion, etc. Figure 12 shows the variation in Cl, R, A 36 Cl and apparent groundwater age along the presumed N-S-directed flow direction from Zambia across the study area. For the given value of R se , the PFM only fits to measured A 36 Cl , but not R. The ADH and ASC models, in contrast, produce similar curves and appear to be suited for simulation of both R and A 36 Cl . The measured 36 Cl/Cl ratios do not decrease significantly south of about 40 km, and instead appear to approach a possible secular equilibrium (Fig. 12b). The ASC is somewhat favored over the ADH, as it is the only model in which R does not approach zero, and A 36 Cl can rise with increasing time and travel distance due to underground production of 36 Cl (Fig. 12c). Figure 12d shows the apparent age of the samples calculated  4) and the ASC model according to Eq. (8). Both show a general increase in groundwater age with travel distance. The ASC produces overall lower ages because aquifer production of 36 Cl is superimposed on radioactive decay of recharged water. Nevertheless, ages also substantially exceed 100 ka for this model. Due to the high uncertainty on both the hydraulic and isotope input parameters, however, absolute groundwater ages cannot be determined with sufficient accuracy. At this stage, the simulations therefore are aimed at providing a better understanding of the flow and transport processes, confirming hydraulic continuity and determining coarse residence time estimates. It is quite plausible that diffusion becomes a dominant process, considering the small advective flow velocities and associated large travel times. The scatter in the data, however, implies that the overall picture is more complex and that several processes, i.e. local differences in recharge concentrations, the amount of dissolution of halite, mixing and diffusion, occur simultaneously at various degrees.

Conceptual model
The position of the two normal Linyanti and Sibbinda faults, together with (relatively sparse) drilling information and geophysical data, reveals that the position of the top of the Karoo basalts varies significantly and abruptly, representing the tectonic structures of a rift zone. Chemical, isotope and hydraulic information indicates that the LKA represents a continuous, extremely old (i.e. >100 ka), but still active hydraulic system within this section of the ORZ. Figure 13 illustrates a conceptual model of how the aquifer system in the ZR might have developed. Accordingly, the LKA originally represents a fluvio-deltaic sedimentary sequence that was deposited prior to rifting by discharge from the northern Kwando and Zambezi rivers. With the onset of rifting associated with the formation of the ORZ-presumably from the Early to Middle Pleistocene-the LKA was tectonically displaced. Rifting  Extensional tectonics in the course of the formation of the ORZpresumably during the Early to Middle Pleistocene-led to the vertical displacement of the LKA. Intermittent lacustrine phases associated with various stages of Paleolake Makgadikgadi probably extended, especially into the southern parts of the study area, thus forming clay-rich confining layers, largely burying the LKA. c As a result of the continuous development of the graben system in the N-S direction during the late Quaternary, the LKA was covered by younger, overall loamier deposits that form the UKA initially occurred along the southern Chobe Fault. This led to the formation of various large paleolake stages in the Makgadikgadi Basin. The paleolake desiccated step by step due to reduction in inflow from its northern tributaries (Moore et al. 2013). Intermittent lacustrine phases probably occurred, especially in the south of the study area, thus forming clayrich confining layers. Repeated desiccation of sections of the paleolakes led to the accumulation of evaporates within the sedimentary sequence. This explains the presence of halite and gypsum and the predominance of sodium at the cation exchange sites. With progressive extension of the graben in the N-S direction and the establishment of the Linyanti and Sibbinda faults, the LKA was further tectonically lowered and covered by younger, overall loamier deposits that form the UKA. To the north of the Sibbinda Fault, the more permeable sands of the LKA are still preserved. This area is characterized by deep-rooting Kalahari woodland. Today, a horst-graben structure is fully developed, representing more or less the final stage of tectonic evolution. The hydrochemical and isotope results provide evidence that the N-S-directed flow is still sustained, which results in a gradual freshening of a more brackish or salty groundwater system. The LKA was probably formed during the Pleistocene and recharged under colder conditions than today, with reduced evaporation. The vegetation was possibly dominated by C4 grasses or macrophytes. Modern recharge of the LKA most likely occurs along the northern graben shoulder via preferential flow through fractures and faults. In a similar way, upward-directed groundwater discharges into the Linyanti Swamps via leakage along the Linyanti Fault. Due to the low gradients and the overall reduced direct recharge under modern climate conditions, however, these processes take place very slowly.

Conclusions
The interdisciplinary investigation of the aquifers in the ZR revealed that the hydrogeological conditions are not simply controlled by modern climatic conditions, but can only be explained by considering changes in the paleoclimatic conditions and the regional tectonic history. The LKA has obviously formed under different paleoclimatic temperature, rainfall and evaporation regimes as well as vegetation cover. Subrecent and recent tectonic rifting largely controls the internal structure, extent and interlinkage with the upper aquifer. Thus, an understanding of the sedimentological and tectonic processes becomes crucial for a more effective and successful exploration and development of such an aquifer. Since characteristically little information on deep groundwater resources is available a priori-due to the risks and costs of drilling deep boreholes-the exploration of deep-seated semi-foss i l o r f o s s i l a q u i f e r s s h o u l d a l w a y s c o m p r i s e multidisciplinary research. As one of the first steps, the use of remote sensing techniques and geophysical surveys should be considered. Of particular importance are the investigation and detection of paleo-drainage conditions, hydraulic boundaries, confining layers and hydraulic windows. Subsurface geophysical data obtained from transient electromagnetic (TEM) surveys are particularly useful for detecting hydraulic boundaries as well as variations in salt content in groundwater. Hydrochemical and isotope results proved very useful in this study, as they provided additional information on groundwater flow and paleo-recharge conditions. The established geochemical evolution down the presumed flow path confirmed hydraulic continuity in a low-gradient basin. Radio-dating revealed that the age of the groundwater of the LKA is much greater than presumed, probably reaching several hundred thousands of years, implying the need for a complete set of different radioactive isotopes in order to cover the appropriate age distribution range. Challenges still exist due to insufficient data and the complex nature of recharge and groundwater flow. Stored groundwater volumes and fluxes of the LKA are still uncertain. Inflow and outflow are assumed to be minor, as modern recharge in the Kalahari Basin is very small overall. Numerical modeling based on the developed conceptual hydrogeological model could provide more insight into current groundwater discharge. Some of the available hydraulic and hydrochemical data, however, may be misleading, because most available boreholes are screened over different sections of the aquifer, and too little is known about the vertical variations in chemical properties within the aquifers. The use of radiocarbon dating in this area is limited due to rock-water interactions and because groundwater ages exceed its inherent dating limit of about 35 ka. Difficulties in dating old groundwater with 36 Cl have included uncertainties about the relative significance of evaporation, mixing and subsurface production in interpreting the measured variations in Cl and 36 Cl. The use of radioactive noble gases such as 39 Ar for younger groundwater near the recharge areas and 81 Kr for very old groundwater is very promising and could substantially improve the interpretation of available 14 C and 36 Cl data.