Estimation of nutrient loads in an agricultural area: Case study of the Lobo Basin in Côte d'Ivoire (West Africa)

1. Introduction

FAO statistics show that agriculture, hunting, fishing and forestry supported 2.57 billion people at the start of the new millennium. This figure represents 42 percent of humanity, including those who work within these sectors and their families. Agriculture is essential for the economy in most developing countries and even in industrialized countries like the USA, where agricultural exports reached $177 billion in 2021[1]. Today, few countries are able to experience significant economic growth without agricultural development.

In Côte d'Ivoire, one estimate about 60% the population living in rural areas, and the agricultural sector occupies 30.7% of the Gross Domestic Product (GDP). Thus, the State of Côte d'Ivoire has set up several Government agencies structures such as National Rural Development Support Agency (ANADER in french), the National Center for Agronomic Research of Cote d'Ivoire (CNRAin french), the Marketing of Food Products Office (OCPV in french) and the Development of Rice Cropping National Office (ONDR in french) to oversee and assist with agricultural management and development for sustainable agricultural growth and self-sufficiency of foods. Faced with the rising demand linked to population growth, the Government of Côte d'Ivoire uses these agencies to encourage traditional farmers and industrialists to use chemicals to increase their agricultural productivity. However, a lack of training linked to insufficient resources in the sustainable use of chemicals and loose policy on environmental protection and quality standards have led to increased usage of agricultural inputs. This irresponsible usage means that these inputs are now polluting the water bodies and river networks of Côte d'Ivoire, causing major eutrophication problems. Several previous published studies [2,3,4,5] demonstrated eutrophication of the Lobo water body based on analyzed biological and physico-chemical parameters. These studies, although important didn’t alow us to understand deeply the eutrophication phenomenon in the Lobo basin because the samples were punctual at the level of the reservoir but failed to relate critical sources in the drainage areas to the water quality problem. To address the lack of on-site data in that context, methodologies based combining hydrological models and Geographic Information Systems (GIS) are supposed to be more efficient. Watershed Hydrological models are commanded by mathematical formulas in link with physically-based process and conceptual models. The basin hydrologic functioning includes an analogy and a concept. In the analogy, soils and groundwater are considered as reservoirs whose rate of emptying depends on the filling. The physically-based models are directly linked to physical processes [6,7]. The principle of calculation is based on precise physical and empirical representations of natural processes. Empirical models [8] such as the USLE (Universal Soil Loss Equation) erosion model [9], derived from laboratory or field experiments. These models are global: this the case of the GR4 model which can represent the watershed in a single entity for global assessment [10]. However, global-scale models are not suitable for studying the hydrological processes of an agricultural watershed because they do not consider the watershed variability spatial of phenomena like runoff or erosion. So, they cannot take into account the land uses effects in processes. Thus, semi-distributed models such as the Soil and Water Assessment Tool (SWAT), come in addition to overcome this lack by now, considering topography, land cover, soil and soil management practices in the basin.

This study aims to assess the nutrient loads from agricultural practices in the Lobo watershed to understand the spatio-temporal evolution of nutrient transport processes using the SWAT model. The assessment is based on the amounts of fertilizer applied to each crop.

2. Materials and Methods

2.1. Materials

2.1.1. Study area

The Lobo watershed at Nibéhibé is located in west-central Côte d'Ivoire between 6°00 and 7°00 West longitude and 6°54 and 8°00 North latitude. It covers about 6,442.66 km² as area (Figure 1). The plateaus about 200 to 400 m altitude compose the relief. This makes it a little bit contrast [11]. The underlying geology is primarily composed of granites. The peaks of precipitation are observed in the month of September with 275 mm and in May with 150 mm. This makes it possible to define two seasons: a dry season and a rainy season. The dry season is from November to February while the rainy season is observed from March to October. The average annual rainfall recorded in the study area was 1335 mm from 1991-2020. This value of rain amount recorded annually is lower than reported [12] in that findings showed 1437,4 mm as annual rainfall from 1943 to 2010 in the same basin using the GR4J model, or a decrease of 7.14%.

The study area's vegetation is composed of dense moist semi-deciduous forest and cleared mesophilic forest. Ferralitic soil, strongly or moderately desaturated, governs this area. The main activity is agriculture based on coffee, cocoa, rubber, oil palm and food crops [13]. the drinking water supply comes from the Lobo reservoir exploited by SODECI.

2.1.2. Data

Data used in this study were Digital Elevation Model (DEM), land use map, soil map, soil physicochemical and biological properties, climate and agronomic data. The Figure 2 shows essential input data for the SWAT model.

Digital Elevation Model (DEM)

The Digital Elevation Model (DEM) of 30 m x 30 m resolution, in the Northern Hemisphere projection UTM Zone 30, downloaded from the website https://search.earthdata.nasa.gov/search has been used to extract the hydrographic network and delineate the watershed.

Land use map

The land use map has been obtained throughout four (4) February 22, 2018 sentinel-2 satellite images treatment using ENVI 4.7. These images include four (4) bands that are 2, 3, 4 and 8 and were obtained from the website https://scihub.copernicus.eu/dhus/#/home. For this study we considered Five classes of land use: forest evergreen (11%), degraded forest (9%), agriculture and fallow (77%), urban (0.5%), and water (2.5%).

Soil data

To run the SWAT model, we need texture, available water content, hydraulic conductivity, bulk density, and organic carbon content of different soil layers. These data are available in the FAO database [15] and have been completed by physicochemical properties obtained after field sampling and in-laboratory analyses. The dominant soil of the basin is the sandy-clay-loam about 38%, followed by the sandy-loam soil (32%) of the basin.

Hydro-climate data

The climate date at a daily time step were obtained from the swat model website http://globalweather.tamu.edu/. These data include minimum and maximum temperatures, precipitation, relative humidity, solar radiation and wind speed. These data were The use of these data is justified by the fact that in the watershed, there are few data available that can allow adequate hydrological modeling of the watershed. Some authors such as Mbungu and Kahaigili [16] believe that global gridded climatological datasets (GGCDs) could be combined with hydrological modelling to understanding hydrological processes in data-scarce environments. Therefore, for this study, we have used Eleven (11) climatic stations over the period 1979-2014. As for the flow data, stream gauge data were collected from the National Drinking Water Office (ONEP in french) for the period from 1981 to 1994 at Nibehibé hydrometric station.

Agronomic data

The types of fertilizers, the doses applied to each type of crop, the types of crops, agricultural practices and agricultural calendars on the watershed were used in this study (http://www.cnra.ci/listefiche.php).

2.1.3. Computer software

The computer tools used in this study consist of QSWAT 1.9 in the QGIS 2.6.1 interface, SWATCUP 2019, ENVI 4.7 and the WGN program. QSWAT 1.9 made it possible to delimit the watershed into sub-watersheds and to extract the hydrographic network and to simulate the water flow and the nutrients transfer. SWATCUP 2019 was used to calibrate the water flow. ENVI 4.7 has been useful in the processing of satellite images for obtaining land use. As for the WGN program, it made it possible to calculate the statistical parameters from the climatic data.

2.2. Methods

The methodology based upon the SWAT (Soil and Water Assessment Tool) model was structured into two main stages: flow calibration and nutrient estimation.

2.2.1. Flow calibration

Model and software description SWAT Model description

SWAT is a semi-distributed agro-hydrological model [17] with a daily time step. SWAT is used for evaluating water quality and quantity, agricultural practices management, surface and subsurface water management and sediment, nutrient and pesticide transfer. SWAT is used worldwide by a large scientists community.

This model is based on the water balance equation which is as follows [18]:

SWt = SWo + Σ (Rday +IRR– Qsurf – Ea – Wseep), (Eq.1)

where SWt is the final soil water content (mm) of the day, SWo is the in soil initial water content (mm), Rday is the daily rainfall (mm), IRR is irrigation volume added to the soil (mm), Qsurf is the surface runoff volume (mm), Ea is the actual evapotranspiration (mm), and Wseep is percolation loss from the soil profile into the shallow aquifer (mm).

SWAT-CUP and SUFI-2 algorithm

Uncertainties may occur while setting up and running the model. Uncertainties come either from measurement errors or from the model itself [19,20,21]. To be used efficiently, the model’s performance must be tested. This study tested the correlation between observed data and simulated outputs statistically and graphically. Thus, Swat-cup (SWAT Calibration Uncertainty Program) was used to optimize the flow model. The program is the most widely used tool by the SWAT community. SWAT-CUP allows to perform calibration, uncertainty and global sensitivity analysis automatically (Table 1). Several methods are integrated into the SWAT-CUP program. We have Sequential Uncertainty Fitting version 2 (SUFI-2) [22,23], Generalized Likelihood Uncertainty Estimation (GLUE) [24], Particle Swarm Optimization (PSO) [25], Parameter Solution (PARASOL) [26] and Markov Chain Monte Carlo (MCMC) [27].

Because of its simplicity in its use, the SUFI 2 method has permitted to perform global parameter sensitivity analysis, uncertainty analysis and calibration [28].

Global sensitivity analysis

The analysis of global sensitivity begins with the definition of the objective function. Objcetive functions wthin SUFI2 include R², Chi2, NS, and R² multiplied by the line regression coefficient b, bR², and SSQR (sum of squared errors) coefficients. But for this study, NS and R² coefficients were the objective functions used. NS is used to analyse the strength of the model predictions and while R² coefficient is related to the correlation between simulations and observed values. Then, the values of the parameters to be optimized, are chosen based on Equation 2 [29]:

b_{j,abs _min} ≤ b_j ≤ b_{j,abs_max} , j=1…m (Eq.2),

where b_j is the j^th parameter and m is the number of parameters to be estimated.

The parameter global sensitivity in each simulation is accounted by the sensitivity matrix J, of g (bj) [29] :

J_ij = ∆_gi/∆b_j i=1…C₂ⁿ (Eq.3),

where bj is the parameter, g is global sensitivity and i is all possible combinations of two simulations and j is the number of columns (number of parameters).

Latin Hypercube sampling is the sampling method used in this study. It makes it possible to calculate the global sensitivity by calculating a multiple regression system of the parameters considered [29]:

$\mathbf{g}=\mathsf{\alpha }+{\sum }_{\mathbf{i}=1}^{\mathbf{m}}\mathsf{\beta }\mathbf{i}\mathbf{b}\mathbf{i}$ (Eq.4),

where $\mathsf{\alpha }=\frac{{\mathsf{\sigma }}_{\mathbf{s}}}{{\mathsf{\sigma }}_{\mathbf{m}}}$, $\mathsf{\beta }=\frac{{\mathsf{\mu }}_{\mathbf{s}}}{{\mathsf{\mu }}_{\mathbf{m}}}$. ${\mathsf{\sigma }}_{\mathbf{s}}$ and ${\mathsf{\sigma }}_{\mathbf{m}}$ are considered to be the respective standard deviations of the simulated and observed data; ${\mathsf{\mu }}_{\mathbf{s}}$and ${\mathsf{\mu }}_{\mathbf{m}}$ are the averages of the simulated and observed data.

The relative impact of each parameter on the flow is identified from a statistical test. This impact is noted bi. The global sensitivity is related to the variation of the objective function which depends on the observed variation of the parameter considered when all parameters change.

Uncertainty analysis

P-factor and R-factor are used to calculate the uncertainties at 2.5% (Xl) and 97.5% (Xu) percentiles of the cumulative distribution at each point. P-factor indicates the percentage of observed data bracketed by 95PPU. R-factor is the average thickness of the 95PPU band obtained from the standard deviation of the observed data [20].

Calibration analysis

We chose NSE, R² and PBIAS as objective functions to analyze the performance and accuracy of the simulations. The PBIAS measures the percentage of bias. When PBIAS is 0, the simulation is said to be accurate. A positive PBIAS shows an underestimation of the model, while negative values indicate an overestimation of the model bias [30]. Its use makes it possible to refine and make more precise the predictions of the model [31] in [32,33,34,35,36,37,38,39,40,41,42,43]. NSE (equation 5) varies from -∞ to 1 (for a strong link between observed and simulated values). Equation (6) is used to calculate PBIAS. Moriasi et al. [35] made the classification of model performance:

PBIAS <± 10: very good performance model;

±10<PBIAS<±15: good performance model;

± 15 <PBIAS < ± 25: satisfactory performance;

PBIAS > ± 25: unsatisfactory performance.

R² is calculated according to equation (7); it varies from 0 to 1, 1 representing a perfect simulation [44]. A value greater than 0.5 is often accepted because it expresses a good restitution of the observed data [32].

$\mathbf{N}\mathbf{S}\mathbf{E}=1-\left[\frac{\sum _{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}\right)}^2}{\sum _{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-{\mathbf{Q}}^{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{n}}\right)}^2}\right]$(Eq. 5)

$\mathbf{P}\mathbf{B}\mathbf{I}\mathbf{A}\mathbf{S}=\left[\frac{\sum _{\mathbf{i}=1}^{\mathbf{n}}\left({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}\right)\mathrm{*}100}{\sum _{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}\right)}^2}\right]$(Eq. 6)

$\mathbf{R}²=\frac{(\sum _{\mathbf{i}=1}^{\mathbf{n}}\left({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}\right)({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}-\overline{{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}}))²}{{\sum }_{\mathbf{i}=1}^{\mathbf{n}}({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-\overline{{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}})²{\sum }_{\mathbf{i}=1}^{\mathbf{n}}({\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}-\overline{{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}})²}$ (Eq. 7)

with Q_i^obs observed flow, Q_i^sim simulated flow, Q^mean average observed flow, $\overline{{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}}$ average observed flow and $\overline{{\mathbf{Q}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}}$ average simulated flow.

Model setup

The model implemented follows the methodology adopted by Koua et al (2019) [45]. The Lobo watershed in Nibéhibé covers an area of 6,442.66 km². The delimitation of the watershed and the extraction of the hydrographic network was possible thanks to the processing of the DEM under QSWAT. The D8 algorithm [46] integrated in QGIS makes it easy to extract watershed boundaries and river network using QSWAT. This step is followed by the integration of meteorological data, land use, pedology and their physico-chemical properties. Then, the hydrological response units (HRU) combining land use type, soil type and slope in a sub-basin were calculated. HRUs are the basis for calculating for the SWAT model. They highlight the differences in evapotranspiration and other hydrological conditions between different types of land cover. the multiple HRU method of the SWAT model made it possible to constitute 163 HRUs. After that, QSWAT established the various input data tables. This made it possible to integrate information relating to the watershed, pedology, hydrometeorology, the HRUs themselves, water resources and their use, watershed management, wetlands, septic tanks, different reservoirs on the watershed, and the watershed master file. At this stage, the hydrological simulation of the Lobo watershed can be done.

Streamflow calibration process

The flow calibration at a monthly time step was possible by using observed data at the Nibéhibé hydrometric station using the SUFI-2 method [22,23] in the SWAT-CUP program [47]. For this study, ten (10) flow parameters have been used. These parameters were the CN2 expressing the humidity conditions II, the quantity of water available in the soil SOL_AWC, the ESCO factor which translates the compensation coefficient of the soil evaporation, RCHRG_DP which is the fraction of percolation water from the deep aquifer, the re-evaporation coefficient of groundwater GW_REVAP, ALPHA_BF (Baseflow alpha factor), the water depth threshold in the shallow aquifer depth GWQMN required for return flow to occur, the water depth threshold in the shallow aquifer REVAPMN for revaporation to occur, the groundwater delay time GW_DELAY, the CANMX (Maximum storage of canopy). Flow calibration was performed over the period 1981 to 1994. A total of 500 iterations were performed during the monthly water flow simulation process. Then, the Nash–Sutcliffe coefficient (NSE), the percentage of bias (PBIAS) and the R² determination coefficient were calculated in order to judge the robustness of the SWAT model on this watershed.

2.2.2. Nutrient loads estimation

While reservoir eutrophication processes are significantly influenced by nutrient contents in the water, there exists little information upon the reservoir water quality. Because of the lack of observation data on nitrogen and phosphorus, the estimation of nutrient inputs to the reservoir was made on the basis of the quantity of NPK fertilizers applied to crops in the watershed because these nutrients come from NPK chemical fertilizers [3,4]. Data relating to fertilizers, such as the amount, duration, and frequency of application, were integrated into SWAT according to crops by sub-basin and HRU. The mineral fertilizers are of the xN-yP-zK formula, x, y and z representing respectively the proportions in percentage (%) of N, P and K (Potassium) in 1 kg of fertilizer used. If the fertilizer formula is known, the quantities of N, P and K are determined (Eq. 8 & 9) [48]:

${\mathbf{N}}_{\frac{\left(\mathbf{K}\mathbf{g}\right)}{\mathbf{K}\mathbf{g}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)}}=\mathrm{\%}\mathbf{N}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)\times 1\mathbf{K}\mathbf{g}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)/100$ (Eq. 8)

${\mathbf{P}}_{\frac{\left(\mathbf{K}\mathbf{g}\right)}{\mathbf{K}\mathbf{g}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)}}=\mathrm{\%}\mathbf{P}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)\times 1\mathbf{K}\mathbf{g}\left(\mathbf{f}\mathbf{e}\mathbf{r}\mathbf{t}\right)/100$(Eq. 9)

For this study, only N and P quantities were estimated.

From the soil, the crops take Nitrogen in the form of nitrates [48]. Nitrate transportation within the watershed is done via surface runoff, lateral flow, or percolation. The amount of nitrate contained in runoff is assessed from the nitrate concentration in moving water. This concentration is a proportion of the nitrate mass in the runoff volume [48]. It is calculated according to Eq. 10:

${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{N}\mathbf{O}3,\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}=\frac{{\mathbf{N}\mathbf{O}3}_{\mathbf{l}\mathbf{y}}.(1-\mathbf{e}\mathbf{x}\mathbf{p}⌈\frac{-{\mathbf{w}}_{\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}}{\left(1-{\mathsf{\theta }}_{\mathbf{e}}\right).{\mathbf{S}\mathbf{A}\mathbf{T}}_{\mathbf{l}\mathbf{y}}}⌉}{{\mathbf{w}}_{\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}},$(Eq. 10),

where ${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{N}\mathbf{O}3,\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}$ is the concentration of nitrate in mobile water for a given soil layer (Kg N/mm H₂O), ${\mathbf{N}\mathbf{O}3}_{\mathbf{l}\mathbf{y}}$ expresses the quantity of nitrate (Kg N/ha), ${\mathbf{w}}_{\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}$ is the amount of mobile water in the layer (mm H₂O), ${\mathsf{\theta }}_{\mathbf{e}}$ is the porosity fraction on which the anion exclusion depends, ${\mathbf{S}\mathbf{A}\mathbf{T}}_{\mathbf{l}\mathbf{y}}$ is the water content when the soil layer is saturated (mm H₂0).

The amount of mobile water in the soil layer is assumed equal to that lost by runoff, lateral flow or percolation and can be estimated by equations 11 and 12.

${\mathbf{w}}_{\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}={\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}+{\mathbf{Q}}_{\mathbf{l}\mathbf{a}\mathbf{t},\mathbf{l}\mathbf{y}}+{\mathbf{w}}_{\mathbf{p}\mathbf{e}\mathbf{r}\mathbf{c},\mathbf{l}\mathbf{y}}$for top 10 mm (Eq. 11)

${\mathbf{w}}_{\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}={\mathbf{Q}}_{\mathbf{l}\mathbf{a}\mathbf{t},\mathbf{l}\mathbf{y}}+{\mathbf{w}}_{\mathbf{p}\mathbf{e}\mathbf{r}\mathbf{c},\mathbf{l}\mathbf{y}}$ for lower soil layers, (Eq. 12)

${\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the depth of surface water runoff on the given day (mm H₂O), ${\mathbf{Q}}_{\mathbf{l}\mathbf{a}\mathbf{t},\mathbf{l}\mathbf{y}}$ is the height of water escaping from the layer laterally (mm H₂O), ${\mathbf{w}}_{\mathbf{p}\mathbf{e}\mathbf{r}\mathbf{c},\mathbf{l}\mathbf{y}}$ is the amount of percolation water in the underlying soil layer on the given day (mm H₂O).

Surface runoff transports nutrients to the top 10 mm above the soil. Equation 13 was used to calculate the amount of nitrate in runoff.

${\mathbf{N}\mathbf{O}3}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}={\mathsf{\beta }}_{\mathbf{N}\mathbf{O}3}.{\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{N}\mathbf{O}3,\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}.{\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}},$ (Eq. 13)

${\mathbf{N}\mathbf{O}3}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the quantity of nitrate retained during surface runoff (Kg N/ha), ${\mathsf{\beta }}_{\mathbf{N}\mathbf{O}3}$ represents the nitrate percolation coefficient, ${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{N}\mathbf{O}3,\mathbf{m}\mathbf{o}\mathbf{b}\mathbf{i}\mathbf{l}\mathbf{e}}$ expresses the concentration of nitrate in the runoff water for the top 10 mm of soil (Kg N/mm H₂O), ${\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the runoff height of the given day (mm H₂O).

Surface runoff allows organic nitrogen mostly attached to sediments transfer to the main channel in the watershed. The amount of organic nitrogen transported with thesediments to the watercourse is calculated using Equation 14 [49,50]:

${\mathbf{o}\mathbf{r}\mathbf{g}\mathbf{N}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}=0.001.{\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{C}}_{\mathbf{o}\mathbf{r}\mathbf{g}\mathbf{N}}.\frac{\mathbf{s}\mathbf{e}\mathbf{d}}{{\mathbf{a}\mathbf{r}\mathbf{e}\mathbf{a}}_{\mathbf{h}\mathbf{r}\mathbf{u}}}.{\mathsf{\epsilon }}_{\mathbf{N}\mathbf{s}\mathbf{e}\mathbf{d}}$(Eq. 14),

where ${\mathbf{o}\mathbf{r}\mathbf{g}\mathbf{N}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the amount of organic nitrogen contained in surface runoff transported in Kg N/ha, ${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{C}}_{\mathbf{o}\mathbf{r}\mathbf{g}\mathbf{N}}$ is the concentration in g N/tonne of soil organic nitrogen in the first 10 mm above soil (g N/metric ton soil), sed is the amount of sediment produced on the given day (metric tons), area is the HRU area (ha), and ${\mathsf{\epsilon }}_{\mathbf{N}\mathbf{s}\mathbf{e}\mathbf{d}}$ is the nitrogen enrichment rate.

The nitrogen enrichment rate is estimated according to Menzel [51]. Thus, the equation 15 was used to calculate the nitrogen enrichment rate:

${\mathsf{\epsilon }}_{\mathbf{N}\mathbf{s}\mathbf{e}\mathbf{d}}=0.78.\left({\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{s}\mathbf{e}\mathbf{d},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{q}}\right)$(Eq. 15),

where ${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{s}\mathbf{e}\mathbf{d},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{q}}$ is the in runoff sediment concentration (Mg sed/m³ H₂O).

The concentration of sediment in surface runoff is calculated according to equation 16 :

${\mathbf{C}\mathbf{o}\mathbf{n}\mathbf{c}}_{\mathbf{s}\mathbf{e}\mathbf{d},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{q}}=\frac{\mathbf{s}\mathbf{e}\mathbf{d}}{{10.\mathbf{a}\mathbf{r}\mathbf{e}\mathbf{a}}_{\mathbf{h}\mathbf{r}\mathbf{u}}.{\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}}$(Eq. 16),

where sed where sed is the amount of sediment on a given day (metric tons), ${\mathbf{a}\mathbf{r}\mathbf{e}\mathbf{a}}_{\mathbf{h}\mathbf{r}\mathbf{u}}$ is the HRU area (ha), and ${\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the amount of runoff on a given day (mm H₂O).

The movement of soluble phosphorus in the soil is by diffusion. Due to its low mobility, transfer by surface runoff only concerns soluble phosphorus stored in the top 10 millimeters of soil [50]. The amount of soluble P transferred is calculated from Equation 17:

${\mathbf{P}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}=\frac{{\mathbf{P}}_{\mathbf{s}\mathbf{o}\mathbf{l}\mathbf{u}\mathbf{b}\mathbf{l}\mathbf{e},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}.{\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}}{{\mathsf{\rho }}_{\mathbf{b}}.{\mathbf{d}\mathbf{e}\mathbf{p}\mathbf{t}\mathbf{h}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}.{\mathbf{k}}_{\mathbf{d},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}}$ (Eq. 17),

where ${\mathbf{P}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the quantity of soluble phosphorus retained during runoff, expressed in Kg P/ha, ${\mathbf{P}}_{\mathbf{s}\mathbf{o}\mathbf{l}\mathbf{u}\mathbf{b}\mathbf{l}\mathbf{e},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the quantity of soluble phosphorus in solution in the first 10 mm (Kg P/ha), ${\mathbf{Q}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the quantity of surface water runoff on the day in question (mm H₂0), ${\mathsf{\rho }}_{\mathbf{b}}$ is the soil density over the first 10 mm expressed in Mg/m³, ${\mathbf{d}\mathbf{e}\mathbf{p}\mathbf{t}\mathbf{h}}_{\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the depth of the surface layer (10 mm), and ${\mathbf{k}}_{\mathbf{d},\mathbf{s}\mathbf{u}\mathbf{r}\mathbf{f}}$ is the partition coefficient of phosphorus in the soil (m³/Mg).

Surface runoff to the main channel allows the organic and mineral phosphorus attached to soil particles transfer. All the nutrients parameters were integrated in SWAT model and nutrients flows were simulated running one time only the model with stream flow calibrated parameters.

3. Results

3.1. Streamflow parameter global sensitivity

The parameter global sensitivity analysis has showed that GW_REVAP, GWQMN, RCHRG_DP, ESCO and CN2 are the most sensitive parameters to the flow (Table 2). It is noted that the most sensitive among these parameters mentioned above are those related to the groundwater flow (GW_REVAP, GWQMN and RCHRG_DP). This sensitivity analysis showed a possible connection between groundwater and surface water within the watershed.

3.2. Streamflow calibration and uncertainty

The analysis of the results showed that the SWAT model is suitable for flow evaluation in the Lobo watershed at Nibéhibé. Indeed, R², NSE and PBIAS had values of 0.63, 0.62 and -8.1 respectively. The values of P_factor and R_factor were 0.48 and 0.52 (Table 3).

The graph in Figure 3 shows the results of the water flow calibration. SWAT model predictions generally follow the trend of observed flow. However, we note a shift in the flow peaks at the end of the simulation period (1981-1994).

From January 1, 1981 to January 1, 1990, the model has globally reproduced the observed fow correctly. But, from January 1, 1990 to December 31, 1994, the opposite is observed. Indeed, the SWAT model could not reproduce the peak flows and overestimated the base flow.

The uncertainties on the predictions are highlighted by P-factor and R-factor. The P-factor and R-factor values obtained were 0.48 and 0.52, respectively. This situation is explained by the fact that the SWAT model could not faithfully restore the peakflow from January 01, 83 to January 01, 85.

3.3. Nutrients fluxes

3.3.1. Nutrients requirements for crops

After integrating the different crops and fertilizers used in Lobo watershed, the SWAT model estimated the average nutrient requirement (N and P) of the main crops grown in the watershed (Table 4). We note that the average needs observed by the CNRA is 41 Kg/ Ha for nitrogen and 28 Kg/Ha for phosphorus. One notices an overestimation of N of + 15% and an underestimation of P of -24% compared to CNRA results.

3.3.2. Mineral nitrogen and soluble phosphorus transferred per sub-basin

The mineral nitrogen and soluble phosphorus spatial distribution on the basis of fertilizers used by crop type is presented in Figure 4A and Figure 4B. In Figure 4A, the variations in mineral nitrogen vary from 0 to 0.049 Kg/Ha. The highest mineral nitrogen fluxes are observed overall in the northern part of the watershed. The sub-basins concerned are 1, 2, 3, 4 and 5. Sub-basin 3 presents the largest contribution (0.049 Kg/Ha) of the entire Lobo catchment area taken in Nibéehibé. The south and southwest parts also show high flows but a little less compared to the north part. These quantities of nitrogen per hectare vary from 0 to 0.037 Kg/Ha. The sub-basins with large quantities are sub-basins 28, 29 (containing the Lobo water reservoir), 30 and 31. It is generally observed that the mineral nitrogen inputs over the entire watershed are weak. Figure 4B shows the spatial variation in soluble phosphorus levels by sub-watershed. The soluble phosphorus flux varies from 0 to 0.31 Kg/Ha. The higher values are observed in the central, eastern and western parts of the basin with rates of up to 0.31 Kg/Ha. It is also noted that the flows of mineral nitrogen are lower than those of phosphorus. This difference stems from the NPK fertilizer formula used.

3.3.3. Organic Nitrogen and organic phosphorus transferred per sub-basin

The simulation of organic nitrogen transport by sub-basin revealed inputs varying from 0.2 to 5 Kg/Ha (Figure 4C). One notes that the sub-basins in the southern part, except sub-basin 16, present a greater contribution of organic nitrogen and organic phosphorus than those of sub-basins located in the central and northern areas. The sub-basins with higher rates of organic nitrogen and organic phosphorus are 12, 13, 14, 24 and 29 with values reaching 3.7 to 5 Kg/Ha. As for organic phosphorus, it is shown in Figure 4D. Unlike the distribution of organic nitrogen, we observe higher amounts only in the northern part of the watershed.

The sub-basins with high inputs are 1, 2, 3, 4, 5, 6 and 17 with rates varying between 0.8 and 1.3 Kg/Ha. The results of this study show that the contributions of nutrients of exogenous and organic origin remain higher than those of mineral origin.

3.3.4. Nitrates and soluble phosphorus concentrations in streams

The monthly mean nitrate concentrations in mg/L were also estimated by the SWAT model. The predictions showed values ranging from 14 to 30 mg/L (Figure 5A). The result showed that the left bank and downstream part of the Lobo river is more vulnerable to an excessive intake of nitrates. This excessive inflow is extended to the outlet of the watershed (Nibéhibé) through sub-basin 29 where the Lobo water reservoir is located. Sub-basin 15 is the most vulnerable with concentrations up to 30 mg/L. Figure 5B presents the monthly average concentrations of soluble phosphorus in mg/L. The values vary from 0.1 to 1.8 mg/L (Figure 5B). The rivers in sub-watersheds 30 and 31 have the highest soluble phosphorus concentrations (1.3-1.8 mg/L). Sub-basins 21, 26, 28 and 29 also recorded high concentrations (0.9-1.3 mg/L). One noted that both for nitrates and soluble phosphorus, the sub-basin 29 containing the Lobo water reservoir recorded a significant contribution. A sampling campaign carried out made it possible to measure nitrate and phosphorus concentrations in surface water in the watershed. The measurement points have been represented on the maps of Figure 5A and Figure 5B.

These measurements revealed nitrate concentrations ranging from 7.4 to 74 mg/L and phosphorus ranging from 2.3 to 66.6 mg/L. In the absence of historical data that can be used to calibrate nutrient inputs (N and P), these measurement points were superimposed on the spatial distribution map in order to verify the robustness of the SWAT model predictions. It appears that high concentration points (dark blue for nitrates and red for phosphorus) coincide with areas of streams where predictions show a high value (orange and red). The points of high concentrations may remain the same over the years. However, the amount of nitrates and phosphorus accumulated can increase significantly over the years.

4. Discussion

The observed flow data used for this study are mostly available for the period 1981-1994. It is assumed that there is no major land use and rainfall trend over this period. In ungauged hydrological catchments where there is a scarcity of historical data, this hypothesis makes it possible to simulate flow and nutrient dynamics [52]. The study revealed that flow parameters such as GW_REVAP, GWQMN, RCHRG_DP, ESCO and CN2 impact the water flow the most. The same conclusion was made by Anoh et al. [53] in the Taabo watershed bordering that of Lobo. Along the same lines, Dakhlalla et al. [54] showed through a study on the evaluation of flow, nutrients and sediments transfer parameters sensitivity and uncertainty that hydrological parameters such as GWQMN, ESCO and CN2 are sensitive to flow. ESCO and CN2 are also parameters considered to significantly influence streamflow in western Mississippi [55]. Moreover, CN2, being a parameter very sensitive to runoff, a decrease in precipitation impacts surface runoff which will modify the values of CN2 [56,57,58] in [59]. However, a parameter can be identified as not sensitive. But, this does not really mean that it does not influence the flow in the watershed [60].

The values of R², NSE and PBIAS of 0.63, 0.62 and -8.1 respectively, highlighted the capacity of the SWAT model according to the evaluation criteria proposed by Moriasi et al. [35]. Upon the Taabo watershed, similar results were obtained by Anoh et al. [61]. Indeed, the average NSE obtained by these authors was 0.7 during the calibration period. A PBIAS value of -8.1 obtained in this study showed a very good performance of the model [35]. For flow, any P-factor > 0.7 and R-factor < 1.5 are considered acceptable. Therefore, the uncertainties on the flow predictions are considered acceptable [62]. P-factor has a value of 0.48 which is therefore less good since the wish is for this value to be close to 1 [29]. On the other hand, R-factor of 0.52 is acceptable. A model is said to be perfect when R-factor is close to 0. Very often, these values are very difficult to reach. This is why some authors recommend a balance between these two parameters [29]. In our case, P-factor is 0.48 and R-factor 0.52. There is a certain balance between these parameters. The prediction uncertainties are therefore acceptable. Despite these shortcomings of the model, the results obtained in this study constitute a solid basis for hydrological modeling in the region with a view to the fight against the eutrophication of surface water facing the population of this region. Compared with the results of the [63], the predictions of SWAT show an underestimation of the N rates of -9.9% and of -5.55% for P rates. The results show low inorganic nitrogen levels. These low inorganic nitrogen levels in the Lobo watershed have also been highlighted by the work of [3,4]. This could be explained by the fact that the recommendations of agricultural advisers to farmers regarding chemical fertilizers are respected. Then, nature, including crops, is able to digest the applied fertilizers, and only a small amount remains in the environment. The flows of mineral nitrogen are lower than those of phosphorus. This difference stems from the NPK fertilizer type used. Indeed, surveys carried out with ANADER (National Agency for Rural Development Support) and CIDT (Ivorian Textile Development Company) revealed that all agricultural activities practiced in the watershed use NPK fertilizer. In fact, NPK-10-18-18 fertilizer is used equally well for Coffee, Cocoa and cotton. But coffee and cocoa are the main industrial crops in this region of Cote d'Ivoire. This fertilizer formula contains 18% phosphorus against 10% nitrogen in 1 kg of this fertilizer, hence the high soluble phosphorus levels relative to nitrogen. Leakage of nitrogen and phosphorus by the SWAT model takes place outside the root zone, i.e., primarily by infiltration for nitrates and runoff for phosphorus, although the model estimates that a fraction of nitrogen migrates by leaching in runoff [64]. The variability of nitrogen and phosphorus leaks between sub-basins is conditioned by the vulnerability of soils to runoff and infiltration. Thus, soils favorable to runoff generate higher transfers of phosphorus and reduce those of nitrogen by reducing infiltration and promoting denitrification [65]. In general, the levels of mineral nitrogen are low compared to the levels of phosphorus. This is because nitrogen is almost always one of the deficient nutrients in cultivated soils with fragile organic status and particularly in ferralitic and ferruginous soils [66]. The contributions of nutrients of exogenous and organic origin remain higher than those of mineral origin. These results are in accordance with the work of [3] who showed that the fertilizer used by farmers throughout the Lobo watershed was manure, mainly used in the Séguéla area located to the north of the basin. According to Maïga et al. (2001) [3], the nutrient losses are of natural origin. Thus, natural fertilizers could provide 10 to 50 Kg/Ha for nitrogen (N) and 0.15 to 0.75 Kg/Ha for phosphorus (P) per hectare in the Lobo watershed. SWAT model simulations gave values lower than measured values in the field. However, the results have proved that the points of high concentrations fit well in the area of high nutrient concentrations predicted by the SWAT model. But, the amount of nitrates and phosphorus accumulated can increase significantly over the years. This explains the concentrations measured in the field, which reached as high as 74 mg/L for nitrates and 66.6 mg/L for phosphorus, the level above water quality criteria set by the World Health Organization (WHO) which stipulates that water with more than 50 mg/L of nitrates is considered polluted. In addition, most unpolluted freshwater contains between 10 and 50 μg/L [67] in [68]. This shows the state of pollution (eutrophication) of the Lobo water reservoir located in sub-basin 29, which is used for drinking water supply to Daloa city and around. These results agree with those of [4] who highlighted the eutrophication phenomenon of the Lobo reservoir. These authors have strived to examine the causes and consequences of such a phenomenon. In the same vein, [69] showed that the Lobo undergoes very advanced permanent eutrophication. The study carried out by [3] highlighted phosphorus values that can reach 0.066 mg/l. The OECD (Organization for Economic Cooperation and Development) threshold values between 0.035 and 0.1 mg/l classify Lobo reservoir in the category of hypereutrophic lakes [70]. The eutrophication phenomenon results from a primary overproduction in the water body, that is, an intense photosynthetic activity causing algal and plant development (autotrophic organisms), more significant than the consumption capacity of heterotrophic organisms (microorganisms, fish and crustaceans). There is then an imbalance between the trophic levels of the ecosystem. So, excess plant matter clutters the water body [70]. When the water body is invaded by floating macrophytes, they cover the surface of the water and limit the penetration of sunlight. The dead plant material settles at the bottom of the basin which it fills. This explains the black color of tap water mentioned by [5].

5. Conclusions

Considering the criteria for evaluating the flow setting, we note that the SWAT model is capable of simulating transfers in the Lobo watershed at Nibéhibé. The evaluation criteria R2, NSE and PBIAS had values of 0.63, 0.62 and -8.1 respectively. P-factor P and R-factor gave 0.48 and 0.52. On the basis of these values obtained, it appears that the SWAT model is indicated for the simulation of flows and nutrient transfers in tropical agricultural watersheds, in particular that of the Lobo.. The estimation of the average quantities of N and P applied to the crops gave 47.24 Kg/Ha of N and 21.25 Kg/Ha of P. The flows of inorganic N and soluble P vary from 0 to 0.049 Kg/Ha for inorganic N and from 0 to 0.31 Kg/Ha for soluble phosphorus. The northern and southern parts of the basin are more vulnerable to a large input of inorganic N when the central part remains vulnerable to the soluble phosphorus transfer. Simulations relating to the organic N transfer revealed values ranging from 0.2 to 5 Kg/Ha, the southern part showing the highest fluxes. As for the flow of organic phosphorus, they vary from 0.3 to 1.3 Kg/Ha. Unlike organic N, the greatest amounts are observed in the north. By comparing the flows of mineral (inorganic) nutrients with those of organic nutrients, one can say that the eutrophication observed in the Lobo water reservoir (sub-basin 29) would mainly have an exogenous origin (domestic, etc.) and organic. Predictions of nitrates concentrations and soluble phosphorus in the stream made using the SWAT model revealed concentrations of nitrates between 14 and 30 mg/L and soluble phosphorus between 0.1 and 1.8 mg/L. Field measurement points with high concentrations usually fall in areas of the stream with high concentrations simulated by SWAT. Apart from the limitations of this study related to the quality and quantity of the input data, it allowed us to gain a spatial understanding of the evolution and distribution of nutrients in the watershed. It also made it possible to identify the major inflow areas that would further contribute to eutrophication in the Lobo reservoir. These results constitute a solid basis which can serve as a basis in the fight against the phenomenon of eutrophication of water bodies in the watersheds of Côte d'Ivoire.