Impacts of Land Use and Land Cover Change on Non-Point Source Pollution in the Nyabarongo River Catchment, Rwanda

Abstract

The Nyabarongo river catchment in Rwanda has experienced significant changes in its land use and land cover (LULC) in recent decades, with profound implications for non-point source pollution. However, there are limited studies on non-point pollution caused by nutrient loss associated with land use and land cover changes in the catchment. This study investigates the spatiotemporal impacts of these changes on water quality considering nitrogen and phosphorus within the catchment from 2000 to 2020 and 2030 as a projection. The SWAT model was used in analysis of hydrological simulations, while the CA–Markov model was used for the future projection of LULC in 2030. The results revealed (1) the important changes in LULC in the study area, where a decrease in forestland was observed with a considerable increase in built-up land, grassland, and cropland; (2) that the R² and NSE of the TN and TP in the runoff simulation in the catchment were all above 0.70, showing good applicability during calibration and validation periods; (3) that from 2000 to 2020 and looking to the projection in 2030, the simulated monthly average TN and TP levels have progressively increased from 15.36 to 145.71 kg/ha, 2.46 to 15.47 kg/ha, 67.2 to 158.8 kg/ha, and 9.3 to 17.43 kg/ha, respectively; and (4) that the most polluted land use types are agriculture and urban areas, due to increases in human activities as a consequence of population growth in the catchment. Understanding the patterns and drivers of these changes is critical for developing effective policies and practices for sustainable land management and protection of water resources.

1. Introduction

Non-point source pollution (NPS) is a pervasive environmental issue with global significance, impacting water bodies and ecosystems across the planet [1]. Unlike point source pollution, which originates from easily identifiable and regulated discharge points such as factories or sewage treatment plants, non-point source pollution arises from diffuse and often unregulated sources. These diffuse sources include various land use activities such as agriculture, urbanization, deforestation, construction, and transportation [2]. The Nyabarongo river catchment in Rwanda stands as a critical nexus of ecological, social, and economic dynamics, embodying both the country’s natural richness and its developmental aspirations [3]. However, Mudahemuka [4] suggested that rapid urbanization, agricultural expansion, and population growth have caused significant alterations in its land use patterns, consequently affecting the hydrological and environmental integrity of the catchment. Central to these transformations is the proliferation of non-point sources of pollution, which lead to intricate challenges in water resource management, ecosystem health, and sustainable development [5].

Nitrogen and phosphorus play a pivotal role in the process of water eutrophication. These essential nutrients are of significant concern when examining the factors that contribute to the excessive enrichment of water bodies [6]. The presence of these pollutants is attributed to the use of both inorganic and organic fertilizers in high-density farming, urban sewage systems, grazing livestock, and farm disinfectants. With the shift in catchment areas towards increased agriculture, deforested lands, and urban development, the amount of contaminants entering the river has escalated [7]. Studying the key areas of non-point source (NPS) pollution and their impact on river contamination as well as how these factors correlate with land use changes within the watershed is a crucial step. This research is vital for establishing strategies to manage land use practices in the watershed and mitigate the decline in river water quality [8].

Recent studies have highlighted the significant impact of land use and land cover (LULC) changes on non-point source (NPS) pollution across various catchments. For instance, a study conducted by Zena [9] investigated the long-term effects of LULC changes on pollution loads in a catchment in central Ethiopia. The research used hydrological modeling and multivariate statistics to analyze the relationship between LULC changes and NPS pollution, revealing a substantial increase in pollutants like nitrate and phosphorus due to human-induced changes such as deforestation and agricultural intensification. The evolution of potential non-point source pollution risks has also been assessed in the Heihe watershed in northern China, which is similarly affected by land use changes [10]. This research found that urbanization significantly contributed to the increase in pollutant loads, with the most recent LULC conditions showing the highest impact. These studies underscore the importance of understanding and managing LULC changes to mitigate their effects on NPS pollution in catchments.

In the Nyabarongo catchment, various studies have been conducted with different purposes, such as soil erosion and flood assessment [11,12,13], but no case of non-point pollution has been assessed in terms of sediment and nutrient contaminants. This research gap attracted the researchers to conduct this kind of study, which endeavors to delve into the intricate interplay between land use cover change and its ramifications for NPS within the catchment. We used the SWAT model, which is a robust hydrological simulation framework, for elucidating the complex dynamics of land usage alterations with cascading effects on water quality and ecosystem health. Several hydrological models have been established to simulate surface runoff, sediment transport, and nutrient distribution in the catchment [14,15,16,17]. In this case, the SWAT model was selected to investigate NPS in the watershed due to its specific design for use in extensive and intricate basins across extended timeframes [18,19]. This tool is capable of modeling and predicting non-point source pollution from its origin to its entry into water bodies. It forecasts the flow and concentration of pollutants at different locations within a watershed over time, using future land use and land cover data projected by the CA–Markov model. Recently, a number of studies have been used this model to project future urban expansion, assess environmental impacts, and support sustainable development planning. For instance, research has shown the model’s capability to predict LULC changes with high accuracy, using Kappa coefficients to validate the reliability of the simulations [20,21,22]. On the other hand, one of the primary advantages of the SWAT model is its adaptable structure, which enables users to segment a vast watershed into numerous smaller sub-watersheds. Additionally, it facilitates the modeling of diverse land management and usage practices through simple adjustments of parameters [23].

This study was therefore conducted to evaluate NPS in the Nyabarongo catchment in terms of nitrogen and phosphorus pollution and spatial and temporal distribution patterns under conditions of land use change using the SWAT and CA–Markov models. Ultimately, the findings of this research will contribute to the preservation of ecological integrity, the enhancement of water security, and the promotion of resilient socio-ecological systems within the Nyabarongo catchment and similar regions grappling with land use-induced environmental challenges.

2. Materials and Methods

2.1. Study Area Description

The Nyabarongo river catchment, spanning latitudes 1°18′ to 2°34′ S and longitudes 29°5′ E to 30°35′ E, is encircled by the central plateau to the south, the eastern plateau and savannah to the east, and a volcanic range to the north with an estimated length of 151.5 km draining a total area of 8478.24 km² [24]. It encompasses roughly one-third of Rwanda’s land area (Figure 1). The Nyabarongo river’s main feeders include the Mbirurume, Rukarara, and Mwogo rivers on the western bank, while the Mukungwa river and Lake Ruhondo are tributaries on the northern bank. The Nyabarongo river and Lake Muhazi are on the eastern bank. Dominated by a tropical climate, the region receives an average annual rainfall of 1231 mm and maintains an average temperature of approximately 17 °C annually. The elevation within the catchment area ranges from 1342 to 4480 m above sea level, averaging at 2911 m. Dominated by agriculture, 74.85% of the catchment’s total land area is composed of croplands. A lack of sustainable land management practices has resulted in a significant soil erosion rate from these croplands, estimated at approximately 618 tons per hectare per year (also containing nutrients), as reported by Karamage and Zhang [25]. Furthermore, the same authors also noted that the Nyabarongo river, flowing into the Akagera river, channels water towards Lake Victoria. As the largest river catchment in Rwanda, it includes different 24 districts from all provinces countrywide and also encompasses industrial areas with significant urbanization, like Kigali city. However, the area experiences various forms of contamination, including domestic waste, agricultural runoff, and emissions from quarrying, mining, and industrial operations, particularly around Kigali and its environs [11].

2.2. Datasets

The digital elevation model (DEM), land use land cover maps, and soil type maps are among some of the spatial data applied for this study. More prominently, meteorological data were served to simulate the catchment model’s climate and precipitation sources. The coordinate system used here is referred to as WGS 1984. The study area’s coordinates were mapped using the UTM projection system, specifically named WGS 1984 UTM Zone 350 S, reflecting the region’s geographical positioning. The basic information of the data used is presented in the table below.

The DEM served as a crucial component of the SWAT model, providing the foundational layer for defining watershed perimeters, stream networks, and sub-basin formation. In conjunction with additional datasets of land use and soil types, the DEM is instrumental in producing hydrological response units (HRUs). For this study, the study area’s topographic features were derived using a DEM with a 30 m resolution, supplied by NASA’s Shuttle Radar Topographic Mission (SRTM). The Landsat image scenes were obtained from the United States Geological Survey (USGS) website, a freely accessible data portal (http://earthexplorer.usgs.gov/ (accessed on 22 June 2022)), as described in

Table 1. Data sources and basic information.

Input Data	Scale	Source
Land use map	Landsat 8.30 m Resolution	Earth exploration (USGS)
Soil map	30 arc-seconds	FAO-UNESCO Soil Map of the World
DEM	30 m	USGS (United States Geological Survey)
Water quality parameters	Monthly, 2010–2020	Rwanda Water Resources Board (RWB) UNILAK-Environmental Laboratory.
Rainfall and temperature	Monthly, 1982–2020	Rwandan Meteorological Agency

.

The soil type distribution data were provided by the HWSD database built by the Food and Agriculture Organization (FAO) of the United Nations, i.e., the FAO 1990 Soil Classification System [26]. Figure 2b depicts the soil type map and

Table 2. Soil type table.

No	Code	Name		Code	Name
1	Nh7-2/3c	Humic Nitosols	7	Fo96-3b	Humic Ferralsols
2	Tm10-2bc	Mollic Andosols	8	Nd39-3bc	Dystric Nitosols
3	Fo42-3b	Humic Ferralsols	9	I-N-c	Lithosols
4	Tm14-1/2c	Mollic Andosols	10	Fh10-3b	Humic Ferralsols
5	Nh5-2/3c	Humic Nitosols	11	Bh14-3c	Humic Cambisols
6	Fo97-3b	Humic Ferralsols

with its classification correspondence.

The categorization of land use significantly impacts the simulation results of the model. To illustrate this, the years 2000, 2010, 2020, and 2030 were chosen as benchmarks, each marking a distinct phase in land use evolution (

Table 3. Images used in this study.

Year	Path/Row	Acquisition Date	Sensor Type	Spatial Resolution (m)	LULC Name	Source
2000	173/62	27 September 2000	Landsat 7 ETM	30	2000 LULC	USGS
2010	172/61	22 July 2010	Landsat 8 OLI	30	2010 LULC	USGS
2020	172/62	27 September 2020	Landsat 8 OLI	30	2020 LULC	USGS

). Understanding these changes is vital for grasping the model’s performance. Historical land use and cover (LULC) patterns from 2010 to 2020 informed the transition probability matrices, which were essential for the model’s initial periods for the CA–Markov model to predict the 2030 LULC. Therefore, to effectively incorporate these land use classifications into the SWAT model, distinct codes were allocated to various land classes, namely FRST for forestland, RNGE for grassland, AGRR for cropland, URLD for built-up areas, WETL for wetland, and WATR for water bodies. This coding facilitates the model’s ability to accurately represent the diverse landscape and its associated ecological functions.

The four stations were selected to collect the variation climate data in terms of precipitation and temperature in the Nyabarongo catchment.

2.3. SWAT Model Description

In this study, SWAT 2012 software was leveraged for hydrological analysis, while ArcGIS 10.8, enhanced with the Arc SWAT as an extension, was utilized for spatial data management. The river basin was analyzed based on soil characteristics, LULC, and topography to delineate sub-basins into hydrologic response units (HRUs). These HRUs are critical for assessing hydrological properties and are central to the simulation, which is driven by the water balance [27]. The simulation process is divided into two main stages: the in-stream stage (which computes the loadings from each sub-basin throughout the stream network) and the terrestrial stage (during which the SWAT model estimates the contributions of flow, sediment, and nutrients from each HRU). These contributions are then cumulatively assessed at the level of the sub-basin [28]. The simulation is based on the following equation:

$${SW}_t={SW}_0+{\sum }_{i=1}^t\left(Rday-Qsurf-Ea-Wseep-Qgw\right)$$

(1)

where SW_t—final content of soil water, SW₀—daily soil water quality, t—time (days), R_day—daily precipitation, Q_surf—surface runoff content per day, E_a—evapotranspiration per day, W_seep—content of water entering the vadose zone from the soil profile per day, and Q_gw—amount of return flow per day.

Calibration and Validation

The SWAT model’s calibration and validation requires the modification of its parameters to enhance correlation between the model’s predictions and the observed data for water variables. After successful calibration, the model becomes a tool for assessing how various land use strategies affect the quality of water and ecosystem vitality [29]. Calibration of the model was performed within the watershed using the SWAT model feature in combination with the SWAT-CUP tool 2012 [30]. The model outcomes and laboratory information were set to calibrate the model with Sequential Uncertainty Fitting Ver.2 (SUFI-2) in SWAT-CUP 2012 [31]. Before the calibration process, we conducted an examination of the sensitivity of the SWAT model’s parameters.

Table 4. Meteorological stations.

No.	Station Name	Latitude	Longitude	Elevation	Year
1	Gatumba	30.1954	−1.9495	1727	1982–2020
2	Ruriba	30.0139	−1.977	1569	1982–2020
3	Nyange	29.6381	−2.109	2082	1982–2020
4	Gasogi	30.1802	−1.9701	1444	1982–2020

displays parameters that were discovered to be the most sensitive. The parameters that were adopted originated from previous modeling studies in the literature [32,33]. The model was calibrated by the chosen sensitive parameters using accessible time-series data of runoff, total nitrogen, and total phosphorus. Therefore, the calibration process took place over the period from 2010 to 2013, followed by a validation phase that occurred monthly from 2014 to 2015.

The assessment of the model’s calibration and validation performance was conducted using two metrics, the coefficient of determination (R²) and the Nash–Sutcliffe model efficiency coefficient (NSE), which served as the benchmarks for evaluation [33]. The coefficient of determination assesses the degree of linear correlation between observed data and the predictions made by a model. Nash–Sutcliffe efficiency (NSE) is a standardized metric that quantifies the proportion of the variance in the observed data that is predictable by a model by comparing the residual variance to the variance of the observed data [2]. Calculations of NSE and R² are shown below:

$$NSE=1-{\displaystyle \frac{{\sum }_{j=i}^n(X_j-Y_j{)}^2}{{\sum }_{j=i}^n(X_j-\overline{X}{)}^2}}$$

(2)

$$R^2={\displaystyle \frac{[{\sum }_{j=1}^n(X_j-\overline{X})(Y_j-\overline{Y}){]}^2}{{\sum }_{j=1}^n(X_j-\overline{X}{)}^2{\sum }_{j=1}^n(Y_j-\overline{Y}{)}^2}}$$

(3)

where n—duration of flow i and the number of errors; Y_j—the value of data simulated at time j; X_j—the value of data observed at time j; and $\overline{\mathrm{X}}$ and Ȳ—the average observed and simulated values.

The R² is in the range of 0 to 1, among which 0 indicates no relationship and 1 indicates perfect correlation. The NSE is within the range of −∞ to 1. If it is negative, the mean of the observed data surpasses the results predicted by the model simulation [29]. The SWAT-CUP 2012 is served by the above two parameters as an objective function in the validation and calibration of the model.

3. Results

3.1. SWAT Model Performances

A total of 25 sensitive parameters, of which 10 referred to runoff and 25 to water quality (nutrients), were selected for the sensitivity analysis, as shown in Table 4.

Throughout the calibration and validation process, adjustments were made to the model’s parameters to reduce discrepancies between the modeled and actual data for runoff, TN, and TP over a five-year period. The R² values ranged from 0.84 to 0.91, demonstrating a robust correlation between the modeled and measured runoff, which is associated with rainfall, total phosphorus, and total nitrogen data, as illustrated in Figure 3, Figure 4 and Figure 5. This suggests that the model captured a significant proportion above 0.72 of the variances in the observed data, demonstrating a good overall fit. However, it is important to note that R² alone does not provide insights into how well the model performs relative to the variability of the observed data. Therefore, the NSE values, which vary between 0.72 to 0.78, were also considered, as shown in Figure 3, Figure 4 and Figure 5. Values above 0.70 for NSE suggest that the model captures approximately 70% of the observed variability, which is considered satisfactory.

3.2. Evaluation of Land Use and Land Cover Changes

The observation of LULC changes in the catchment was performed using data from three years (2000, 2010, and 2020); supervised classification of Landsat imagery resulted in six classes. The main land types in the catchment were cropland, grassland, and forest, as shown in

Table 5. List of used parameters for sensitivity analysis.

No	Parameter Name	Description	Min_Value	Max_Value	Fitted_Value	Taget
1	a_CN2.mgt	SCS runoff curve coefficient	−0.2	0.2	−0.189	Runoff
2	a_ALPHA_BF.gw	Base flow α coefficient	0	1	0.330	Runoff
3	a_GW_DELAY.gw	Groundwater hysteresis factor	30	450	275.000	Runoff
4	a_GWQMN.gw	Groundwater re-evaporation coefficient	0	2	1.500	Runoff
5	a_ESCO.bsn	Soil evaporation compensation factor	0	1	0.043	Runoff
6	a_SOL_AWC.sol	Soil water availability	0	1	0.203	Runoff
7	a_SOL_BD.sol	Wet capacity of surface soil	0.5	2.5	2.080	Runoff
8	a_SLSUBBSN.hru	Average slope length	10	100	76.900	Runoff
9	a_OV_N.hru	Manning factor for slope diffuse flow	0	100	37.000	Runoff
10	a_LAT_TTIME.hru	Soil flow measurement delay index	0	100	40.333	Runoff
11	a_NPERCO.bsn	Nitrogen permeability coefficient	0	1	0.257	Water Quality
12	a_PPERCO.bsn	Phosphorus permeability coefficient	10	17.5	13.225	Water Quality
13	a_PHOSKD.bsn	Soil phosphorus partition coefficient	100	200	186.333	Water Quality
14	a_PSP.bsn	Index of phosphorus effectiveness	0.01	0.7	0.127	Water Quality
15	a_N_UPDIS.bsn	Nitrogen absorption distribution parameters	20	100	61.333	Water Quality
16	a_P_UPDIS.bsn	Phosphorus absorption distribution parameters	20	100	51.733	Water Quality
17	a_FIXCO.bsn	Nitrogen fixation factor	0	1	0.263	Water Quality
18	a_CH_ONCO_BSN.bsn	Concentration of organic nitrogen in the river	0	100	58.333	Water Quality
19	a_CH_OPCO_BSN.bsn	Concentration of organic phosphorus in the river	0	100	83.667	Water Quality
20	a_ORGN_CON.hru	Organic nitrogen concentration in runoff	0	100	27.000	Water Quality
21	a_ORGP_CON.hru	Organic phosphorus concentration in runoff	0	50	19.167	Water Quality
22	a_BIOMIX.mgt	Biomixing efficiency	0	1	0.757	Water Quality
23	a_ERORGP.hru	Organic phosphorus enrichment rate	0	5	2.450	Water Quality
24	a_POT_NO3L.hru	Nitrate decay rate in potholes	0	1	0.043	Water Quality
25	a_ERORGN.hru	Enrichment rate of organic nitrogen	0	5	4.217	Water Quality

. In the three refence years, cropland occupied 5223.812 km² (62.4%), 5000.115 km² (59.7%), and 6021.218 km² (72.0%), respectively, meaning the catchment is considered an agriculture area. The area of forest has been steadily decreasing from 2603.76 to 1255.32 km², resulting from many activities related to agriculture; the area of built-up land has also been increasing (Table 5).

Table 6. Changes in land use areas from 2000 to 2020.

	2000	2010	2020	2000–2010		2010–2020		2000–2020
Classes	Area_Sqkm	Area_Sqkm	Area_Sqkm	Variation	%	Variation	%	Variation	%
Forestland	2603.76	2777.50	1255.32	173.74	6.67	−1522.18	−54.80	−1348.44	−51.79
Grassland	235.47	251.38	638.05	15.91	6.75	386.67	153.82	402.58	170.97
Cropland	5223.81	5000.12	6021.22	−223.70	−4.28	1021.10	20.42	797.41	15.26
Built-up	78.39	97.03	210.25	18.64	23.78	113.22	116.68	131.86	168.20
Wetland	85.54	106.62	106.02	21.09	24.65	−0.60	−0.56	20.49	23.95
Water	134.77	129.04	130.79	−5.74	−4.26	1.75	1.36	−3.98	−2.96
Total area	8361.75

depicts that forestland underwent a loss of −51.79% and water bodies a loss of −2.96% in recent years, although grassland land, cropland and built-up areas have increased by 170.97%, 15.26% and 168.20%, respectively. The most significant land use change in the period 2000–2020 was observed in grassland and built-up areas, which is linked to the increase in the population of the catchment (Figure 6). This observation goes hand and hand with the exploitation of the forest and water.

3.3. Total Nitrogen Trends in Different Land Use Scenarios

3.3.1. Annual Range Simulation Analysis of Total Nitrogen

This study was carried out using land use data from three separate periods (2000, 2010, and 2020). The nitrogen levels discharged from the watershed exhibited various trends in alterations in random periods between 2010 and 2020. The land use map included data from 2000 onwards, where the lowest count was in 2010 (with the value 15.36 kg/ha) and the highest was in 2020 (with a value of 79.61 kg/ha). The 2010 map showed a lowest value in 2011 (39.8 kg/ha) and a highest value in 2020 (84.51 kg/ha); the last map of 2020 showed a clear increase in nitrogen load in the catchment, with a lowest value of 28.61 in 2010 and a highest value of 145.71 kg/ha in 2020 (Figure 7). The outlet point concentrations revealed a progressive increase, which also caused the ground water quality limit set by World Health Organization (WHO) to be exceeded in 2019 and 2020.

3.3.2. Spatial Distribution Characteristics of Total Nitrogen

The SWAT model was fed data from three land use maps to determine the average annual load for each sub-basin during the years 2000, 2010, and 2020 (Figure 8).

The variation in TN in 2000, 2010, and 2020 was between 2.47 to 4.55 kg/ha, 3.06 to 6.44 kg/ha, and 6.05 to 13.31 kg/ha, respectively, in all sub-basins. Therefore, the highest distribution was observed in sub-basins 8, 11, 12, 13, 16, 20, and 23 with an increasing tendency during the years 2000 to 2020. These sub-basins are located mostly in urban and agricultural areas wherein domestic wastewater and fertilizers can actively contribute to water pollution in the catchment.

3.4. Total Phosphorus Trends in Different Land Use Scenarios

3.4.1. Simulation Analysis of Total Phosphorus in Annual Range

This study used land use data spanning three distinct time frames to analyze the exported phosphorus loads from the watershed. This analysis revealed a varied pattern of phosphorus fluctuations during random intervals between 2010 and 2020. Specifically, the 2000 land use data indicated a lowest phosphorus load in 2010, at 2.6 kg/ha, and a highest load in 2020, at 8.38 kg/ha. Meanwhile, the 2010 data depicted a minimum load of 4.46 kg/ha in 2011, escalating to a peak of 10.8 kg/ha in 2020. The final set of data from 2020 demonstrated a clear upward trend in phosphorus load within the catchment area, starting from 7.76 kg/ha in 2010 and surging to 15.47 kg/ha by 2020. Thus, the concentration levels at the outlet points demonstrated a consistent upward trend, ultimately surpassing the groundwater quality standards established by the World Health Organization (WHO) from 2018 to 2020 (Figure 9).

3.4.2. Spatial Distribution Characteristics of Total Phosphorus

The SWAT model was utilized to process three land use maps, resulting in the computation of the mean annual load for each sub-basin during the years 2000, 2010, and 2020, as depicted in Figure 10. TP levels varied across all sub-basins, with readings ranging from 0.30 to 0.54 kg/ha in 2000, 0.41 to 0.86 kg/ha in 2010, and 0.73 to 1.65 kg/ha in 2020. Notably, sub-basins 3, 4, 8, 11, 12, 13, 16, 20, and 23 exhibited the most significant increases over the two-decade span. These particular sub-basins, predominantly situated in regions of urban development and agricultural activity, are prone to heightened water contamination due to the active contribution of domestic sewage and agricultural fertilizers.

3.5. Estimated Future Loss of TN and TP with the Projected LULC in 2030

The projected LULC in 2030, as provided by the CA–Markov model, displayed increases in built-up and cropland areas, with values of 14.5% and 2.45%, respectively. This expansion is expected to have a consequential impact on existing natural landscapes, with grasslands and forests projected to decrease by 13.87% and 7.63%. These changes are primarily attributed to the demand for land to support agricultural activities and urban development (

Table 7. Changes in land use area from 2010 to 2030.

	2010	2020	2030	2010–2020		2020–2030
Classes	Area_Sqkm	Area_Sqkm	Area_Sqkm	Variation	%	Variation	%
Forestland	2777.4988	1255.3198	1159.4289	−1522.1789	−54.8039	−95.8909	−7.6388
Grassland	251.3760	638.0492	549.5345	386.6733	153.8227	−88.5147	−13.8727
Cropland	5000.1154	6021.2184	6169.3167	1021.1030	20.4216	148.0983	2.4596
Built-up	97.0321	210.2529	241.3488	113.2208	116.6839	31.0959	14.7897
Wetland	106.6249	106.0250	110.3860	−0.5999	−0.5627	4.3611	4.1132
Water bodies	129.0359	130.7897	131.6326	1.7538	1.3591	0.8429	0.6445
Total area	8361.75

).

Therefore, this scenario involved an increase in nutrient loss, as shown by the SWAT model results. The projected LULC showed an increase in concentrations of TN and TP in the context of land use and land cover changes by 2030, which is a significant environmental concern. The estimated future loss, as indicated by an increase of 11.8 to 19.1 kg/ha in TN and 2.2 to 3.1 kg/ha in TP suggests a substantial alteration in nutrient loading, which can have profound implications for ecosystem health and water quality.

Thus, the spatial and temporal analysis from 2020 to 2030 showed an alarming increase in nutrient concentrations, most notably in sub-basins 4, 8, 11, 12, 13,15, 16, and 20 (Figure 11). This trend suggests that ongoing land management methods are in conflict with ecological preservation efforts. Notably, the areas most affected are those with significant human settlement and agricultural activity.

4. Discussion

The research highlighted the dynamic nature of NPS (nitrogen and phosphorus) as nutrients loaded across a landscape over time and influenced by land use changes. By integrating the SWAT model, this study analyzed land use data over three distinct periods, providing a detailed examination of the spatial and temporal variations in phosphorus and nitrogen loads within a watershed [34]. This approach revealed a complex pattern of nutrient changes, which is crucial for understanding the environmental impacts of land management practices not only in the present but also in the future by utilizing projections.

The data in Figure 7 and Figure 9 indicated a significant fluctuation in nutrient levels, with the lowest recorded levels of nitrogen and phosphorus being in 2000 (at 15.36 kg/ha and 2.6 kg/ha) before a stark increase to 79.61 kg/ha, 8.38 kg/ha, respectively. This trend suggests a correlation between land use practices and nitrogen exportation. It is important to consider the types of land use during these periods; for instance, agricultural practices often contribute to higher nitrogen and phosphorus levels due to fertilizer application. Urban development can also impact such pollutant loads through increased runoff and reduced natural filtration systems, as reported by Yazdi [35]. The data from 2010 further supported this, showing an initial low at 39.8 kg/ha and 4.46 kg/ha, escalating to 84.51 kg/ha and 10.89 kg/ha. The LULC 2020 provides a clear visualization of rising nitrogen and phosphorus loads, starting from 58.61 kg/ha and 7.68 kg/ha in 2010 and more than quadrupling to 145.71 kg/ha and 15.47 kg/ha in 2020. Finally, a more concrete situation was observed when using the future projected LULC in 2030, as also used by Gasirabo and Xi [20], which showed a gradual increase in the estimated loss of nutrients from 67.2 kg/ha and 9.3 kg/ha to 158.8 kg/ha and 17.43 kg/ha for nitrogen and phosphorus, respectively. Furthermore, the temporal analysis from 2000 to 2020 with projections for 2030 revealed a worrying escalation in nutrient loads, with the highest distribution in sub-basins 3, 4, 8, 11, 12, 13, 16, 20, and 23 (Figure 8, Figure 10 and Figure 11). This pattern is indicative of the cumulative effects of sustained land use practices that are not in harmony with environmental conservation. The spatial distribution of these loads is predominantly in areas of human habitation and crop production [28].

These findings were alarming as they not only reflect changes in land use but also potential risks to water quality and ecosystem health due to eutrophication [36]. Studies like the one conducted on the Ashe river basin in China have used the SWAT model to assess the impacts of land use changes on nitrogen and phosphorus pollution as non-point pollutants, revealing that land conversion to urban areas and the associated increase in precipitation have led to varying nitrogen and phosphorus exports over the last two decades [33]. Similarly, research in the Chesapeake Bay basin in the USA has shown that agricultural and urban areas have significantly impacted nitrogen transport, affecting water quality and potentially altering the role of nitrogen in sediments as a source or sink of the water body [29]. These findings underscore the importance of integrated land and water management strategies in mitigating the adverse effects of land use changes on nutrient loads and, by extension, on environmental and human health. Understanding the patterns and drivers of these changes is crucial for developing effective policies and practices for sustainable land management and protection of water resources.

Thus, the simulation results indicated that the increase in human settlements and the population directly correlates with a rise in contamination of the Nyabarongo catchment. To mitigate this environmental impact, it is essential to implement specialized sewage treatment facilities in key locations. These plants will play a crucial role in safeguarding the watershed from pollution, thereby actively reinforcing the commitment to environmental conservation.

5. Conclusions

The integration of the SWAT model to assess the impact of land use and land cover change on non-point source pollution has proven to be a valuable tool in environmental management. This study aimed to investigate the influence of LULC from 2000 to 2020 (and predicted conditions in 2030) on the concentration of surface runoff and two important pollutants, total nitrogen and phosphorus, in the Nyabarongo river catchment, which is largely made up of agricultural, urban, grassland, and forest areas. The selection of the watershed was attributed to the dynamic nature of land cover, which has changed in tandem with population expansion and development over time. The findings highlighted increases in pollutants from 15.36 to 158.8 kg/ha kg/ha for TP and from 2.46 to 17.43 kg/ha for TN during the study period. According to above results, both agricultural and urban zones were the areas most vulnerable to NPS pollution, particularly load nutrients. The findings of this research can assist decision makers and land custodians in crafting precise and effective conservation tactics. These strategies will be customized to unique land use and land cover categories and the intensity of non-point source pollution, with the goal of safeguarding water resources and maintaining ecosystem viability. During the study, the researchers met with uncertainties and limitations related to data variability and availability for model calibration and validation due to low quantity of hydrological stations within the catchment. Therefore, the spatial and temporal resolution of the data used may significantly influence the model’s performance, with coarse data possibly leading to underestimation or overestimation of pollution levels. Furthermore, it is recommended that in future research, LULC data are integrated within current best management practice (BMP) allocation frameworks to reflect the varied spatial presence of nonpoint source pollutants across landscapes. Particularly in agriculture and urban locales, where NPS hotspots are widely dispersed, it is essential to examine BMP implementation patterns and pinpoint successful community-based measures for managing NPS pollution. Additionally, assessing the societal effects of BMP implementation may contribute to devising impactful community strategies. These suggestions are intended to deepen the comprehension of NPS pollution management and foster the creation of more efficient BMP distribution methods.