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

**Rudy Rossetto 1,\*, Alberto Cisotto 2, Nico Dalla Libera 2, Andrea Braidot 2, Luca Sebastiani 1, Laura Ercoli <sup>1</sup> and Iacopo Borsi <sup>3</sup>**


**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.

**Keywords:** water pollution; solute transport modelling; Geographic Information System (GIS); pollution prevention plans; pharmaceuticals; carbamazepine; longitudinal dispersion coefficient; decay rate coefficient

#### **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–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–9], especially by means of analytical solutions simplifying system description and reducing complexity [10–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].

**Citation:** Rossetto, R.; Cisotto, A.; Dalla Libera, N.; Braidot, A.; Sebastiani, L.; Ercoli, L.; Borsi, I. ORGANICS: A QGIS Plugin for Simulating One-Dimensional Transport of Dissolved Substances in Surface Water. *Water* **2022**, *14*, 2850. https://doi.org/10.3390/w14182850

Academic Editor: Cristina Di Salvo

Received: 5 August 2022 Accepted: 8 September 2022 Published: 13 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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 opensource 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.

**Figure 1.** Schematic draw of contaminant movement along a surface water reach.

#### **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 \mathcal{C}}{\partial t} + \frac{\partial \upsilon\_{\mathcal{X}} \mathcal{C}}{\partial \mathbf{x}} = E \frac{\partial^2 \mathcal{C}}{\partial \mathbf{x}^2} - k \mathcal{C} \tag{1}$$

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

*vx*: 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 [L2/T].

*k*: is a first-order decay rate [T<sup>−</sup>1].

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

$$\mathcal{C}(\mathbf{x},t) = \frac{\mathbb{C}\_{0}}{2} \left\{ \exp\left[\frac{\mathsf{U}\mathsf{x}}{2E}(1-\varGamma)\right] \text{erfc}\left(\frac{\mathbf{x}-\mathsf{U}t\varGamma}{2\sqrt{Et}}\right) + \exp\left[\frac{\mathsf{U}\mathsf{x}}{2E}(1+\varGamma)\right] \text{erfc}\left(\frac{\mathbf{x}+\mathsf{U}t\varGamma}{2\sqrt{Et}}\right) \right\} \tag{2}$$

$$
\Gamma = \sqrt{1 + 4 \frac{kE}{\mathcal{U}^2}} \tag{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*<sup>0</sup> = 0 [T] the initial condition is *C*(*x*,0) = 0 [M/L3] 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:

	- starting date and time, in YYYY-MM-DD HH: MM: SS format;
	- average flow velocity in the surface water body, in m · <sup>s</sup>−1. 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<sup>−</sup>1);

(d) the value of the coefficient of longitudinal dispersion (m2/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.


**Figure 2.** The ORGANICS plugin main window.

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:



**Figure 3.** The ORGANICS plugin *Plot results* section.

**Figure 4.** Example graph: solution at a selected point as function of time.

**Figure 6.** Example graph: solution at a selected time as function of distance from the inlet point.

#### **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* =0s−<sup>1</sup> (no decay is simulated) *E* = 5.0 m2/s where: *C*(0,0) = *C*0, and *C*(*x*,0) = 0 at *x* > 0

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

**Figure 7.** Computed simplified analytical solution implemented in Organics in time at *x* = 100 m and *x* = 1000 m.

#### *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:


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. The .csv file and the .shp file are provided in the Supplementary Material.

**Table 1.** Parameters values used in the example problem.


#### 3.2.1. *C*<sup>0</sup> Mass Injection, Constant over Time

In this test, we simulated the release of a mass with concentration *C*<sup>0</sup> = 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.

**Figure 8.** *C*<sup>0</sup> mass injection, constant over time: concentration values at selected points of the reach.

#### 3.2.2. Time-Limited Pulse *C*<sup>0</sup> Mass Injection

In this test we simulated a time-limited pulse *C*<sup>0</sup> = 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*0; 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 *C2h* = 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.

**Figure 9.** Time-limited pulse *C*<sup>0</sup> mass injection: concentration values at 500 m from the inlet.

3.2.3. *C*<sup>0</sup> Mass Input, Variable over Time (Multi-Pulse Boundary Condition)

In this test, we simulated *C*<sup>0</sup> mass input, variable over time (multi-pulse input condition) according to data presented in Table 2. 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:


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

**Figure 10.** Simulated concentration at *x* = 400 m, x = 800 m, and *x* = 1200 m from the inlet point with 20 min time step length.

**Figure 11.** Simulated concentration at *x* = 400 m, *x* = 800 m, and *x* = 1200 m from the inlet point with 10 min time step length.

**Figure 12.** Simulated concentration at *x* = 400 m, *x* = 800 m, and *x* = 1200 m from the inlet point with 5 min time step length.



#### *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 *PSMw* = 0 m) and the outlet (point PSMz = 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.

**Figure 13.** Case study location, investigated channel reach, and monitored points.

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


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, PSMz, 420 m from the inlet point. The *Simulated value (PSMz\_sim)* column in Table 3 presents the simulated value against the measured ones (*Outlet (PSMz) column*). In Figure 13 CBZ simulated concentrations are themed with color scales depending on the concentration value.

**Figure 14.** Simulated concentration at *x* = 20 m, and *x* = 420 m from the inlet point with 10 min time step length.

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 (R2 = 0.95) between simulated and measured concentrations with values of the longitudinal dispersion coefficient of 35 m2/s and decay rate equal to <sup>3</sup> · <sup>10</sup>−<sup>5</sup> <sup>s</sup><sup>−</sup>1. Good fit (R2 > 0,9) was also obtained varying these two parameters within the range of 30 and 35 m2/s for longitudinal dispersivity and between 2.5 · <sup>10</sup>−<sup>5</sup> and <sup>3</sup> · <sup>10</sup>−<sup>5</sup> <sup>s</sup>−<sup>1</sup> 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 m2/s and 3 · <sup>10</sup>−<sup>5</sup> <sup>s</sup>−<sup>1</sup> 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–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.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/w14182850/s1. The Supplementary Material folder contains: folder organics: contains the plugin to be installed in QGIS v3.xx; folder example\_problems: contains the files used to run the Example problem described in Section 3.2 of the paper.

**Author Contributions:** Conceptualization, R.R. and I.B.; methodology, R.R., I.B., L.S., L.E., A.C. and N.D.L.; software, I.B. and R.R.; validation, R.R., N.D.L. and L.S. and L.E.; formal analysis, R.R. and L.E.; resources, A.B. and R.R.; data curation, R.R. and I.B.; writing—original draft preparation, R.R. and I.B.; writing—review and editing, R.R., A.C., N.D.L., L.S. and L.E.; visualization, R.R. and I.B.; funding acquisition, A.B. and R.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper presents the ORGANICS plugin. The development of the ORGANICS plugin was financed through the research agreement between Distretto Idrografico Alpi Orientali and Istituto di Scienze della Vita—Scuola Superiore Sant'Anna: *CONVENZIONE DI RICERCA per lo sviluppo di un modello con metodi numerici integrati in software per la gestione dei dati spaziali (GIS) che descriva i meccanismi di trasporto, abbattimento e circolazione in acque superficiali di fitofarmaci e composti azotati* (2019– 2021), co-funded within the GREVISLIN INTERREG Italia-Slovenia Project (Green infrastructures for the conservation and improvement of the condition of habitats and protected species along the rivers, CUP number: G76I18000080007) [https://www.ita-slo.eu/en/grevislin (accessed on 1 August 2022)].

**Data Availability Statement:** Data used in the section "Case study application" come from research run within the Italian-Israeli bilateral project PHARM-SWAP MED (co-funded by the Italian Ministero degli Affari Esteri e della Cooperazione Internazionale).

**Acknowledgments:** The authors wish to thank Sara Pasini and Matteo Bisaglia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

