Forecasting water quality using seasonal ARIMA model by integrating in-situ measurements and remote sensing techniques in Krishnagiri reservoir, India

In recent years hydrologic and water quality mathematical models are developed and used by hydrologists for planning, managing and predicting water quality effects of changed conditions (Ziemińska-Stolarska & Skrzypski 2012; Loucks & Beek 2017; Liu 2018; Olowe 2018; Mbuh et al. 2019; Najafzadeh et al. 2021). Spatio-temporal change in water quality parameters serves to be an indicator in many studies (Abdollahi et al. 2017; Saturday et al. 2021) for environmental assessment through its trophic nature. Araujo et al. (2011) studied Funil Reservoir and stated that wet periods tend to bring in higher nutrient loadings, which can fuel algal blooms and degrade water. There is a direct linkage between nutrient loading and eutrophication (Yan et al. 2021). Nutrients such as nitrogen and phosphorus cling to sediments entering the reservoir feed the algae and proliferate its growth (Siriwardhana et al. 2019). Krishnagiri reservoir in Tamil Nadu, India, a hypereutrophic reservoir, was subjected to investigations on nutrient loads (Arunbabu et al. 2014), with the conclusion a systematic monitoring programme is essential to manage the reservoir water quality. The conventional approach to mapping the reservoir water quality is labour intensive and time-consuming. Also, it may not meet the desired standard of accuracy as it considers limited sample points and limited duration of monitoring (Kuha et al. 2020). With the advent of space technology, there are many satellites (IKONOS – Sawaya et al. 2003; IRS-LISS III – Mahato et al. 2004; LANDSAT – Wang et al. 2006; Georgas et al. 2009; Akbar et al. 2010; Wu et al. 2010; Modis Aqua – Ortega et al. 2010; Somvanshi et al. 2012; HICO- Keith et al. 2014) available on earth's orbit that can be used for water quality monitoring purposes (Gholizadeh et al. 2016).

Earlier in the Krishnagiri reservoir, a research work (Elangovan & Murali 2020) used a combination of bands 1, 3, 6, 7 of LANDSAT 8 OLI to develop a regression model for chlorophyll-a (chl-a; R²=0.635) with water quality parameters tested in the lab incorporating a time delay in testing with sparse sample points. The present study is also carried out in Krishnagiri reservoir, to improve the accuracy of chl-a estimation, lab testing of samples was replaced with the in-situ measurement using sophisticated chlorophyll sonde and multiparameter probes. This will eliminate the time delay in testing which eventually improve the performance of the regression equations for the optically active parameters such as chl-a, total dissolved solids (TDS) and electrical conductivity (EC). Further, the secondary goal of the study was to forecast the chl-a and EC of the reservoir to develop management strategies. A review of major studies in this area confirmed that the validated regressive model can be a useful tool in reservoir water quality management studies with a suitable forecast.

There are many time-series tools such as seasonal auto-regressive integrated moving average (SARIMA) least-square SVM (LSSVM) (Kaytez et al. 2015), multiple linear regression (Panklib et al. 2015), ARIMA-ANFIS (Barak & Sadegh 2016), ANFIS (Nokar et al. 2017), first-order fuzzy time-series (Tay et al. 2018) and ARIMA (Jain et al. 2018) and software available to find a fit for non-stationary data and develop a forecast. But most of the time-series data that remain non-stationary with a trend and seasonal variability such as water quality prediction (Katimon et al. 2018; Xu et al. 2019; Dastorani et al. 2020), dengue occurrence (Martinez et al. 2011), foreign tourism at the airport (Rusyana et al. 2016), precipitation (Geetha & Nasira 2016), temperature forecast studies, electrical consumption (Sim et al. 2019) and climate studies (Dimri et al. 2020), used SARIMA.

Historical data (2014–2020) of chl-a, and EC in Krishnagiri reservoir generated from the regression model exhibits seasonal variability and the SARIMA model was chosen for the forecast. Accordingly, the present study focused on integrating satellite data with time series to forecast the quality of inland water bodies (reservoir). Therefore the present work is the first attempt in Krishnagiri reservoir to predict the water quality parameters such as chl-a and EC for future years. The specific objectives of the study are (i) to develop a regression model integrating LANDSAT 8 OLI imageries with the in-situ measured optically active parameters such as chl-a, TDS and EC (ii) to validate the regression model with measured data and (iii) to forecast the water quality parameters (chl-a and EC) of the reservoir using the SARIMA model for the development of a better management strategy.

Krishnagiri dam (Figure 1) was constructed across the Ponnaiyar River in 1957 to sustain agriculture through irrigation. It also serves as a place for fish culture. It is located between 77 °41′56.4″ and 78 °16′51.6″ East longitude and 12 °23′45.6″ and 13 °2′27.6″ North latitude with a designed capacity of 68.2 MCM (Mohanakrishnan 1988). The dam receives water from a catchment area of 5,365 km² predominately has agricultural activities. The bathymetry survey conducted in 2012 showed that the reservoir has lost as much as 52% of its capacity since 1957 due to sedimentation as a result of soil erosion from the catchment area (Arunbabu et al. 2014). The inflow into the reservoir is generally maximum between October and November. It has two main canals, namely the Left Main Canal (LMC) and the Right Main Canal (RMC), running parallel to the Ponnaiyar River and irrigates an area of 3,642 hectares. The length of the dam is 1,000 m with eight spillways of span 12.2 m each, three river sluices and two canal sluices. The water spread area of the reservoir is 12.32 km² at full reservoir level (FRL). A hot climate prevails in the study region with a minimum and maximum temperature range of 22–24 °C and 34–37 °C respectively. The average annual rainfall in this region is about 800 mm.

Available evidence supports (Arunbabu et al. 2014; Elangovan & Murali 2020) a higher incidence of fish kills in the Krishnagiri Reservoir Project (KRP) attributed to increased cyanobacterial blooms resulting from decreased dissolved oxygen concentrations in the reservoir. The KRP reservoir was subjected to several investigations (Arunbabu et al. 2014; Elangovan & Murali 2020; Saha et al. 2021) in the aspects of capacity loss and nutrient loadings but there has been no study on the forecast of water quality.

A small motorised boat was deployed to conduct the field measurements using a multi-parameter probe (Aquaread) and a field chlorophyll sonde (Exo) accomplished for logging the concentration of parameters such as TDS, salinity, EC, chl-a, temperature and dissolved oxygen. Both the sensors recorded the geographical locations and the depth of measurements. All the parameters were measured in-situ at a point below 30 cm from the free surface of the water to determine the spatial variation in the reservoir. Sampling was carried out during October, November and December of 2019 and February 2020. The spatial and temporal distribution of chl-a was mapped in the GIS platform (ArcGIS 10.1 version).

LANDSAT 8 OLI imageries of spatial resolution 30 m × 30 m and seven spectral bands on row 143, path 51 captured the study reservoir once every 16 days. The images are first taken in units of absolute radiance, using 32-bit floating-point computations, which are then converted into 16-bit integer values in level 1 products, as represented by digital numbers in the image. These level 1 imageries were downloaded from the USGS website (https://earthexplorer.usgs.gov/). A total of 8 images between October 2019 and December 2019 and February 2020 were collected. Of this total, only the images obtained for October 2019, November 2019 and February 2020 could be used for the present study since the other imageries were not fully cloud-free. Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) in ENVI 5.3 package (ENVI 2009) was used to correct the imageries for atmospheric errors and produce the atmospherically corrected reflectance image.

NSE ranges between −∞ and 1.0 (1 inclusive), with NSE = 1 being the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas values <0 indicate that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance (Moriasi et al. 2007).

A regressed model with good accuracy is used for extracting historical data of chl-a and EC from reflectance band values of LANDSAT 8 OLI from 2014 to 2021. Not all months are free from cloud cover but one or two months in every quarter were left behind with minimal cloud cover. LANDSAT 8 imageries with 20% cloud cover were downloaded and processed with atmospheric correction. Sample points over the Krishnagiri reservoir for every imagery were established as shapefile in ARC GIS and band reflectance values were extracted. These reflectance values of bands as independent variables are substituted in developed empirical equations to estimate dependent variables chl-a and EC. The mean value of all parameter concentrations extracted from each sample point in a month was taken as the monthly average concentration of that parameter. Similarly, the average concentration of two or three months for every quarter was taken as a quarterly concentration of that parameter. Time series analysis and prediction can be carried out for a fixed time interval. For the present study, a constant time interval of three months with periodicity 4 (four quarters every year) is taken for analysis.

The IBM SPSS Statistics package offers time series analysis with SARIMA model using Box- Jenkins method (Sim et al. 2019). Auto correlation factor (ACF) and partial auto correlation factor (PACF) correlogram were used for selecting the model (auto-regressive (p) or moving average (q) or ARMA model (p, q) (Table 2) and seasonal model seasonal auto-regressive (P) or seasonal moving average (Q) or SARMA model (P, Q) (Table 3).

ACF and PACF correlograms developed with historical time series data for the parameters chl-a and EC were checked for their stationarity. If the ACF and PACF have large values (positive) that decrease very slowly with time, this indicates that the integrating factor (d) is greater than zero, i.e., differencing should be done, whereas if the ACF and PACF have lower values (negative) indicating the integrating factor (d = 0). Non-stationary series have an ACF that remains significant for half a dozen or more lags, rather than quickly declining to 0. Differencing such a series until it is stationary is required before identifying the process (IBM SPSS Forecasting 20 User's Guide).

Autoregressive processes have an exponentially declining ACF and spikes in the first one or more lags of the PACF. The number of spikes indicates the order of the autoregression. Moving average processes have spikes in the first one or more lags of the ACF and an exponentially declining PACF. The number of spikes indicates the order of the moving average. Mixed (ARMA) processes typically show exponential declines in both the ACF and the PACF. Seasonal processes show these patterns at the seasonal lags (the multiples of the seasonal period).

A suitable model for chl-a and EC was identified among all possible ARIMA (p, d, q) and SARIMA (P, D, Q) combinations through parameter estimates and Bayesian information criteria (BIC) values. A model with a low BIC value suits the best fit (Wali et al. 2017). ACF and PACF residues of chl-a and EC were used for model validation (Sim et al. 2019).

The forecasted results of chl-a and EC can be used to suggest suitable management practices to control the water quality degradation of the reservoir.

The water level in the dam was 12.79 m (Figure 2) at its highest in the year with a water spread of 10.96 km². Thus the 15 sampling locations (Figure 3) in October were well spread across the entire water spread. The release of water in the irrigation canals decreased the water level of the dam in November, December and February 2020. The lowest water level (7.68 m) was during February sampling (only 5 locations) with a water spread of 1.344 km². As there was a reduced inflow to the dam from the catchment, the water level in the dam decreased from October 2019 to February 2020, eventually decreasing the water spread area and hence the sampling sites and the numbers varied during each measurement.

All the measured data were then imported into the ArcGIS platform to map the spatial distribution of chl-a in the reservoir. During the period from October 2019 to February 2020, the chl-a concentration ranged between 14 mg/m³ and 150 mg/m³. There exists a spatial difference in chl-a concentration in all sampling months (Figure 4). There is a higher concentration in the inflow region, medium concentration in the middle and lower concentration near the dam structure during October (20.93 mg/m³–43.75 mg/m³), whereas November (14 mg/m³–52 mg/m³) showed an increased concentration at inflow point, lesser in the middle area and uniformly distributed near the structure. A similar trend is seen during December (19.26 mg/m³–39.80 mg/m³) with uniform distribution of chl-a concentration from middle reach towards the structure and greater concentration at inflow points. The north-western part of the reservoir, especially the river mouths (inflow region), showed higher chl-a concentration, especially in the monsoon periods. The reservoir receives its peak inflow from the catchment in October, which brings sediment load to the reservoir as a result of soil erosion from the catchment. The residual fertilizers on topsoil from the farming fields reach the reservoir in addition to the urban waste load from the nearby city (Saha et al. 2021). The release of nutrients and organic matter associated with fine sediments has deleterious effects on water quality and specifically can decrease oxygen levels of the reservoir. The chl-a concentration peaked in February 2020 up to 150 mg/m³ when the reservoir recorded the lowest water level of 7.63 m due to a significant reduction of inflow.

It is observed that the concentrations of chl-a were found to be decreasing from October to November and started to increase from December (Figure 5). Inflow and outflow volume from the reservoir has a direct relation to the change in concentration of chl-a (Elangovan & Murali 2020). Figure 5 shows the variation of chl-a with changes in inflow and outflow. Inflow volume into the reservoir showed that the increased inflow to the reservoir in October induced self-dilution of the chl-a concentration and hence there is a decreased chl-a concentration in November. Similarly, the decrease in inflow volume during November increased the concentration of chl-a in December.

TSI provides the level of eutrophication based on the chl-a concentration. The TSI value ranges from 62.93 in October 2019 to 80.38 in February 2020 (Table 4). The maximum value is recorded in February 2020 when the reservoir was shallow with no inflow and outflow. Also in the Krishnagiri reservoir for October and November, the TSI value was found to be at its minimum because of increased inflow to the reservoir. In the subsequent months water from the reservoir is released for irrigation and also there is a decrease in inflow into the reservoir. With decreased inflow, algae proliferation starts, thus increasing the TSI value from December 2019.

The in-situ test was performed using NAQUA READ multi-parameter probe and YSI chlorophyll sonde in Krishnagiri reservoir exactly on the date of the LANDSAT 8 satellite pass on 16th October (2019), 17th November (2019), 19th December (2019), 3rd February (2020), and measured the concentration of optically active water quality parameters chl-a, TDS, EC with a count of 20 sample locations every month. But points with cloud cover are not considered for analysis. Cloud-free sample points taken for analysis numbered six for October, nine for December, and 13 for February (2020) (Table 5). As the LANDSAT 8 image of December 2019 showed maximum cloud cover, it is not considered for regression analysis.

The correlation matrix for the chl-a, TDS and EC showed a significant correlation with the Coastal, Red, Green, NIR, SWIR 1, SWIR2, indicating that changes in the chl-a, TDS concentration and EC in the reservoir water would result in significant changes in the reflectance of significantly correlated bands (Table 6). Thus, it is possible to map the changes in chl-a, TDS concentration and EC with these bands.

Significant correlation exhibited with change in concentration of optically active water quality parameters and change in reflectance in the Coastal, Red, Green, NIR and SWIR 1 bands, indicating the quantity of energy reflected from the water surface containing algal matter in these bands will exhibit significant changes, compared to other bands. Accordingly, multiple regression equations for estimating the chl-a, TDS and EC were developed with the reflectance in the Coastal, Red, Green, SWIR 1, NIR and SWIR 2 bands (as independent variables) and chl-a, TDS concentration and EC (as dependent variables). Two equations for each parameter were developed as a result of the regression analysis (Table 7), using the combination of different bands.

Among the regressed model with different band combinations, one with greater R² and minimum standard error and RMSE were used for selecting the best model. From different regression operations, it is clear that bands Green, Red, NIR and SWIR1 are common in all models with greater R² (0.812 for chl-a, 0.972 for TDS, 0.980 for EC) and their band combination alone is enough to estimate the optically active water quality parameters such as chl-a, TDS, EC of Krishnagiri reservoir. The previous model developed (Elangovan & Murali 2020) for chl-a with the combination of OLI sensor bands Coastal, Green, NIR, SWIR1, SWIR2 for the Krishnagiri reservoir had a coefficient of determination of 0.635. The present study showed an improved coefficient of determination of 0.761 with the same band combinations eliminating the lag time of sample testing by in-situ measurements. Thus it is proved the lag time of sample tests had a significant influence on the estimation of chl-a. Even though the samples are stored and transported to the lab for testing, in-situ measurements at the time of satellite pass would improve the model accuracy. Further, the inclusion of the Red band ((Table 6) 75% correlation with chl-a) in the combination also improved the model accuracy (R² = 0.812).

The NSE for the regression model for the three parameters chl-a (0.812), TDS (0.830) and EC (0.960) showed better performance (Figure 6) as NSE is close to 1.

The regressed model with good accuracy with greater R² and minimum standard error (Table 7) was used for extracting historical data of chlorophyll-a and EC from reflectance band values of LANDSAT 8 OLI from 2014 to 2021 (Table 8).

Analysis of time series invariably involves the evaluation of trends and seasonality in the data. Trends are the long-term increase or decrease in the time series, whereas seasonality refers to the variations in the data at regular short intervals such as weekly, monthly, biyearly, quarterly, etc. (Wang et al. 2013). Time series data of chl-a shown in Figure 7 has no trend in its flow but has some seasonality variation. As the data is already detrended differencing parameter ‘d’ has no value.

Also, from the ACF and PACF correlograms (Figures 8 and 9, respectively) it is evident that the data is detrended as the ACF and PACF do not have large values (positive) that decrease very slowly with time (Dimri et al. 2020). autoregressive (AR; p), moving average (q), seasonal AR (P) and seasonal MA (Q) can be estimated with the help of ACF and PACF plot with first-order seasonal differencing D = 1. From the correlogram with a small significance in seasonal periodic lag 4 of ACF (Figure 8), a peak significance in periodic lag 4 of PACF (Figure 9) and no significant peak in lag 1, 2, 3 of ACF and PACF it is inferred that the parameter ‘p’ and ‘q’ has no value while P = 1 also with a less chance of Q = 1. So our possible models will be S(P,D,Q) of S(1,1,1) and S(1,1,0) otherwise ARIMA (1,1,1)₄ and ARIMA (1,1,0)₄

The SPSS platform is used to run the identified models (1,1,1)₄ and (1,1,0)₄ and the best fit is selected based on BIC and parameter significance test (Tables 9 and 10, respectively). BIC is a criterion for model selection among the finite set of models where the model with the lowest BIC is preferred (Rahman & Hasan 2017; Wali et al. 2017). The SARIMA model generated a forecasted time series plot of chl-a and enlarged forecasted plot of chl-a for Krishnagiri reservoir is shown in Figures 10 and 11 respectively.

From the model statistics (Table 9), it is inferred that model (1,1,0)₄ is found to be the best fit with a low BIC value and model parameter significance <0.05 (Table 9), accepting the hypothesis that the model is significant, and it can be accepted. Sim et al. (2019) identified a significant model using the parameter significance test. The residual plots of ACF and PACF shown in Figure 12 can be used to check the adequacy of the model. It shows a random variation from the origin zero (0), the points below and above are all uneven, hence the model fitted is adequate.

Excess of salt increases the osmotic pressure of the soil solution, a situation that can result in a physiological drought condition. Thus, even though the soil in the field appears to have plenty of moisture, the plants will wilt. This occurs because the plant roots are unable to take up soil water due to their high osmotic potential (Zaman et al. 2018).

The total soluble salts (TSS) content of irrigation water is measured either by determining its EC, reported as micro siemens per centimetre (μS cm⁻¹), or by determining the actual salt content in parts per million (ppm). According to USSL Staff (1954) salinity class C4 is not suitable for irrigation (Table 11). Time series data of EC (Table 8) showed non-stationarity by the Dickey-Fuller test (p = 0.614 > 0.05), and Phillips-Perron test (p = 0.80 > 0.05). Jalil & Rao (2019) defined a better non-parametric test for stationarity in the Phillips-Perron test.

It is inferred that there are unit roots in data from both tests and the p-value is significant for non-stationarity (greater than 0.05). Non-stationarity in data may be due to seasonality. From the time series of EC (Figure 13) and PhD thesis of Karunakaran (2004), it is evident that the EC value has an increasing trend over two decades with a mean value of EC as 0.707 milli siemens per centimetre in 2002, thus incorporating difference factor (d, D) with a time series plot of EC. Possible model parameters estimated from the ACF and PACF plot with d = 1 and D = 1 model shown in Figures 14 and 15 are

• (1) (2,1,0) and (1,1,0)

• (2) (2,1,1) and (1,1,0)

• (3) (2,1,0) and (0,1,0)

From the model statistics (Table 12) it is inferred that EC model_3 (2,1,0) (0,1,0)₄ is found to be the best fit with a low normalised BIC value of −3.521 and model parameter significance <0.05 accepting the hypothesis that the model is significant (Table 13), and it can be accepted.

The residual plots of ACF and PACF shown in Figure 16 can be used to check the adequacy of the model. It shows a random variation from the origin zero (0); the points below and above are all uneven, hence the model fitted is adequate.

From the forecast (Table 14), it is found that by the start of 2032 first quarter (Q1) Electrical conductivity is nearing 2.25 mS/cm. Irrigation water above 2.25 mS/cm is classified under C₄ (Table 11). It takes almost a decade from 2022 for the Krishnagiri reservoir to degrade in its salinity from class C3 (high salinity) to class C4 (very high salinity) thus making it unsuitable for irrigation purposes. SARIMA model generated Forecasted time series plot of EC is shown in Figure 17.

The present study mapped the spatial and temporal distribution of chl-a in Krishnagiri Reservoir. TSI derived with chl-a provided the level of eutrophication. It ranged from 62.93 in.October 2019 to 80.38 in.February 2020 indicating the higher nourishment of nutrients. The spatial distribution of chl-a indicated higher concentrations at the inflow regions in the monsoon periods whereas the chl-a distribution was uniform in February when the reservoir recorded its lowest water level in that season. Regression results of LANDSAT 8 OLI images and optically active water quality parameters (chl-a, TDS, EC) can provide adequate information on the chl-a, TDS, EC concentration variations for Krishnagiri reservoir.

• (i)

  The present study concluded that the time difference between the sample collection in the site and the lab testing should be an influencing factor for the model accuracy. This is proved with the in-situ measurement of chl-a showing a higher coefficient of determination (R² = 0.761) than the sample tested at lab (R² = 0.635) (Elangovan & Murali 2020).

• (ii)

  Inclusion of the Red band, which showed 75% correlation with chl-a, had significantly improved the model accuracy from 0.635 to 0.812.

• (iii)

  The estimate of the chl-a, TDS and EC concentration based on the developed regression model exhibited a high R² value of 0.812, 0.945 and 0.960 respectively.

• (iv)

  Forecast results of chl-a with the SAR model (1,1,0)₄ showed a clear indication that the reservoir will remain hypereutrophic for the next two years (up to 2023) (any value of chl-a >22 mg/m³ is hypereutrophic).

• (v)

  The EC forecast with significant model parameters (2,1,0) (0,1,0)₄ alerts the decision-makers that the salinity of the Krishnagiri reservoir will be shifting from class C₃ (high salinity) to class C₄ (very high salinity) in the next decade.

It is evident from the study that integrating remote sensing data with in-situ measurements will significantly help in the monitoring and management of reservoirs. The developed regression model is an effective tool in extracting historical water quality data where historical data is not available and the SARIMA model can be used to predict seasonal water quality transformation over years. The results of the study can be used by the watershed managers to create the necessary framework for watershed conservation programs to reduce the sediment and nutrient loads to the reservoir.

The authors are grateful to the Centre For Water Resources, Anna University for providing wet chemistry lab facilities for performing various analyses.

The authors declare that they have no competing interests regarding the publication of this article.

Abdul Wahid performed the conceptualization methodology, writing, collected data and performed the analyses with the help and supervision of Arunbabu. Both authors contributed to the writing and review of the manuscript.

All relevant data are included in the paper or its Supplementary Information.