Investigation of scarce input data augmentation for modelling nitrogenous compounds in South African rivers

Fresh water is a scarce resource in South Africa, with an average annual daily rainfall of 490 mm (WWF-SA 2016) being around half the global average. A recent local annual climate report shows that precipitation patterns have remained unchanged in South Africa (SAWS 2020). There has therefore not been any increase in rainfall while the demand for water increased steadily. In addition to water scarcity, South African freshwater resources are subjected to strain through pollution from underperforming water treatment facility effluents. Eutrophication, a result of excess phosphoric and nitrogenous nutrients in a river system, has been previously highlighted by van Ginkel (2011) and Harding (2015) as a serious problem in the country. Therefore, effective water resource management is crucial to managing the quantity and quality of available in-demand levels.

Water resource management solutions require a reliable representation of the state of water in the regions of the country where water is monitored and managed. This can be achieved through the use of established sampling networks that record the quantity and quality of available water for various uses. Numerical models can be further used to complement actual measurements to achieve a similar objective.

Numerical models can describe processes (such as geo-hydrology, climate properties, sources, sinks) recorded in the field to provide a more reliable and accurate representation. This is possible when adequate measured boundary condition data as well as data for calibration and validation are available. When data are scarce, augmentation techniques can be applied to adjust the recorded information for the model input, which is necessary for a more accurate simulation of the environment. Data availability can therefore limit the choice of reasonable modelling approaches (including items of structure, complexity, and spatial and temporal resolution) that can be applied (Slaughter et al. 2017).

Spatial and temporal resolution are examples of necessary requirements for the model study that seeks to simulate a dynamic process such as nitrification in rivers. This study focuses on the impact of augmented model input data generated by applying simple (basic interpolation) and advanced (artificial neural networks (ANNs)) augmentation methods applicable to modelling nitrogenous compounds in the river system.

In instances where the input data are of unsatisfactory resolution; data augmentation techniques may be applied for gap filling and disaggregation to meet the model requirements (Baffaut et al. 2015). Examples of recent successful applications of data augmentation through interpolation to water quality modelling can be obtained in the literature (see Yang et al. 2020 and Kim et al. 2021). Blöschl & Sivapalan (1995) details how upscaling and downscaling of data can be used to align the scales of available data with model requirements. In this context, upscaling refers to the transfer of information from a small scale to a large scale and downscaling refers to transferring information to a small scale. Downscaling consists of two steps (disaggregation and singling out) to transfer information to a smaller scale (Blöschl & Sivapalan 1995). Input hydrodynamic data may be disaggregated from monthly to daily flows to accommodate the variation in daily concentration during water quality modelling. This was demonstrated when flow duration curves and mathematical relationships were applied to flow data for the WQSAM model to study water quality (see Slaughter (2017)). Disaggregation in hydrology is one of the steps used for downscaling hydrology data to meet the scales required for meeting the modelling objective.

Time series for hydrological processes can also be generated using either deterministic (basic interpolation methods) or stochastic (ANNs) models. Koutsoyiannis et al. (2008) detail the differences in the application of stochastic models against deterministic models. The important finding in this study is that good stochastic models are those that are linked with an understanding of the natural behaviours of the system. Furthermore, deterministic models such as the analogue model can be a good simplistic analytical tool, however, good prediction from this model does not necessarily represent consistency with natural processes. Since this study focuses on nitrogenous compounds which can be closely linked to hydrogen cycle processes through river stream flow, some stochastic models may not be ideal for generating time-series analysis depending on whether autocorrelation structure is for short range dependence or long-range dependence (Dimitriadis et al. 2021).

There is a variety of gap-filling methods in the literature that have been applied successfully in water quality modelling studies. Table 1 lists studies on this topic.

In this study, different augmentation methods applied towards modelling a nitrogenous compound in the river system were investigated. This was done to inform augmentation method choice when dealing with scarce data. First, the impact of interpolation method choice as well as whether the level of gaps in the input data impacts the outcome of each augmentation method differently was investigated. Second, an advanced interpolation method (machine learning) was applied towards gap filling while exploring the best method for partition training and validation given limited data. Finally, the output of the models was compared for scenarios where different levels of artificial gaps are introduced in the input data.

The outcome of this study provides guidance towards augmenting data for water quality modelling under scarce data conditions, which are similar to real-world data availability in the study area. Improvement of model fit to the measured data is beyond the scope of this paper, however there is evidence that the selection of the appropriate machine learning method can translate to better model fit (see Rozos et al. 2022). This study is limited to observing the significance of the differences between the augmentation method output.

Two scenarios were investigated. First, the effect of using deterministic interpolation for data augmentation. Second, the effect of using deterministic and machine learning interpolation to infill weekly gaps at locations where maxima occur in the ammonia input was also investigated.

Two data operations (gap filling and temporal disaggregation) were required to eliminate the gaps in the original dataset and to upscale the temporal resolution of the input data from weekly to daily frequency. This is to meet the model requirement for dynamic processes according to the recommendations by Moriasi et al. (2012).

The chosen basic interpolation methods (linear, quadratic, cubic, spline (first and second order), polynomial, piecewise polynomial, and derivatives) in the Python Pandas library (The pandas development team 2021) were applied to upscale data to a daily resolution. A detailed discussion of these interpolation methods can be found in Virtanen et al. (2020). The scope of this article only covers the application of each method; whereas method algorithm details can be found in the references listed in Table 2.

Each of these interpolation methods was applied to generate an augmented dataset with a daily frequency; as recommended by Baffaut et al. (2015) for modelling the dynamic processes that affect nitrogenous compounds in river systems. Upsampling in this study refers to when the frequency of the samples is increased such as from weekly to daily. For this to be done, the generation of a time series is required. Resampling, in this case, is required to transform the available irregular frequency of data to a regular frequency and to increase the number of samples to create more data on which the neural network can be trained. It is important to note that resampling methods have the limitation of destroying the long-range dependence that appeared in the data.

Two different simulation models were included in this study.

In Equation (1), V represents the reactor volume, c is parameter concentration in the reactor, W(t) represents the lumped loading, t is the time, k is the reaction rate constant, A_s is the cross-sectional area, and v is the flow velocity.

These equations were solved using Python's fourth-order Runge-Kutta solver for numerical integration because it was simple to apply to the system of differential equations in this study. Ammonia concentration was computed on the selected checkpoints in the river reach.

Water Quality Analysis Simulation Programme (WASP) software was additionally used to model the river reach. WASP (Wool et al. 2020) is an open-source dynamic compartment-modelling programme for aquatic systems, including both the water column and the underlying benthos. It was developed and distributed by the Environment Protection Agency in the United States (Wool et al. 2020). It allows the user to investigate 1-, 2-, and 3-dimensional systems as well as a variety of pollutant types and processes such as eutrophication. This programme was selected because it is a recognised software that is capable of modelling nitrogenous compounds in a river system. Details about the development of this programme and capabilities can be obtained from the work of Wool et al. (2020).

Measured on-site data sets were available, but with varying spatial consistency, which required cleaning where data values were absent in addition to gap filling for missing data. Each data sample (upstream and downstream from the WWTP) showed varying information. Table 3 provides a list of raw data features.

Hydrodynamic data measurements on the reach were not available. Estimates of flow rates were derived using measurements from a nearby station at a reach with similar features (catchment size and climate) as recommended by Daggupati et al. (2015) when dealing with data scarcity. The available hydrodynamic data covered a full year (between 1 October 2017 and 1 October 2018). The reader is reminded here that the focus of this research was not to test model fit, but rather to investigate differences in model outcomes when using different water quality data augmentation methods. This data application was therefore accepted as a realistic representation of hydrodynamic data for the system.

As mentioned previously, the study area consisted of a single river reach with one WWTP discharging into the stream. The upstream boundary of the study consisted of observed water quality data with weekly temporal resolution for the period 2012–2020. The water quality parameters relevant to this study were monitored nitrogenous compounds (ammonia and nitrates). The data sets upstream of the WWTP were used as a boundary condition for the river models; the downstream observed data sets were used for output comparisons. The augmentation study focused on the impact that the boundary condition data had on the model outputs as observed at the downstream boundary. The focus was to discern whether applying different augmentation techniques would yield significant differences in the model outputs.

The input ammonia concentration data were divided into four categories:

• No gaps – the raw measured upstream boundary data as measured.

• Low gaps – raw data with 10% random artificial gaps.

• Medium gaps – raw data with 25% random artificial gaps.

• High gaps- raw data with 50% random artificial gaps.

Each data set was subjected to interpolation for the data gaps and upsampled to daily concentration data through applying the linear, quadratic, derivatives, cubic, piecewise polynomial, 1st order spline and 0th order spline interpolation methods. The augmented data were used as input to the Basic Model (CSTRs in series) and the WASP model in turn to simulate a 5.9 km long single reach river system for nitrogenous compounds with ammonia selected as a proxy parameter to nitrogenous compounds to represent the changes brought about through the nitrification process. The simulation models were driven by flowrates as estimated using transference as explained above. The output of each model for each input data set as generated through the use of the interpolation methods was compared to determine statistically significant differences.

As previously stated, two scenarios are investigated. First, we investigate the effect of using deterministic interpolation for data augmentation. Second, we investigate the effect of using deterministic and machine learning interpolation to infill weekly gaps at locations where maxima occur in the ammonia input.

The log transform is needed here to prevent the minority of very large ammonia values to influence our model training too significantly: by applying a log transform action, the small and large ammonia values become comparable. Normalising the data generally speeds up learning and leads to faster convergence.

To analyse the dependence (or correlation) of consecutive measurements on each other, an autocorrelation plot was made of the upsampled data. Autocorrelation plots of the data are shown in Figure 4.

Figure 4 shows a significant correlation for a lag of 3 days or less. It should be noted that this can only be used as a guideline when selecting the number of timesteps to use to predict the next sequence of timesteps; since it only provides information on the linear- and not the non-linear relationship of the measurements. The pre-processing stage was then concluded by dividing the data into a training- and a testing set for which 70% of the upsampled data were used for training the neural network and 30% was used for testing.

The validation plot did not increase towards the training line, which is an indication that the model was not over-trained.

To visualise the effectiveness of the model to predict gaps that contain local maxima, artificial 7-day gaps were created within the test set by selecting positions that contain local maxima.

In this case, the machine learning method seems to predict high peaks compared to the rest of the methods (linear, Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)) except for the polynomial method at the local maxima.

This part of the investigation was done to determine the effect of time-series data disaggregation by interpolation on the boundary of a river simulation method. The design layout involved simulating the selected model for each of the interpolation methods. The output of the models for each method was compared, within a model, and then between the models.

To evaluate the comparisons between the models, statistical techniques were applied to determine if the difference between the simulated output of the interpolation model was statistically significant. If the difference is significant, this will indicate that the choice of an interpolation method has an impact on the outcome of a simulation model. Therefore, the choice of the augmentation model should be considered a crucial variable when a model is developed for a scarce data system.

The opposite of this outcome would suggest that for this dataset, the choice of an interpolation method has no significant impact on the output of the models. This would imply that when setting up a model with inadequate data, the choice of interpolation method from the list discussed in this study, does not make a significant difference to the model output.

To quantify this, two statistical methods were applied to determine the differences caused by the model boundary interpolation methods. The differences between the interpolated boundary conditions were determined using quantitative statistics and inferential statistics. For this study, the T-test and the ANOVA test were applied to the model results.

A T-test is a type of inferential statistic used to determine if there is a significant difference between the means of two groups (model output from different interpolation methods), while the ANOVA test does the same for more than two groups. The null hypothesis for both tests assumes that there is a significant difference between the groups, which is an assumption made about the groups. The T-test and ANOVA test used t values and F values, respectively, to determine whether the null hypotheses pass or fail. A probability value (p-value) is also calculated for each test that suggests whether the hypothesis should be accepted (chosen p < 0.05).

The results of the model simulation for the Basic Model and the WASP model were analysed using the T-test and the ANOVA test. This was done to evaluate how each method compares for both models. This was done to determine the differences between all the methods in each model. This was followed by the T-test for methods output comparison across models (Basic Model output against WASP model output) to evaluate the differences between the models.

As was found for the Basic Model results, the WASP model results showed seemingly small variations. The 25% IQR value for this model yielded 0 meaning there was no discernible variability between the methods.

The model outputs from the Basic Model and WASP produced similar trends for each interpolation method for all input data with varying artificial gaps. The results from the low gaps dataset (10%) and medium gaps dataset (25%) in Figures 9 and 10 infilled the ammonia concentration with a higher value than the actual measured values at the original dataset. This is because the random gaps were introduced at a critical peak towards the end of the time series. The dataset with medium gaps (25%) (Figure 10) resulted in the largest deviation between interpolation method outputs for both models. This indicates that the location of gaps may be of importance as opposed to the level of gap sizes. Furthermore, the results indicate that in the case that random gaps were introduced at critical peaks, the choice of interpolation method may be important. This notion is supported when looking at the case in which high gaps (50%) were introduced, which indicated few differences between the interpolation methods, contrary to intuition.

Further statistical investigation of the simulation results was done to determine and to quantify differences between methods if they exist. The ANOVA test was applied to the Basic Model output data to compare the different outputs when applying the different interpolation methods, with the varying levels of random gaps on the input datasets. Table 4 provides a summary of the test results.

Table 4 lists a small F-value as well as the p-value (PR) > 0.05 for all datasets except for the one with 25% artificial gaps (PR = 0.0001), which suggests that for these datasets the differences between the method outputs were not statistically significant. The contrary result for the case where 25% gaps were included (statistically significant differences in results) may have been due to the introduction of most of the gaps at the peak values in the original dataset. This once more indicated that the location of data gaps is of concern.

Additionally, a post hoc (Tukey) test was performed to compare the outputs from the application of the different interpolation methods only on the dataset without artificial gaps. This was to observe if there were any significant differences between the methods themselves. The results of the test are listed in Table 5.

The p-values for each individual comparison are shown to be above 0.05, which further confirms the rejection of the null hypothesis. These results indicate that there were no statistically significant differences between model outputs when applying the different interpolation methods.

The ANOVA test was similarly applied to the WASP model output data. The results are listed in Table 6.

The results from this test also support the rejection of the null hypotheses (PR > 0.05) except, once more, for the input dataset with 25% artificial gaps. Therefore, as was found for the Basic Model, differences in the output values when applying the different interpolation methods were not statistically significant for the WASP model. This is further confirmed by the post hoc test result listed in Table 7, for the input dataset without gaps.

Method-by-method comparison yielded a rounded-up p-value > 0.05 for each case, solidifying the rejection of the null hypothesis. This corresponds to the results from the outputs obtained on the data set generated by the Basic Model.

The study results indicated that there was no statistically significant difference between the outcomes of each interpolation method applied to the full dataset, however, the introduction of random artificial gaps resulted in significant differences in outcomes between interpolation methods for the case where 25% of the gaps were introduced to the original dataset. Machine learning approaches produced reasonably accurate results. However, upsampling was necessary to obtain the recommended minimum data size, required for learning.

The following was concluded:

• Application of the different interpolation methods applied to input water quality data did not produce statistically significantly different augmented datasets with low gaps.

• Increasing the gaps in the original data sets did not always yield greater differences between augmented datasets for each method.

• The locations of the artificial gaps created statistically significant differences between the augmented datasets for each interpolation method when compared to the effect of high gaps.

• The selected machine learning methods to infill real and artificial gaps were successful in upsampling the original dataset and the dataset with artificial gaps.

• There was no significant difference between the simulated output of WASP and the Basic Model.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.