ORGANICS: A QGIS Plugin for Simulating One-Dimensional Transport of Dissolved Substances in Surface Water

Abstract

Surface water in streams and rivers is a valuable resource and pollution events, if not tackled in time, may have dramatic impacts on aquatic ecosystems. As such, in order to prepare pollution prevention plans and measures or to set-up timely remedial options, especially in the early stages of pollution incidents, simulation tools are of great help for authorities, with specific reference to environmental protection agencies and river basin authorities. In this paper, we present the development and testing of the ORGANICS plugin embedded in QGIS. The plugin is a first attempt to embed surface water solute transport modelling into GIS for the simulation of the concentration of a dissolved substance (for example an organic compound) in surface water bodies including advection dispersion and degradation. This tool is based on the analytical solution of the popular advection/dispersion equation describing the transport of contaminants in surface water. By providing as input data the concentration measured at the entry point of a watercourse (inlet boundary condition) and the average speed of the surface water, the model simulates the concentration of a substance at a certain distance from the entry point, along the profile of the watercourse. The tool is first tested on a synthetic case. Then data on the concentration of the pharmaceutical carbamazepine monitored at the inlet and outlet of a vegetated channel, in a single day, are used to validate the tool in a real environment. The ORGANICS plugin aims at popularizing the use of simple modelling tools within a GIS framework, and it provides GIS experts with the ability to perform approximate, but fast, simulations of the evolution of pollutants concentration in surface water bodies.

1. Introduction

Surface water in streams and rivers is a valuable resource, and pollution events, if not tackled in time, may have dramatic impacts on aquatic ecosystems [1,2,3,4]. As such, in order to prepare pollution prevention plans or to set up timely remedial options, especially in the early moment of pollution incidents, simulation tools are of great help for authorities, with specific reference to environmental protection agencies and river basin authorities. Modelling of solute transport in surface water is then a valuable and common option [5,6,7,8,9], especially by means of analytical solutions simplifying system description and reducing complexity [10,11,12], providing fast, even approximate, answers, particularly in cases of insufficient data availability for the implementation of more complex numerical models.

In order to advance the use of modelling tools and to support the digitalization of the technical sector these tools must be user-friendly, and built around free and open-source codes. An open code may guarantee the reproducibility and the reliability of the analyses performed [13,14] and their early deployment and impact [15,16].

Geographic Information Systems (GISs) are worldwide mainstreaming tools and methodologies for storing, managing/analyzing and visualizing large temporal, spatial and non-spatial datasets (for both geometric and alphanumerical data). Hence, they are widely applied to support environmental modelling [17,18]. The possibility of including modelling tools in Geographic Information Systems (i.e., via plugins), and the existing community of users and developers in constant growth, is potentially the most relevant strength of such tools [18,19,20]. For thirty years GIS and water modelling codes/applications have been successfully integrated. Vieux [21] presented an application of the GIS, ARC/INFO and the finite element solution to the kinematic wave equations to process the spatially variable terrain in a small watershed using a Triangular Irregular Network for the solution of overland flow. Oliveira et al. [22] presented ArcGIS-SWAT, a geodata model and GIS interface for the Soil and Water Assessment Tool (SWAT; [23]). Becker and Jiang [24] developed a computationally efficient method for predicting contaminant mass flux to a specified boundary, carrying out the method in a GIS and taking full advantage of widely available digital hydrologic data. Akbar et al. [25] showed a GIS-based modelling system called ArcPRZM-3 for the spatial modelling of pesticide leaching potential from soil towards groundwater. Rossetto et al. [26] integrated a suite of groundwater modelling tools in the gvSIG GIS application [27]. Oliveira and Martins [28] developed an application for the preliminary characterisation of the river boundary condition for a MODFLOW [29] finite difference groundwater flow numerical model. Bittner et al. [30] developed the LuKARS GIS-based model for simulating the hydrological effects of land use changes on karst systems.

Among GIS desktop applications, QGIS [31] is probably the most popular free geospatial software. Rosas-Chavoya et al. [20] conducted a bibliometric analysis on the acceptance of this application on documents published in Scopus from 2005 to 2020, considering 931 manuscripts. They observed a favorable trend in the acceptance of QGIS across the world and the development of large collaborative networks.

Several plugins have been developed for QGIS within the hydrological or aquatic domain. In particular, Nielsen et al. [32] developed the Water Ecosystems Tool, a workflow implemented (as a plugin) in QGIS, for the application and evaluation of aquatic ecosystem models. Ellsäßer et al. [33] developed the QWaterModel as an easy-to-use tool to make evapotranspiration predictions available to broader audiences. The QWaterModel is a QGIS plugin compatible with all versions of QGIS3. Dile et al. [34] developed an open-source user interface for the SWAT [23], QSWAT, using various functionalities of the QGIS application. Rossetto et al. [18] presented the FREEWAT plugin for managing the groundwater resource, including tools for the management of hydrochemical data [35] and nitrate leaching assessment [36].

In this paper, we present the development and testing of the ORGANICS plugin as a first attempt to embed surface water solute transport modelling into GIS. This tool allows users to simulate the concentration of a dissolved substance (for example an organic compound) in surface water bodies by applying an analytical solution of the advection dispersion equation, which includes also a first-order degradation term. By providing as input the stepwise time-variant concentration measured at the entry point of a watercourse, along with the related average water velocity, the concentration of a substance at a certain distance from the entry point along the profile of the watercourse, is simulated. A sketch of the presented problem is shown in Figure 1.

After presenting the theoretical and modelling approach, we show an example application of the plugin (to be used as a tutorial), and then a real case study application to simulate carbamazepine concentration in a vegetated channel collecting poorly treated wastewater.

2. Materials and Methods

The model used in the ORGANICS plugin is based on the analytical solution of the popular advection/diffusion/decay equation, in one-dimensional form, taken from [10,37]: where:

$$\frac{\partial C}{\partial t}+\frac{\partial v_{x}C}{\partial x}=E\frac{{\partial }^{2}C}{\partial x^2}-kC$$

(1)

C: is the solute concentration expressed as mass per unit volume of water [M/L³],

v_x: longitudinal fluid flow velocity is the input velocity [L/T].

x: is the longitudinal coordinate [L],

t: is time [T],

E: is the longitudinal dispersion coefficient accounting for the combined effects of ionic or molecular diffusion and hydrodynamic dispersion [L²/T].

k: is a first-order decay rate [T⁻¹].

Using a constant concentration boundary condition at x₀ = 0 [L] (the inlet point of a surface water body reach) at the initial time (t₀ = 0) [T], for each x > 0 and t > 0 the Equation (1) may be reduced to the following analytical solution (2) [10,37].

$$C\left(x,t\right)=\frac{C_0}{2}\left\{\exp \left[\frac{Ux}{2E}\left(1-\Gamma \right)\right]\mathrm{erfc}\left(\frac{x-Ut\Gamma }{2\sqrt{Et}}\right)+\exp \left[\frac{Ux}{2E}\left(1+\Gamma \right)\right]\mathrm{erfc}\left(\frac{x+Ut\Gamma }{2\sqrt{Et}}\right)\right\}$$

(2)

$$\Gamma =\sqrt{1+4 \frac{kE}{U^2} }$$

(3)

and U is the velocity module [L/T].

This approach entails the following assumptions:

(a) the flow is one-dimensional and oriented according to the main direction of the flow in the surface water body;

(b) the concentration at the inlet remains constant for the specified simulated time interval (first kind boundary condition);

(c) chemical interactions between different dissolved substances are not considered, nor are reactive geochemical processes simulated with other components (i.e., riverbed matrix);

(d) the morphology of the bed of the surface water body does not affect the solution;

(e) no sorption processes or production terms are considered;

(f) at t₀ = 0 [T] the initial condition is C(_x,₀) = 0 [M/L³] along all the simulated domains.

The developed code uses Equation (2) to calculate the concentration value along the line input provided by the user to represent the selected surface water body. This input can be given as a linear vector layer (for instance the common ESRI. shp file). The code calculates the solution (concentration value) at nodes at homogeneous lengths from the inlet point as specified by the user, and at selected homogeneous time-steps also defined by the user.

2.1. Software Development

The plugin was developed in Python3 [38] language, using the Qt5 [39] graphics libraries and the QGIS [40] Application Program Interface. The solute transport code is embedded in the QGIS desktop application. This means a full integration at programming language level, with models using GIS data format included as full component of the host GIS application [41,42].

The plugin can be used as an add-on to the QGIS desktop application, version 3.4 or higher, once it has been installed on a PC. During the development phase, some specific Python libraries were used, but each of them is already included in the official QGIS desktop distribution. This choice allows the user to use the plugin without requiring further software updates.

2.2. Data Needs

In order to use the plugin, the user has to prepare a set of input files. These are:

• a *.csv file specifying the water average longitudinal velocity (U) and concentration values at the inlet of the watercourse (constant concentration boundary condition), and the time these data refer to. The file must comply with the template format defined for the plugin. In particular, the file must contain data relating to (at least) one dataset, specifying:

  -

  starting date and time, in YYYY-MM-DD HH: MM: SS format;

  -

  average flow velocity in the surface water body, in m · s⁻¹. This value will be used at all the node of the surface water body;

  -

  the concentration of the source at the inlet point.

When considering time-varying boundary condition (that is the concentration input changing with time), the user must specify for each different time all the above information on consecutive lines of the. csv file.

• an ESRI linear Shapefile representing the surface water body. The file may consist of one or more segments. The line must be digitized towards the flow direction. When more segments are used, the topology must be respected (all the lines must be connected).

An example of the required files to run a first test are provided in the template_files folder of the plugin itself as Supplementary material.

2.3. Model Implementation and Run

Once data are prepared in the form of the required files, the first operation consists in loading the Shapefile geographic layer of the line into the QGIS view.

By clicking on ORGANICS in the Plugin menu, the main window opens. This is divided in four sections: (1) Run; (2) Plot Results; (3) Help; (4) About. By entering the Run section, the user input the following data (Figure 2):

(a) the *.csv file. Upon the file selection, the drop-down menus will automatically update. Through these menus the user must select the name of the fields in the *.csv file corresponding to the required information;

(b) the linear shapefile representing the watercourse. At this step the user must specify the length of the homogeneous reaches at whose ends the concentration values will be calculated;

(c) the value of the first-order decay rate coefficient (s⁻¹);

(d) the value of the coefficient of longitudinal dispersion (m²/s);

(e) the length of the timestep (in minutes) at which the solution will be calculated over a time interval (in minutes from the start of the simulation) at each point of the reach;

(f) the name of the output file (*.shapefile) and the directory where the file will be saved. Should this field be left blank, the output will be saved as a temporary layer (memory layer) in QGIS.

Although time data are input both in seconds and in minutes, all the calculations are internally run in seconds, while results are provided in minutes or within hours.

Once the simulation is run, the code will: (i) divide the line in several nodes, according to the length specified by the user, and (ii) calculate the solution at each node at the end of a reach, at each time step as specified by the user. The result consists in a *.shapefile point layer in which the simulated concentration values at different times are saved at each point. This layer can be used to visualize the results using the tools in QGIS. For example, it can be themed with color scales depending on the level of the selected solute concentration. Animations to visualize the evolution of the concentration in the various points may be produced by applying the TimeManager plugin (plugin, downloadable from the QGIS PluginManager).

A *.csv file of the output will also be saved in the previously defined destination folder. This file can be used by the user to conduct further analyses externally to the plugin and/or from QGIS (i.e., using spreadsheet).

Graphs of the solutions may be drawn opening the Plot Results section, where a number of options for producing solution plots or further customizing the draw and to save it in image format (i.e., as *.png file) are provided (Figure 3).

The user must choose the output layer to process. Graph drawing can be performed at any time, even after the execution of the model, by selecting the output file from the drop-down menu. However, in order for the desired layer to appear in the ORGANICS menu, the layer must be loaded in the QGIS layer panel. Once the layer is selected, the drop-down menus below will automatically update.

The following options for creating the graphs are available:

• Select a position (distance in m from the entry point): this option will create a graph displaying the concentration trend in a point defined by the user at a certain distance from the starting point, as a function of time (Figure 4);
• Use the selected position on the layer: this option allows to view the same result as above, but in this case the position is provided by selecting, using the classic selection tools on the map, one or more points of the output layer (Figure 5);
• Select a time: this graph will display the concentration values, at a given simulated time, as a function of the distance from the inlet (Figure 6). The times available for selection correspond to the discretization obtained with the time step chosen in input by the user.

3. Results and Discussion

3.1. Model Validation

We successfully validated the model implemented within the tool by simulating the same case described by the analytical solution shown in [10]. The simulation results (Figure 7) are obtained using the following parameters:

U = 1.0 m/s

k = 0 s⁻¹ (no decay is simulated)

E = 5.0 m²/s

where:

C_(0,0) = C₀, and

C_(x,0) = 0 at x > 0

The solution (concentration value) is provided at 100, and 1000 m from the source.

3.2. Example Problem

In this section, three tests are presented in order to show the behavior of the calculated solution in different scenarios. The cases tested are:

(a) C₀ mass injection, constant over time;

(b) time-limited pulse C₀ mass injection;

(c) C₀ mass input, variable over time (multi-pulse input condition).

Concerning the cases with time-limited or time-varying input (cases b and c, respectively), we implemented a solution using the principle of superimposition of the individual analytical solutions at each stress period, as reported in [10,43,44]. The parameter values used in the example problem are presented in

Table 1. Parameters values used in the example problem.

Parameter	Value	Units
Total length of the reach	1200	m
Simulation length	50	m
Time step	10	min
Velocity (U)	0.1	m/s
First order decay rate (k)	0.00005	s−1
Longitudinal dispersion (E)	5.0	m2/s
Initial time	28 May 2018 00:00	dd/mm/yyyy hh:mm

. The .csv file and the .shp file are provided in the Supplementary Material.

3.2.1. C₀ Mass Injection, Constant over Time

In this test, we simulated the release of a mass with concentration C₀ = 100 ng/L constant over time. The solution is presented at the end of the simulation time, which is the time needed for the mass to reach the outlet.

Figure 8 displays the solution at the beginning of the reach (x = 0 m), at x = 500 m, and approximately at the end of the reach (x = 1100 m). In the middle of the reach (x = 500 m) the solution tends to an asymptotic value, which is less than 100 ng/L due to the effect of the simulated decay process.

3.2.2. Time-Limited Pulse C₀ Mass Injection

In this test we simulated a time-limited pulse C₀ = 100 ng/L mass injection for a 2 h duration. This time-limited pulse case, at constant concentration, can be implemented by defining in the *csv file an initial period (of known duration, with concentration C₀; first line in the *csv file) followed by a second period with zero concentration (second line in the *csv file). In this test, the second input is then two hours long with concentration set at C_2h = 0 ng/L. The solution is then displayed in Figure 9 at an infinite time (that is, the time needed for the dissolved substance to reach the outlet of the water course considered). The solution is presented at x = 500 m (Figure 9). At this distance, mass arrival is recorded after 30 min from the beginning of the simulation.

3.2.3. C₀ Mass Input, Variable over Time (Multi-Pulse Boundary Condition)

In this test, we simulated C₀ mass input, variable over time (multi-pulse input condition) according to data presented in

Table 2. Data used in the test presented in Section 3.2.3.

Date and Time	C0 (ng/L)
28 May 2018 0:00	0
28 May 2018 0:20	100
28 May 2018 2:00	50
28 May 2018 2:30	25
28 May 2018 3:30	75
28 May 2018 4:00	0

. The global solution works as the superposition of the several pulses, each one having a constant condition for a specified time interval. Superposing each “pulse-solution” makes the model able to consider time-dependent boundary conditions. Results are shown at x = 400 m, x = 800 m, and x = 1200 m from the inlet point for time step length of:

-

20 min (Figure 10);

-

10 min (Figure 11);

-

5 min (Figure 12);

from the beginning of the simulation, in order to present the impact of the different time discretization on the solution (Figure 10, Figure 11 and Figure 12).

3.3. Case Study Application

The ORGANICS plugin was then applied to compute the concentration of the pharmaceutical compound carbamazepine at a reach of a vegetated channel receiving poorly treated wastewater, in low-flow conditions, in the Pisa municipality (Tuscany, Italy, Figure 13). Carbamazepine (CBZ) is an anticonvulsant or anti-epileptic drug commonly found in poorly treated wastewater and consequently in surface- and ground-water [45,46]. Removal efficiency in secondary wastewater treatment is typically less than 10% for CBZ [47]. Composite samples (2 volumes of 0.5 L every 30 min representative for one hour) were collected approximatively every two hours during an experiment run on 28 May 2018 at the inlet (point PSM_w = 0 m) and the outlet (point PSM_z = 420 m) of the channel reach (Figure 13) in low-flow conditions. Analytical determinations were performed following the method described in [48]. Mean longitudinal flow velocities were measured by means of an acoustic digital current meter (OTT Messtechnik GmbH, Kempten; Germany). Data for CBZ and mean longitudinal flow velocities are presented in

Table 3. Carbamazepine concentrations at the inlet (PSMw) and outlet (PSMz) points of the vegetated channel on 28 May 2018, and simulated results.

Time	Inlet (PSMw) (ng/L)	Outlet (PSMz)(ng/L)	Flow Velocity (m/s)	Simulated Value (PSMz_sim) (ng/L)
07:20	123	-	0.025	-
09:50	181	-	0.026	-
11:00	-	105	-	112
12:20	162	-	0.026	-
13:10	-	112	-	116
14:20	162	-	0.021	-
15:40	-	115	-	119
16:50	150	-	0.029	-
18:00	-	125	-	122

.

A multi-pulse boundary condition was prepared by exploiting data from the five monitoring points. Figure 14 shows the simulated carbamazepine concentration at the outlet point, PSM_z, 420 m from the inlet point. The Simulated value (PSM_{z_sim}) column in Table 3 presents the simulated value against the measured ones (Outlet (PSM_z) column). In Figure 13 CBZ simulated concentrations are themed with color scales depending on the concentration value.

The first-order decay rate and the longitudinal dispersion coefficients are relevant parameters in our analyses, and, more in general, in surface water quality modelling [49,50]. No experimental data were available for these parameters. As such, the model was calibrated to get the best fit (R² = 0.95) between simulated and measured concentrations with values of the longitudinal dispersion coefficient of 35 m²/s and decay rate equal to 3 · 10⁻⁵ s⁻¹. Good fit (R² > 0,9) was also obtained varying these two parameters within the range of 30 and 35 m²/s for longitudinal dispersivity and between 2.5 · 10⁻⁵ and 3 · 10⁻⁵ s⁻¹ for the decay rate coefficient. The values of the longitudinal dispersion coefficient are coherent with values found in [49,50] for similar open channels. The calibrated values of the decay rate are slightly higher than calculated values from half-life time data for CBZ reported in [51].

4. Conclusions

The developed open-source and free plugin allows simulating transport of dissolved substances in water courses following advection/dispersion, and degradation processes. The present formulation combining a simple analytical solution of the advection dispersion equation and GIS tools guarantees intuitive spatial data management. Authorities may also benefit from the ease of use of such tools in order to set in place pollution prevention measures. This solution, because few parameters are needed, could hence be applied to data-scarce environments. Furthermore, using this tool values for longitudinal dispersivity and first-order decay rate coefficients may be derived. In our case study application, we estimated the longitudinal dispersion coefficient and the decay rate coefficient to be 35 m²/s and 3 · 10⁻⁵ s⁻¹ for the pharmaceutical compound carbamazepine. Another potential application could be in the feasibility stage of the design of water-related green infrastructures for the improvement of water quality [46,52]. On the other hand, the increasing number of integrated geographical databases (including surface water bodies characteristics) along with the increasing availability of sensors gathering and distributing quasi real-time monitoring data (such as surface water heads) may allow for the combined use of monitoring and modelling to set up early warning systems to track pollution events [53,54,55]. To this aim, further research and pilot experimental sites are needed.

The open and free characteristics of the developed code guarantee reproducibility of the run analyses, and also use of the tool, one of the most important aspects of science, making free software an ideal framework for scientific work. In this context, the use of free software is consolidating in our societies in a gradual, but constant way [56]. Integration in the FREEWAT plugin [18] and inclusion within the list of the official plugin of QGIS will guarantee the dissemination and potentially the application of this research product.

In the present formulation, the ORGANICS tool does not allow users to simulate the transport of substances under conditions where flow may increase/decrease downstream not only by tributaries, but also by continuous groundwater drainage from the surrounding domain. Nor does it allow for the spatial or time variability of the degradation rate. The latter could be beneficial to differentiate, for example, the relative importance in time of biodegradation from photodegradation processes [50,57,58]. Future development may include integration of more complex analytical solutions, comprising also source terms.

The ORGANICS plugin is a first attempt to popularize the use of simple modelling tools. For more complex solutions and the inclusion of time-varying source/sink terms, a wide range of numerical tools exist [7,10]. We wish therefore to stress that the main element of this tool resides in its simplicity. Finally, the ORGANICS plugin provides GIS experts the ability to perform approximate, but fast simulations of the evolution of pollutants concentration in surface water bodies at selected targets.