Effects of Land Use and Land Cover Change on Catchment Water Balance in Lagha Bor Catchment, Wajir County Kenya ()
1. Introduction
Land use and land cover changes can contribute significantly to land surface and hydrological alterations. In arid and semi-arid areas, with no reliable surface water systems, groundwater remains the single most important source of this scare resource. However, human activities continue to influence changes through utilization of land and related resources in an attempt to meet their growing needs. While some of the changes can have positive effects, others could lead to irreversible and far-reaching negative effects, such as drying and non-use of boreholes, high salinity of borehole water, degradation of soil and vegetation as well as extinction of palatable vegetation and emergence of foreign and invasive vegetation types. As a result, a plan for timely and accurate collection of land uses and land cover information is necessary to inform water resource managers regarding historical and potential changes in hydrologic fluxes and volumes, and strengthen the positive effects while preventing or reducing the effects of negative effects (Albhaisi et al., 2013) [1] .
Among the most recent methods developed to facilitate accurate and timely collection of land use and land cover change information has been the application of remote sensing in the form of time series landsat imageries (Kumar and Mutanga, 2018) [2] . The change detection and land use monitoring process has been further been simplified and made more efficient through development of google earth engine which is a GIS-web based platform that facilitate archiving of remote sensing data, retrieval, analysis, visualisation and exporting of final products (Sidgu et al., 2018 [3] ; Rahmi et al., 2022 [4] ). The advantage of cloud computing has not only reduced computation time but also limitations associated with computer machine requirements (Rahmi et al., 2022) [4] . The timely access to accurate land use and land cover change information has further reduced the time taken to carry out other activities such as predicting the influences of land use changes on catchment water resources (Albhaisi et al., 2013) [1] .
In ASAL areas, land use and land cover change highly attributed to overgrazing around livestock watering points and growth in settlement which affect vegetation cover and bare-land among other land covers (Egeru et al., 2015 [5] ; Jawuoro et al., 2017 [6] ). Many efforts have been made towards studies on effects of land use and land cover changes on catchment hydrology (Lei et al., 2022 [7] ; Ayivi and Jha, 2018 [8] and Shawul et al., 2019 [9] ). The changes have often been observed to alter the interactions between water balance components including surface flows, infiltration, and percolation among other hydrological components. The incremental and overall effect of the land use changes on catchment hydrology cannot be underestimated. According to Lei et al., (2022) [7] , land use changes were observed to have a high influence on the dynamics of surface and ground water quantity within Stör catchment in of Germany. The changes in areas under pasture, settlements and agriculture negatively affected the quantity of surface and groundwater in the catchment.
Lagha-bor catchment which falls in Wajir County in Northern part of Kenya has been experiencing spatial and temporal land use and land cover changes with potential influence on hydrological variations within the catchment. Considering that pastoralism is the main livelihood activity in the catchment, quantification of the effects of land use changes on the integrated surface and groundwater resources was necessary to avert irreversible changes. The study sought to determine the land use and land cover change dynamics in the catchment and their corresponding effects on water quantity. This was achieved with application of GEE and integrated SWAT + and gwflow model with a main aim of providing bases for policy direction on the best management interventions.
2. Materials and Methods
2.1. Location, Extent and Description of the Study Area
The catchments slopes in a northwest-south east direction with elevation ranging from 1381 m ? 397 meters above sea level with the predominant slope of 0.5% - 2% covering over 90% of the catchment. Lagha-bor is the main ephemeral stream draining the entire catchment from Moyale through Griftu to Diff covering an area of approximately 22,730 Km2. The catchment area accounts for 40% of Wajir County land mass. The area falls within an ASAL environment where the main livelihood activity is nomadic pastoralism where dominant land covers include sparsely distributed shrubs and bare ground. An average annual rainfall of 350 mm is experienced. The upstream area receives the highest rainfall at 500 mm while the lowest parts of the catchment receive as low as 250 mm annually. The maximum and minimum temperatures are 36˚C and 21˚C respectively with an average annual temperature of 28˚C. Figure 1 shows the location of the study area.
2.2. Methodology
2.2.1. Land Use and Land Cover Change Analysis
The study was carried out for the period 1986-2018. It had been established that from 1986, most of the water points initially developed and managed by government were handed over to the community for operation and use. In the period 1986-1990, the only water projects operated by the government were within
![]()
Figure 1. Map of study area. Source: Author―generated from SRTM DEM and ESRI datasets on international boundaries, Counties boundaries, roads, rivers and towns.
Wajir town. The shift in management of water points, led to the change in vegetation management with areas that were considered either wet or dry grazing put under free grazing. The land use and land cover changes could hence bet traced to the changes in water management and arising settlements.
Land use land cover change (LULCC) detection and analysis was carried out in Google Earth Engine (GEE) platform. The GEE platform provides for both open access datasets through retrieval from the archive, cloud computing, analysis and visualisation of results (Kumar and Mutanga, 2018) [2] . The platform operates with the help of javascript code programming.
Appropriate javascript codes were developed to facilitate import of the project area shapefile; retrieval of landsat images along with cloud cover filtering and clipping to the project area; development of classifiers, undertaking supervised classification and carrying out image differencing for land use and land cover change analysis (Kumar and Mutanga, 2018 [2] ; Rahmi et al., 2022 [4] ). Training data and classifiers for the four land use and land cover types namely dense shrubs, sparse shrubs, grassland, and bare-ground/built-up were developed using high resolution google earth embedded in the GEE. In the classification, random forest (RF) classifiers were used to extract land cover types for each landsat image (Rahmi et al., 2022) [4] .
The four land use and land cover types corresponded to the descriptions provided in Table 1. Validation and accuracy assessment was also undertaken in GEE to ensure high level of classification. In the accuracy assessment, the study adopted a minimum of 85% for producer, user and overall accuracy and 0.80 kappa values to achieve the most reliable outcomes (Hütt et al., 2016 [10] ; Mohajane et al., 2017 [11] ).
2.2.2. Predicting the Effects of Land Use and Land Cover Change on Catchment Water Balance Components
The integrated SWAT + gwflow model was utilized to simulate hydrologic processes in the study catchment. The combination of SWAT + and gwflow modflow module helped improve the overall simulation by supporting and better accounting for both surface, sub-surface and deep percolation processes.
![]()
Table 1. Description of the land use and land cover types for the study.
SWAT is considered as one of the few widely and globally applied program with over 20 years of continued development. In the recent past, it was modified and restructured into SWAT + and made available as QSWAT + plugin in QGIS. To overcome limitation and better present the spatial interactions between the various catchment hydrologic processes, the SWAT has seen modifications and restructuring to the current SWAT + version. The restructured SWAT + version also offers better user interface (Abbaspour et al., 2017 [12] and Bailey et al., 2020 [13] ).
While simulating the catchment hydrologic process, SWAT + makes use of digital elevation model (DEM), soil, land use and climatic datasets. The process involves delineation of sub-catchments and creation of hydrological response units that form basic units for simulating changes and management interventions within the catchment. Similar to the original SWAT, the SWAT + simulates for both land and soil hydrologic fluxes such as evapotranspiration, surface flow, percolation, lateral flow and aquifer/groundwater processes among others (Abbaspour et al., 2017) [12] . Its limitation towards deep percolation necessitated integration of gwflow in SWAT + (Bailey et al., 2020) [13] . Figure 2 shows simulation components and interactions between the various catchment hydrologic processes in a SWAT + model. The components are generated as outputs during simulation in a SWAT model.
To improve the representation and simulation of groundwater storage and flow in catchment systems, Bailey et al., (2020) [13] developed the gwflow module for SWAT+. The module is imbedded within the SWAT + code as a subroutine, and links with other hydrologic objects (HRUs, channels, reservoirs) within the catchment to simulate groundwater interactions such as recharge, groundwater-stream exchange, groundwater-reservoir exchange, tile drainage, and saturation excess flows. The gwflow module uses grid of cells (250 m by 250 m was used in this study) to represent individual aquifer control volumes. For each control volume, the groundwater storage are updated during each daily time step
based on the balance of groundwater inflows (e.g., recharge, stream seepage) and outflows for example groundwater discharge to streams, tile drainage, groundwater pumping (Bailey et al., 2020) [13] .
2.3. Integrated Model (SWAT+ and Gwflow)
2.3.1. SWAT + Model Set-Up
A SWAT + model was developed using QSWAT + plugin in QGIS. The interface guides the users in developing the SWAT + model from loading the required input datasets, catchment delineation, creation of HRUs along writing input files and running the SWAT + model in SWAT + Editor 60.04 (Abbaspour et al., 2017 [12] ; Bailey et al., 2020 [13] ). Figure 3 shows the SWATPlusEditor interface in the process of model development. The input spatial datasets including DEM, Land use and Soil were projected into the same co-ordinate system and the look-up tables for land cover and soil maps prepared (Abbaspour et al., 2017 [12] ; Abbaspour et al., 2015a [14] ).
Four different models were set-up based on different land-use/cover maps of 1990, 2000, 2010 and 2020. On the DEM, a default drainage threshold of 620 Km2 was used to discritize the catchment in to 23 sub-basins and a 5% threshold
![]()
Figure 3. Model definition and simulation in SWAT Plus Editor.
applied for on the land-use, soil and slope to generate the HRUs. This resulted in 1160, 1137, 1137 and 1100 HRUs for the 1990, 2000, 2010 and 2020 SWAT + models respectively.
The HRUs are normally filtered based on landuse, soil and slope by applying a threshold for land soil and slope factors. Previous studies have deduced that while thresholds of 5%, 10%, 15% and 20% are applied in order to aggregate small HRUs and increase of computation time, such increase in threshold often lead to over-prediction of daily surface flows by 0.45%, 1.25%, 2.78% 3.83% and 4.81% for 5%, 10%, 15%, 20% and 25% thresholds respectively (Jiang, et al., 2021) [15] . This contributes to reduction of R2 and NSE coefficients hence reducing model accuracy and increasing uncertainty (Jiang, et al., 2021) [15] . While the effect of threshold is considered less on monthly flow measurements as compared to daily flows, higher effects are observed on water quality and sediments (Jiang, et al., 2021) [15] .
To maintain a balance between reduction in computation time and model accuracy for water quality and quantity outputs, the current study chose to adopt a 10% threshold for the land, soil and slope. This threshold was considered optimal owing to the ASAL nature of the catchment where rainfall is highly sporadic. The application of this threshold led to reduction of HRUs through aggregation of small HRUs from 3800 to 1100 which greatly help on improving the computation time. For each of the four models corresponding to the land use maps of 1990, 2000, 2010 and 2020, the simulation time was set from 1986 to 2015 with a warm-up period of 3 years. The warm-up period was used to initiate model stability and hence no results are provided for the period.
2.3.2. Gwflow Model Set-Up
The gwflow model for the study was prepared using the procedure and data types outlined by Bailey et al., (2020) [13] . A 250 m grid cell spacing used resulted to 990,945 Cells and 74 permeability zones. The aquifer saturated hydraulic conductivity (K) for most of the zones were 0.005, 1.007 and 19.2 m/day while the aquifer porosity ranged 0.23 - 0.29. The specific yield was considered at 60% of the aquifer porosity, which gave an estimated value of 0.18. Figure 4 and Figure 5 represent part of the datasets used in the development of the gwflow module.
All the necessary data for the running of gwflow model was prepared and entered in the three main files: gwflow.input, gwflow.hrucell, gwflow.cellhru. As a result, a gwflow module was prepared for each of the corresponding SWAT + model and run with a warm up period of 3 years for which the outputs were not printed and a time step of 0.25 days used to ensure stability of the groundwater balance equation (Bailey et al., 2020) [13] . The simulated groundwater heads were compared to the measured values. The results showed that most of the boreholes with low yields had been sited or drilled on areas with low groundwater volumes.
![]()
(a)![]()
(b)
Figure 4. HRU_grid cells and River_grid cells intersection.
![]()
(a)![]()
(b)
Figure 5. Hydraulic conductivity and Bedrock data.
2.4. Calibration, Sensitivity and Uncertainty Analysis
2.4.1. Model Calibration and Validation
Model calibration is a routine practice in modelling works and involves adjustment of model parameters to a level where the simulated results best matches the observed data. For large catchments, the range for most parameters are not well known owing to spatial variability, measurement errors and complexity of the hydrological process (Moges, et al., 2022) [16] . It is for these reasons that SWAT + models were calibrated before the inclusion of the groundwater flow module (gw) for overall simulation. This was carried out to ensure that the hydrological fluxes simulated well mimicked the observed data.
The traditional practice in calibrating many hydrological models including SWAT model has been the use of measured stream discharge either from a single or multiple hydrological stations within the catchment. Where measured data such stream flow is either unavailable or unreliable, the application of remote sensing derived data has increasingly gained prominence for use in calibration of hydrological models. Remote sensed data especially actual evapotranspiration (AET) and soil moisture has been relied on in the recent past for calibration and validation of SWAT models (Tobin and Bennett, 2017 [17] ; Odusanya et al., 2019 [18] ).
For the current study, once the SWAT models were built, the simulation results were compared with AET data from different remote sensed actual evapotranspiration products which included Global Land Evaporation Amsterdam Model (GLEAM), Moderate Resolution Imaging Spectrometer (MODIS16) and Global Land Data Assimilation System (GLDAS). Among the three datasets considered, Actual Evapotranspiration (AET) from the GLEAM data provided the best comparison between observed and simulated data and was hence selected and used for simulation. In advantage, the AET GLEAM data has previous been used and proved reliable in calibration and validation of SWAT models (Odusanya, et al., 2019 [18] ; Tobin and Bennett, 2017 [17] ).
Calibration for the QSWAT + model was undertaken using Sequential Parameter Estimation (SPE) in SWATPlus-CUP using GLEAM satellite derived Actual Evapotranspiration (AET) in line with the guidelines described within the documentation for the SWATPlus-CUP software (Abbaspour et al., 2017 [12] ; Chunn et al., 2019 [19] ). A ten year data period was used for calibration and validation where two-thirds of the data, 2001-2007 was used for calibration and the remaining one-third, 2008-2010 was used for validation (Abbaspour, et al., 2015b) [20] . SWATPlus-Cup model performs automatic calibration where parameter changes happen internally as the model runs. The SPE program compares the simulated values to the observed values using objective functions such as NSE, KGE, R2 and bR2 statistics. The model was considered sufficiently calibrated when R2 ≤ 1 and 0.5 ≤ bR2 ≤ 0.7 (Arnold, 1998 [21] ; Moriasi et al., 2007 [22] ). Validation was undertaken after calibration period.
The gwflow module model was calibrated using measured groundwater level data from 50 boreholes available from the catchment. The project setup and initial runs were carried out under steady-state conditions, to establish an initial condition before running transient simulations. After the initial run, variable recharge data from corresponding SWAT model was added to the model (Chunn et al., 2019) [19] . The recharge data was subdivided by SWAT HRUs and added to the top layer of the MODFLOW model in the corresponding locations (Chunn et al., 2019 [19] ; Guzman, et al., 2015 [23] ).
2.4.2. Sensitivity and Uncertainty Analysis
In order to determine the effects of land use change on groundwater quality and quantity, a determination was made on the sensitivity of the parameters to be calibrated. This was achieved by running each model using the initial range of 10 parameters provided in Table 2 for 300 simulations. The results of the parameters sensitive to GLEAM AET data are provided in Table 2. The most sensitive parameters were considered as having a p-value lower that 5% and a high absolute value of t-stat (Bennour et al., 2022) [24] . The results showed that parameters most sensitive to AET were moist soil bulk density (bd), soil depth (z), Soil evaporation compensation factor (esco), Available water capacity of soil layer (awc) and Moisture condition II―SCS curve number (CN2) as shown in Table 2. The most sensitive parameters were adjusted to significantly improve the models ability to simulate monthly catchment evapotranspiration (Musyoka, et al., 2021 [25] , Odusanya, et al., 2019 [18] and Moges, et al., 2022 [16] ). Table 2 provides further details on the parameter ranges for the parameters selected for calibration.
In the SWATPlus model, the CN2 governs the generation of surface runoff from HRUs based on the relationship between land use, soil hydrologic group and precipitation (Donmez, et al., 2020 [26] , Moges, et al., 2022 [16] ). Higher values of CN2 signifies reduction in infiltration and hence increased runoff while lower values translated to reduction in surface runoff due to increase in percolation (Moges, et al., 2022) [16] . The results of the calibration showed that an increase of 15% - 25% in CN2 was necessary for all models. This was attributed to the reduction sparse scrubs and increase in bare ground which consequently contribute to increase in surface runoff, less percolation and recharge.
Soil Depth (z) was the second most sensitive parameter for the models. Soil depth determines the level of soil moisture achieved through infiltration and groundwater evaporation. Deeper soils hold more soil moisture while shallow soils offer least soil moisture storage. Soils in the ASAL areas are characterised by shallow depths. The calibration results showed tendency to reduce the soil depth by 40% - 90%. The Harmonised World Soil Database (HWSD) shows that the soil depth of near 1.0 m which indicates that the catchment soils are considerably shallow (Wieder, et al., 2014) [27] . The reduction in vegetation is most likely leaving the soil vulnerable to increased erosion thereby contributing to reduction in soil depth, reduction in infiltration and storage with increase in surface runoff.
![]()
Table 2. Parameter sensitivities at calibration stage-2020 model.
r = relative percent change, v = replacement/absolute change.
The Soil bulk density (bd) defines the relative amount of pore space or porosity in the soil layer (Odusanya, et al., 2019 [18] and Moges, et al., 2022 [16] ). The increase in bulk density decreases infiltration rate of the soil. The high sensitivity of this parameter means that it highly affects the hydrological processes in the catchment. The calibration results for this parameter showed an increase of the bulk density across the catchment by 27% - 33%. The increase in bulk density called for decrease in porosity which consequently led to decrease in lateral and vertical movement of water in the soil layers (Moges et al., 2022) [16] . According to the Harmonised World Soil Database (HWSD), the major soils in the area include Luvisols, Arenosols and Solonetz (Wieder et al., 2014) [27] . These soils are associated with high clay levels which has higher porosity and lower bulk density of 1.23 - 1.49 g/cm3 (Wieder et al., 2014) [27] . The need to increase the bulk density may have been necessitated by reduction in vegetation or increase in bare ground within the Luvisols and Solonetz soil zones.
The available water capacity (AWC) in the soil denoted the capacity of the soil to hold water. Higher AWC meant that the soil has higher ability to hold water hence leading to less surface runoff and percolation. The calibration results showed the need to decrease the awc by 20% - 25% thereby leading to reduced ability of the soil in the catchment to hold water. This also contributed to reduced percolation and ground water recharge in addition to increase in runoff from the catchment. The shallow soil depths and higher bulk density may have contributed to reduction in soil moisture. Similar findings were deduced in related analysis (Moges et al., 2022) [16] .
ESCO was also identified as the third most sensitive parameter for the study models and controls the soil depth needed to meet evaporation demands of the model. A change in this parameter therefore affects the water balance in the catchment. High ESCO values means a smaller depth meeting the evaporation demand leading to increase in surface runoff while low ESCO values allows the model to extract more water from the lower layers in order to meet water deficit in the upper soil layers (Malagò et al., 2015 [28] and Moges et al., 2022 [16] ). Results of the calibration process estimated an ESCO range of 0.85 - 0.92 which was within the recommended range of 0.75 - 1.0. The medium to high ESCO values for the catchment means lower evapotranspiration demand from lower soil layers a factor supported by sparse scrubs in the semi-arid environment.
2.4.3. Assessment of SWAT Model Accuracy
The monthly GLEAM evapotranspiration data was used to calibrate the models. The accuracy of the model was determined by comparing between the simulated and observed results using three approaches namely graphical, statistical and analysis of model water balance (Musyoka et al., 2021) [25] . From a statistical perspective, the model accuracy and uncertainty were evaluated using p-factors and r-factor respectively (Abbaspour et al., 2017) [12] . The p-factor defines the number of observed data points bracketed by the 95PPU while r-factor is defined by a ratio of average width of the 95PPU and standard deviation of the observed/measured data (Odusanya et al., 2019) [18] . Figure 6 shows the p-value and r-factor for one of the four models. Analysis of the p-factor and r-factor values for all the calibrated models deduced that p-factor values ranged between 55% - 70% which meant that close to 70% of the observed data was within the 95PPU while the values for R-factor ranged 1.1 - 1.5 showing a good level of the model’s ability to predict uncertainty.
The use of individual objective function gives calibration results that are conditioned to the respective objective function (Abbaspour et al., 2017) [12] . From the post processing of multi-objective function, the weights for the individual
![]()
Figure 6. Variation of p-factor and r-factor in some target HRUs of the 1990 SWAT model.
objective functions are provided in Table 3. An assessment of the individual objective function observed that the model performance was more than satisfactory under monthly assessment with coefficient values of R2 ≥ 0.4, NSE ≥ 0.45, PBIAS ± 20% and KGE ≥ 0.40 (Moriasi et al., 2015 [29] and Odusanya et al., 2019 [18] ).
In order to incorporate and enhance the advantage of each objective function, a multi-objective function is recommended. For this study, a multi-objective function was used by combining four statistical coefficients namely coefficient of determination (R2), Nash-Sutcliffe Efficiency (NSE), Percent Bias (PBIAS) and Kling-Gupta Efficiency (KGE) (Abbaspour et al., 2015b [20] ; Moriasi et al., 2015 [29] ) as follows:
Multi-objective function = w1R2 + w2NSE + w3KGE − w4PBIAS
where w1, w2, w3, and w4 corresponds to the weights of R2, NSE, KGE and PBIAS provided after post-processing of calibration results for multi-objective function. W1 is taken as equal to 1 while the weight for the rest individual function were achieved by dividing the weight of R2 by that of the specific objective function. An overall multi-objective goal of 1.08086, 1.02002, 1.05315 and 1.0004 were achieved for 1990, 2000, 2010 and 2020 models.
As shown in Figure 7(a) and Figure 7(b), the visual analysis of the resulting graphs for the four models showed that there was good agreement between the model simulations and observed HRU evapotranspiration data. The visual assessment observed that the simulation results well matched the monthly GLEAM evapotranspiration observation data. However, not all peaks in the observed data well match the modelled simulation. This could be attributed some of the land management such as sparse and transitional settlements which are characteristic to the nomadic pastoralism nature of the communities.
1) Simulation and analysis on effects of land use land cover change on groundwater quantity in the catchment
2) Upon setting up the SWAT + gwflow model for simulation, the effects of land use land cover change on groundwater recharge, depth and volume was assessed by varying land use land cover imageries. Borehole pumping data was further incorporated in order to assess the effects of increased abstractions. The simulated groundwater depth was checked against the current water depth in the boreholes.
3) Effects of land use change were assessed by analysing the changes in groundwater quantity under changing land use and land cover. The model outputs in
![]()
Table 3. Weights for the individual objective functions.
![]()
Figure 7. Simulation and observed flows for the calibration and validation phase―Year 2010 model for HRUs 323 and 588.
terms of catchment water balance including surface flows, groundwater recharge, percolation, groundwater transfer to the soils profile among others were assessed from trend and spatial image differencing.
2.4.4. Assessment of Gwflow Model Accuracy
All the four SWAT models namely 1990, 2000, 2010 and 2020 were calibrated for the period 2001 to 2007 and validated for the period 2008 to 2010. The results from SWAT + model formed an input into the gwflow modflow module (groundwater model). This meant that four groundwater models were developed corresponding to the respective SWAT + models due to difference in HRUs.
The results from (gwflow module) run were only considered successful and meaningful after assessing the results from water depth and/or groundwater volume outputs. This resulting curve was assessed to ascertain that a stable point had been attained from which the observed results were both consistent and constant (Bailey et al., 2020) [13] . For, example, Figure 8 shows that the model results had attained a stable state from year 2001 when water depth at the shallow aquifer was 7 m and corresponding groundwater volume was estimated 12.36 meters annually. All the models were observed to have attained trend similar.
The accuracy of the SWATPlus-gwflow model was assessed by determining the error in the catchment and groundwater balance. For the catchment balance, the error was determined by computing the difference between the total water input (precipitation and groundwater transfer) and the total water output (surface runoff, lateral flow, percolation, evapotranspiration, and saturation excess flow). As shown in Table 4, the results of the model accuracies from a catchment water balance perspective revealed that the model error ranged between 0.2%, 0.8%, 0.9% and 1.4% with corresponding accuracies of 99.8%, 99.2%, 99.1% and 98.4% respectively for years 1990, 2000, 2010 and 2020. The high level of model accuracy confirmed that the simulated results were considerably reliable (Bailey et al., 2020) [13] .
For the groundwater balance, the model accuracy was assessed by comparing the results of the difference between groundwater after and before and the sum fluxes including recharge, groundwater evapotranspiration, groundwater flow to surface runoff, surface runoff percolation to groundwater, saturation excess flow (satex) and groundwater flow to soil profile. The results of the model errors were as given in Table 5. Analysis of the model errors showed that the error for the groundwater flow module (gwflow) ranged 0.3% - 3.2%. From Table 5, the model errors of 1.8%, 3.2%, 0.3% and 1.4% with corresponding accuracies of 98.2%, 96.8%, 99.7% and 98.6% restively for years 1990, 2000, 2010 and 2020. This high accuracies confirmed that the model performed well and exhibited good level of simulation accuracy (Bailey et al., 2020 [13] ; Yimer et al., 2022 [30] ).
![]()
Figure 8. Groundwater volume and water depth variations―1986-2015 (2020 model).
![]()
Table 4. Model accuracies for the catchment water balance.
![]()
Table 5. Model accuracies from a groundwater balance perspective.
3. Results and Discussion
3.1. Land Use and Land Cover Change
Figure 9 shows the temporal variation for the predominant land use and land cover types mainly sparse shrub and bare ground in the catchment. The analysis observed that bare-ground cover increased from 35% to 46% in the period 2000-2020. The analysis of the temporal variation showed that area under bare ground cover recorded gradual increase from 1985 to 2020 while area under sparse shrubs depicted a decreasing trend over the years. According to Figure 9, the catchment has recorded an overall decline in area under sparse shrubs by approximately 2680 Km2 which is equivalent to 11.8% in favour of bare ground cover which increased by 2470 km2 corresponding to 10.8%.
While Figures 10(a)-(d) shows the spatial variations in the four land use/cover classes namely dense shrubs, sparse shrubs, grassland and bare-ground for years 1990, 2000, 2010 and 2020 respectively. The most conspicuous classes were sparse shrubs and bare-ground covers. Visual observations of the spatial variations show that much of the changes were in the central part of the catchment especially between Buna and Abaqdere. The most notably and visible change was observed between the imageries of 2010 and 2020 when the increase in bare ground was highest.
![]()
Figure 9. Temporal variation in land cover (bare ground and sparse shrubs).
![]()
Figure 10. Land cover types for the period 1990-2020.
3.2. Spatial-Temporal Variation of Surface and Groundwater Fluxes
The results of main components in the catchment and groundwater balance files were extracted and analysed to determine the trends. These included surface flows, percolation, lateral flows, evapotranspiration, groundwater transfer from shallow aquifer or water table to the soil profiles (gwtranq), saturation excess flow (satex) and groundwater recharge. Figure 11 shows the variations of these fluxes for the period 2001-2013 during which the model was in a stable state. The influence of land use land cover change on surface and groundwater fluxes was assessed by relating the changes in fluxes shown in Figure 11 to the changes land use and land cover as shown in Figure 12.
Figure 11 shows that precipitation and evapotranspiration (ET) accounted for the highest proportion of the fluxes. For the 13 years period (2001-2013), the analysis of the fluxes showed that precipitation remained constant at 9600 mm. This was occasioned by the fact that for each of the four models (1990, 2000, 2010 and 2020), the weather conditions were held constant while the land use and land cover was varied. As per Figure 11, evapotranspiration was observed to have declined from 7748 mm to 6654 mm between year 1990 and 2000 after which it increasing to 8860 mm in 2010 before reducing to 7400 mm in 2020.
The analysis revealed that evapotranspiration decreased with decrease in area under sparse shrubs and only increased if the sparse shrubs were found to increase. The changes in bare-land/ground cover were also found as having an influence on the level of evapotranspiration. Increase in area under bare-ground cover led to a reduction in evapotranspiration where increase in bare ground/land contributed to decline in evapotranspiration and vice versa.
Similar trend was observed while examining change in groundwater transfer to the soil profile. The transfer of groundwater to the soil profile was highest at 5679 mm in year 2010 which corresponded to the same period when the area
![]()
Figure 11. Variation in surface and groundwater fluxes―2001-2013.
under sparse shrubs cover was highest. Groundwater transfer to the soil profile was lowest at 1850 mm in year 2000 during which the catchment recorded the lowest area under sparse shrubs. This deduced that the presence of vegetation greatly influenced the transfers of groundwater from the shallow water tables to the soil profile. Within Lagha Bor catchment, reduction of vegetation cover increases the area under bare-ground leading to reduction in groundwater transfer from the shallow aquifer/water table to the soil profile to meet the evapotranspiration needs.
Land use and land cover changes were also found to influence the amount of surface flows. The least surface flow (1655 mm) was observed in year 2010 when the catchment was observed to have record the highest area under sparse shrubs and lowest area under bare-ground cover. It is during the same period that the highest levels of percolation (3110 mm), lateral flows (480 mm) and saturation excess flow (1280 mm) were high. This confirmed that the amount of vegetation cover played a great role in influencing these fluxes.
Land use and land cover change in lagha bor catchment was also found to have influence on groundwater recharge. The total recharge results provided in Figure 13 corresponds to the deep groundwater recharge results obtained from simulation using the 2010 gwflow modflow model. The results showed that the highest groundwater recharge was reported by the 2010 gwflow model which utilized the outputs from the 2010 SWAT + model that had been set-up using the 2010 land use land cover map.
Groundwater recharge was estimated at 2030 mm, 1900 mm, 2860 mm and 1950 mm in 1990, 2000, 2010 and 2020 respectively. This showed that highest recharge was observed in 2010 when area under sparse shrubs was highest and area under bare-ground cover was lowest. This confirmed that indeed land use and land cover changes affected recharge in the catchment. Similar findings were reported by Shamsuddaha et al., (2011) [31] and Siddik et al., (2022) [32] who observed that land use and land cover change was one of the main anthropogenic activities affecting groundwater recharge in northwestern Bangladesh.
![]()
Figure 13. Variations of groundwater with land use and land cover changes.
Figure 12 shows the spatial variations in groundwater recharge between 1990 and 2020. It is clear that over 80% of the catchment experience groundwater recharge of between 1 - 50 mm, with the rest area reporting recharge values of between 51 - 200 mm. this was characteristic of ASAL where rainfall is low, coupled with sparse shrubs and high bare-ground cover. According to MacDonald et al., (2021) [33] annual groundwater recharge in most semi-arid areas in Africa range between 60 - 200 mm hence confirming that the groundwater recharge results obtained were within acceptable range. In other areas, studies by Shamsuddaha et al., 2011 [31] and Siddik et al., 2022 [32] reported that groundwater recharge varied spatially between 200 mm to 512 mm in Banglashes. Assessment further observed that the spatial variations in groundwater recharge varied greatly depth to the bedrock layer. Areas with shallow depth recorded low levels of recharge while area with higher depths recorded more recharge.
Analysis of changes in groundwater volume for the period 2001-2013 showed that groundwater volume increased with increase in area under sparse shrubs while at the same time reduced with increase in bare-ground cover. Figure 13 shows that groundwater volume was estimated at 58,670 m3 in 1990 when area under sparse shrubs was 12,690 Km2 but was found to have reduced to 58,335 m3 in year 2000 when the area under sparse shrubs reduced to 11,745 Km2. A reduction of sparse shrubs by 945 Km2 led to reduction in groundwater volume by 335 m3.
In 2010, the groundwater volume increased to 58,750 m3 when the area under sparse shrubs increased to 13,065 Km2 before reducing to 58,555 m3 in 2020 in line with reduction of area under sparse shrubs to 11,130 Km2. These changes are explained in Figure 14. This meant that incorporation of measures that lead to improved percolation would consequently improve recharge and groundwater volume. Similar trend was observed when evaluating influence of land use and land cover change on groundwater recharge and percolation.
Spatial changes in recharge between 1990 and 2020 were obtained by image differencing. An optimized hot spot analysis was carried out on the resulting images to determine areas with significant changes as reflected in Figure 14(a) and Figure 14(b). The assessment showed that there was notable decline in significant recharge hot spots at 99% confidence levels to fewer recharge hot spots at 95% confidence levels.
Further visual inspection of Figure 14(a) and Figure 14(b) showed more significant hot spot areas at 95% confidence level as a result of increase in areas exhibiting significant decrease in recharge especially between Wajir town, Eldas and Ogralle. The diminish of significant recharge hot spot areas and the increase of significant recharge cold spot area could only be associated with changes in land use and land covers.
Figure 15 shows the spatial variation of groundwater volume in the catchment. Overly of the borehole location data on groundwater volume maps revealed that over 80% of the drilled were sited on areas with low groundwater
![]()
(a)![]()
(b)
Figure 14. Changes in groundwater recharge between years 1990-2020.
![]()
Figure 15. Spatial variations of groundwater in the catchment for the period 1990-2020.
![]()
Figure 16. Pumping effect on groundwater heads.
volume. This hence informed the reason why most of the drilled boreholes recorded low yields especially for the areas around Dambas, Sarman to Buna. Low groundwater volume was also influenced by the depth to the bedrock which consequently affects the size of the saturation depth in addition to specific yield. Low groundwater volumes were observed in areas with shallow depths to the bedrock.
The saturation depth is computed as difference between the groundwater head and depth to the bedrock. Results further show existence of small pockets of varying low to medium groundwater volume potentials. Similar finding were reported through investigations by Water Resources management Authority in that the groundwater system in the catchment comprises of small isolated groundwater pockets. The system is controlled by localised weathered basement systems surrounded by non-weathered basement hence allowing for no interaquifer movement (WRA, 2020 [34] ; Kuria, 2013 [35] ). According to Olago, (2019) [36] , presence of local variations in the character of the aquifers within Wajir which highly contributed to variance in borehole depths, water levels, yield/groundwater volume and water quality. The challenge of low groundwater volumes is compounded by the fact that the largest part of the catchment experience recharge rates of between 5 - 10 mm as shown in Figure 16.
Lagha-bor has been experiencing increased and continuous abstractions through pumping. With exception of some few areas around Dambas where the saturation depth and groundwater volume are relatively low (less than 50,000 m3), most of the areas within the catchment would have groundwater potential volumes ranging between 900,000 m3 to 1.36 million cubic meters of water (see Figure 16). Less impact is likely to be felt even with increased pumping. This is due to the nature of the aquifers in Wajir. With more isolated small aquifers, pumping is only likely to affect the individual small aquifers other than the impacts of the same over the entire catchment. While investigating the impacts of intensive groundwater abstraction on regional aquifer system, Shamsuddaha et al., 2011 [31] arrived at similar findings in that while the groundwater levels varied temporally and spatially in the north-eastern region of Bangladesh, intensive abstractions along with land use and land cover changes largely influenced the groundwater depths.
4. Conclusions
This study provided useful information on how the land use and land cover changes are influencing groundwater quantity and quality in lagha bor catchment. It hence forms important scientific baseline information towards planning for management of water and vegetation cover as well as informing future developments within the catchment. The assessment established that evapotranspiration, percolation, recharge and groundwater transfer from the soil profile to the root zone reduced with decrease in vegetation cover or with increase in bare-ground cover. The least surface flow (1655 mm) was observed in year 2010 when the catchment was observed to have recorded the highest area under sparse shrubs and the lowest area under bare-ground cover. It is during the same period that the highest levels of percolation (3110 mm), lateral flows (480 mm) and saturation excess flow (1280 mm) were realised.
This confirmed that the percent of vegetation cover played a great role in influencing these fluxes. For example, the transfer of groundwater to the soil profile was highest at 5679 mm in year 2010 which corresponded to the same period when the area under sparse shrubs cover was highest. Groundwater transfer to the soil profile was lowest at 1850 mm in year 2000 during which the catchment recorded the lowest area under sparse shrubs. This deduced that the presence of vegetation greatly influenced the transfers of groundwater from the shallow water tables to the soil profile. An assessment of changes in groundwater volume between 1990 and 2020 revealed that variations in land use and cover contributed to changes in groundwater volume. The groundwater volume was found to increase with increase in area covered with sparse shrubs while it consequently decreased with increase in bare-ground cover.
Acknowledgements
We are grateful to Wajir County government and more particular to the water department for the immense support in sharing information related to state of water resources in the county, developments and management of water resources in the county. Your support in planning and taking part during the field visits was most appreciated. I want to most sincerely thank the operators at various water resource points especially the boreholes for their invaluable time, discussions and inputs during the field visits. We appreciate your support.