Biomathematical Model for Water Quality Assessment: Macroinvertebrate Population Dynamics and Dissolved Oxygen

Abstract

Sustainable water management is important to ensure its availability for future generations. The study of water quality is fundamental for this purpose. Assessing the health of aquatic ecosystems through bioindicators has been shown to be reliable and inexpensive. The objective of this work was to evaluate water quality through a biomathematical model that involves environmental stress indicator organisms and their close relationship with dissolved oxygen. In this direction, a system of differential equations describing the population dynamics of aquatic macroinvertebrates under the influence of dissolved oxygen is proposed. The model is validated by its application in the Coyuca Lagoon, Mexico. Likewise, population changes over time were represented, which allowed us to deduce that the increase or decrease in the aeration/oxygenation rate significantly affects the population dynamics of the bioindicator organisms. In addition, to classify water quality, a one-to-one correspondence was established between water quality and the equilibrium points of the system of differential equations. The results obtained allow inferring that the proposed techniques are useful for the study of water quality since they can predict significant changes in the ecosystem and provide researchers and water managers with tools for decision making.

1. Introduction

One of the sustainable development goals for 2030 is to guarantee the availability of water and its sustainable management [1]. To achieve this goal, the health of aquatic ecosystems must be continuously studied and monitored to know their current and future conditions, ensure better quality, and prevent degradation. In this direction, the assessment of the quality of freshwater ecosystems is fundamental, and the development of prospective methods contributes to this purpose.

Some methods for assessing water quality make use of indicators that include biological, physicochemical, and geomorphological parameters (see [2,3,4,5,6], and references therein). In addition, it has been shown that such indicators directly or inversely relate to some parameters [7]. Theoretical representations of observable phenomena in different contexts have also been developed [8,9,10]. Beyond the advantages of these methods, their main limitation is that they only provide a temporal state of the ecosystem, but do not integrate changes in water quality over time. As a solution, differential mathematical models are used, which are more representative of the state of water quality over time.

Dissolved oxygen (DO) is a very important element that has been used in the evaluation of water quality in aquatic ecosystems [11,12]. It is also fundamental for life and as an indicator of ecosystem degradation [13,14,15,16]. Likewise, DO is key to changes in aquatic communities because some organisms respond negatively or positively to different levels in their concentrations [17,18]. For example, eutrophication causes a decrease in water quality [19] and an increase in hypoxia [20], which implies the absence or presence of environmental quality indicator organisms, particularly some aquatic macroinvertebrates (AMs) taxa [21,22]. Furthermore, Galic et al. [23] in their meta-analysis of the adverse impacts of hypoxia on aquatic invertebrates infer that reproductive and respiration rates decrease when DO availability decreases and those organisms rapidly reduce their reproductive capacity in the face of environmental stress.

AMs as biological indicators have proven to be reliable in the study of water quality, habitat disturbances, and restoration of aquatic ecosystems. Likewise, AMs are sensitive or tolerant to high organic pollution loads and increased or decreased DO [24,25,26,27]. Moreover, they have different ecological requirements and degrees of tolerance to natural and anthropogenic disturbances [28,29]. In addition, AMs are used to study the ecological health of natural areas [30] and the effectiveness of the use and application of artificial barriers, which fragment the landscape, as a method of reclamation [31]. Consequently, the application of biotic indices, qualitative and quantitative, using AMs has become popular worldwide (see [32] and references therein), even though such methods are retrospective, and do not help to determine or have a close idea of what may occur in the future. On the other hand, to contribute to the reliability of AMs applications, they are included in the construction of multimetric and topological indices [33,34], and in some prospective methods based on a set of environmental parameters.

Models for assessing water quality involve one or several environmental parameters, with different methods and theories [34,35,36,37,38,39]. These models evaluate water quality in a timely manner or pose possible scenarios over time. For example, the Streeter–Phelps model relates DO deficit and biochemical oxygen demand [40], through a system of two coupled differential equations. In addition, this model is taken as a starting point in a historical review of models for water quality assessment [41]. According to the context, the models are either very simple or complex depending on the number of equations, input data, or time in the sampling effort [16]. Recently, as reported by Costa [42], AQUATOX is one of the most applied models in the last five years, because it predicts the behavior of various pollutants and their effects on the ecosystem. AQUATOX within its set of mass balance equations contains differential equations to observe the change, with respect to time, in benthic biomass and dissolved oxygen concentration, independently [43].

The objective of this work is to study water quality in freshwater ecosystems through time and population fluctuations of environmental stress indicator organisms. In this direction, a differential mathematical model describing the population dynamics of AM indicators of water quality and their relationship with DO is proposed. To this end, Section 2 describes the state variables, parameters, and the system of four coupled nonlinear differential equations that represent the dynamics between AMs and DO; Section 3 shows the main qualitative theoretical results in the ecological context. Finally, Section 4 shows the numerical simulations. In addition, we validate the model by applying it to the water quality assessment of a real case: the Coyuca Lagoon, Mexico.

2. Model Conceptualization

Disturbances in freshwater ecosystems can be studied through the population dynamics of environmental stress indicator organisms, such as AMs, and since the population fluctuations of AMs have a direct or inverse relationship with DO concentration, in water quality assessment, and restoration of aquatic ecosystems, certain taxa can be classified according to their oxygen requirements and tolerance to organic pollution. Thus, we define the sets $A=\{a_1,\cdots ,a_n\}$, $B=\{b_1,\cdots ,b_m\}$, and $C=\{c_1,\cdots ,c_s\}$ as Class I, II, and III organisms, respectively, where each class corresponds to the AM families with higher/lower, regular/regular, and lower/higher oxygen requirement and tolerance to contamination. Some images, common names, and stress tolerance classifications (Classes) of AMs can be consulted in [44,45].

The cardinal of A, denoted by $\#A$, is the number of Class I AM families. Similarly, the cardinals for sets B and C are defined. The cardinal of $a_i\in A$, for $i=1,\cdots ,n$, denoted by $\#a_i$, is the abundance of the aquatic macroinvertebrate family $a_i$. The cardinal of $b_j\in B$ and $c_k\in C$, for $j=1,\cdots ,m$ and $k=1,\cdots ,s$, are defined analogously to $\#a_i$. Taxonomic richness is the sum of the cardinals of sets A, B, and C. Total abundance is the sum of the cardinals of $a_i$, $b_j$, and $c_k$.

The structure of aquatic macroinvertebrate communities is defined by their abundance, richness, and diversity [46,47,48]. Richness can often be assessed at higher taxonomic levels, e.g., genus, family, and order. In addition, richness measures reflect the diversity of the aquatic assemblage. On the other hand, relative abundance, defined as the ratio of a taxon to all others contained in an ecosystem, determines how rare, common, or dominant a taxon is. A population is defined as a group of organisms of the same species, each population has several characteristics, such as its population density or relative size in an available space. In this study, each population is a conglomerate of AM populations that have similar characteristics in the face of environmental impact, particularly in terms of their DO requirements.

The research on the evaluation of water quality, previously analyzed, has allowed us to infer the basic principles for the development of this research. Thus, the differential mathematical model proposed in this study considers the intraspecific competition between AM families and as a limited resource: the concentration of DO. In addition, the oxygenation/aeration of the water body, either natural or mechanical, reaches a maximum saturation level. Moreover, regarding the DO consumption and population growth of the different classes of AMs, the concentration of DO in the water is diminished by pollution associated with biotic components, particularly AMs, and the amount of energy absorbed by living organisms is greater than that which they use.

Based on the above, the state variables represent: $o\left(t\right)$ the concentration of DO in the water and $x_1\left(t\right),x_2\left(t\right)$, and $x_3\left(t\right)$, the sum of the cardinals $a_i,b_j$, and $c_k$, i.e., the sum of the abundances of the AM families of class I, II, and III, respectively. Note that, the state variables depend on t, but in denoting them, we will dispense with this. The flow diagram of the model is shown in Figure 1.

Mathematical Formulation

The mathematical representation of the above is described by the following system of non-linear differential equations:

(1)$\begin{array}{c}\hfill \begin{array}{cc}\hfill o^{\prime }& =r\left(1-{\displaystyle \frac{o^2}{m^2}}\right)-o\sum _{i=1}^3{\beta }_i{x}_i,\hfill \\ \hfill x_i^{\prime }& =r_i{x}_i\left(1-{\displaystyle \frac{m{x}_i}{k_{i}o}}\right)+{\gamma }_i(m-o)x_i,\hfill \end{array}\end{array}$

for $i=1,2,3$, in the positive invariant region

$\Omega =\{(o,x_1,x_2,x_3)\subset {\mathbb{R}}^4:0<o\le m,0\le x_i\le k_i,\mathrm{with}i=1,2,3\},$

and the parameter set $\omega =\{m,r,{\beta }_i,r_i,k_i,{\gamma }_i\in \mathbb{R}:m,r,{\beta }_i,r_i,k_i,{\gamma }_i>0,\mathrm{with}i=1,2,3\}$.

The differential equations associated with system (1) represent the rate of change over time in DO concentration and AM populations. The equations for the AM families make use of the logistic equation since this allows modeling phenomena with limited resources and populations of organisms that reach their maximum density, which is limited by the availability of resources. Moreover, its application in the study of populations has been widely extended [49,50]. All parameters are positive constants and their interpretation is as follows: r is the constant oxygenation/aeration rate of the water body, m is the maximum value to which the DO tends asymptotically in the absence of AMs, and the DO consumption for the different classes of AMs is given by $o{\sum }_{i=1}^3{\beta }_i{x}_i$. For the different classes of organisms: ${\beta }_i$ is the average respiration rate, $r_i$ is the average reproduction rate, $k_i$ is the average carrying capacity i.e., the average upper limit of population growth [51], ${\gamma }_i$ is the average rate of DO utilization, for $x_i$ with $i=1,2,3$. Furthermore, it is established that ${\gamma }_i<{\beta }_i$ because organisms consume some amount of energy, but do not fully utilize it. Note that if $o=m$, then $x_i$ asymptotically approaches $k_i$, for $i=1,2,3$. The parameter values or ranges, which will be used in the numerical simulations in Section 4, are shown in Table 1.

3. Qualitative Analysis

Although the system of differential Equation (1) has 23 equilibrium points (two without the presence of AMs, nine with the presence of only one class of AMs, nine with the presence of two classes of AMs, and three with the presence of the three classes of AMs), our study will focus only on those that are in the region of interest and that describe the absence or presence of different classes of AMs; more specifically, we focus on the equilibrium point of coexistence of the three AM classes for the information provided to us.

3.1. Absence, Presence, and Coexistence of AMS

Proposition 1. 

Given the system of the nonlinear differential Equation (1), there is $o^{*}\in (0,m)$, such that the following points of ${\mathbb{R}}^4$ are the equilibrium points of the system and lie in the region of interest Ω.

• (i) 

  AM-free equilibrium: $E_0=(m,0,0,0)$.

• (ii) 

  Presence equilibrium: $E_1=(o^{*},x_1^{*},0,0)$, $E_2=(o^{*},0,x_2^{*},0)$, $E_3=(o^{*},0,0,x_3^{*})$,

  $E_4=(o^{*},x_1^{*},x_2^{*},0)$, $E_5=(o^{*},x_1^{*},0,x_3^{*})$ and $E_6=(o^{*},0,x_2^{*},x_3^{*})$.

• (iii) 

  Coexistence equilibrium: $E_7=(o^{*},x_1^{*},x_2^{*},x_3^{*})$ with

  (2) $x_i^{*}={\displaystyle \frac{k_i}{m{r}_i}}{o}^{*}\left(r_i+{\gamma }_i(m-o^{*})\right)fori=1,2,3.$

Proof. 

To find the equilibrium points, we equal to zero the right-hand sides of system (1) and, thus, obtain the following homogeneous system of nonlinear algebraic equations

(3)$\begin{array}{ccc}\hfill 0& =& {\displaystyle r\left(1-{\displaystyle \frac{o^2}{m^2}}\right)-o\sum _{i=1}^3{\beta }_i{x}_i,}\hfill \end{array}$

(4)$\begin{array}{ccc}\hfill 0& =& r_i{x}_i\left(1-{\displaystyle \frac{m{x}_i}{k_{i}o}}\right)+{\gamma }_i(m-o)x_i,\mathrm{for}i=1,2,3.\hfill \end{array}$

Equation (4) can be written as follows

$0=x_i\left[r_i\left(1-{\displaystyle \frac{m{x}_i}{k_{i}o}}\right)+{\gamma }_i(m-o)\right],\mathrm{for}i=1,2,3.$

Which are satisfied if $x_i=0$ or

(5)$\begin{array}{ccc}\hfill x_i& =& {\displaystyle \frac{k_i}{m{r}_i}}o\left(r_i+{\gamma }_i(m-o)\right),\mathrm{for}i=1,2,3.\hfill \end{array}$

When substituting $x_i=0$, for all $i=1,2,3$, in Equation (3), and solve it; it is obtained that $o=m$. Thus, an equilibrium point of the system is $E_0=(m,0,0,0)$, and as $m>0$, then $E_0\in \Omega$, so the proposition has been proven for (i). On the other hand, if $x_i\ne 0$, substituting expression (5) in Equation (3) we have

(6)$\begin{array}{ccc}\hfill 0& =& {\displaystyle \frac{1}{m}}\sum _{i=1}^3{\displaystyle \frac{{\beta }_i{\gamma }_i{k}_i}{r_i}}{o}^3-\left({\displaystyle \frac{r}{m^2}}+{\displaystyle \frac{1}{m}}\sum _{i=1}^3{\beta }_i{k}_i+\sum _{i=1}^3{\displaystyle \frac{{\beta }_i{\gamma }_i{k}_i}{r_i}}\right)o^2+r.\hfill \end{array}$

Note that Equation (6) is a polynomial of degree three in the variable o and we intend to calculate its roots. Then, we define the function

(7)$\begin{array}{ccc}\hfill f\left(o\right)& =& {\displaystyle \frac{1}{m}}\sum _{i=1}^3{\displaystyle \frac{{\beta }_i{\gamma }_i{k}_i}{r_i}}{o}^3-\left({\displaystyle \frac{r}{m^2}}+{\displaystyle \frac{1}{m}}\sum _{i=1}^3{\beta }_i{k}_i+\sum _{i=1}^3{\displaystyle \frac{{\beta }_i{\gamma }_i{k}_i}{r_i}}\right)o^2+r.\hfill \end{array}$

Since $f\left(o\right)$ is a polynomial, then it is a continuous function, particularly in $\left[0,m\right]$; furthermore, $f\left(0\right)=r>0$ and

$f\left(m\right)=-m\sum _{i=1}^3{\beta }_i{k}_i<0,\mathrm{so} \mathrm{that}f\left(0\right)·f\left(m\right)<0,$

consequently, by Bolzano’s Theorem, there exists $o^{*}\in (0,m)$, such that $f\left(o^{*}\right)=0$; thus, we guarantee the existence of a positive real solution in $(0,m)$ for $f\left(o\right)$. In addition, when $o\to +\infty$, it follows that $f\left(o\right)\to +\infty$; on the other hand, if $o\to -\infty$, then $f\left(o\right)\to -\infty$; therefore, the function $f\left(o\right)$ has two positive real solutions, namely $o^{*}\in (0,m)$ and $o^{**}>m$; and one negative real solution, $o^{***}<0$. Note that for roots $o^{**}>m$ and $o^{***}<0$, the equilibrium points with first coordinate $o^{**}>m$ or $o^{***}<0$, are not in $\Omega$.

Finally, for the root $o^{*}\in (0,m)$, substituting it in the expression (5), we obtain

(8)$x_i={\displaystyle \frac{k_i}{m{r}_i}}{o}^{*}\left(r_i+{\gamma }_i(m-o^{*})\right)>0,\mathrm{for}i=1,2,3.$

Thus, any equilibrium point having $o^{*}$ as the first coordinate, and at least one of the coordinates $x_i$, as in (5), will be an equilibrium point of system (1), since it satisfies (3) and (4), and will also be in $\Omega$. If we denote as $x_i^{*}$ the equilibrium values of the expression (8), the equilibrium points $E_1,E_2,\dots ,E_7\in \Omega$, and, thus, ($ii$) and ($iii$) are proven. □

The solution $E_0$ and absence of all classes of AMs represent oxygen depletion or supersaturation in aquatic ecosystems, when the DO tends to zero or when it tends to m, respectively. From expression (3) we can see that, in the absence of AMs, since $o=m$ is an equilibrium value of this scenario, the DO tends asymptotically to the value m, and in the presence of at least one AM family, the level of DO will be less than m.

According to expression (5), we have that

(9)$\underset{o\to m}{lim}{x}_i=\underset{o\to m}{lim}{\displaystyle \frac{k_i}{m{r}_i}}o\left(r_i+{\gamma }_i(m-o)\right)=k_i\mathrm{para}i=1,2,3,$

that is, if o tends to m, the equilibrium values $x_i$ will tend to $k_i$. Furthermore, if $o\in (0,m)$ and $x_i\ne 0$ for some $i=1,2,3$, then

(10)${\displaystyle \frac{k_i}{m{r}_i}}o\left(r_i+{\gamma }_i(m-o)\right)<k_i.$

Expression (9) indicates that the equilibrium value of each $x_i$ tends to the value $k_i$ for $i=1,2,3$, while relation (10) suggests that such equilibrium values will not be able to reach the value of $k_i$; summarizing these two results, we can say that if o tends to the value m, then the equilibrium values $x_i$ tend asymptotically to the values $k_i$ below, for each $i=1,2,3$.

The equilibrium point $E_7=(o^{*},x_1,x_2,x_3)$, having all AM classes is associated with water bodies with significant richness and abundance, because of this; water quality, as a function of $E_7$, accepts the following classifications shown in Table 2.

The strict inequalities, described in Table 2, between the different classes of organisms implies that there is a dominant class; therefore, the one that determines the quality of the water. On the other hand, if the abundances are approximately equal ($x_1\approx x_2\approx x_3$), then the water body is diverse and uniform, which implies that it tends to recover much faster than one that is not. On the other hand, the system requires the presence or absence of different classes of AMs and their relationship with the DO, i.e., the loss of diversity and abundance. Note that each equilibrium solution, explicit in the presence equilibria of the previous proposition, represents the loss of diversity of AM families, in the function of the DO concentration, which allows establishing, in a more specific way, the stress level of a freshwater ecosystem. For example, the absence of Classes II and III families in the equilibrium solution $E_1$ implies an ecosystem with sufficient DO to sustain the Class I organism communities, the absence of pollution, and excellent water quality. In $E_2$, the absence of Class III organisms implies two cases: (a) if $x_1<x_2$, then there is regular pollution with a tendency to decrease and the water quality is good and (b) if $x_2<x_1$, then there is low pollution with a tendency to increase and the water quality is good. Note that the absence of Classes I and II AMs in the equilibrium solution $E_6$ allows inferring that there is high pollution, and the water quality is very poor. The degree of contamination and water quality for the equilibrium solutions $E_3$, $E_4$, and $E_5$, can be inferred analogously.

The importance of analyzing the stability of equilibrium solutions lies in knowing under what conditions an ecosystem is altered, it returns to its original state. In addition, it is possible to infer which parameters have the greatest influence on the dynamics of the system. This analysis is developed in the following section.

3.2. Stability Analysis

Holling defines ecological resilience as the capacity of a system to absorb disturbances and reorganize while undergoing changes in order to essentially preserve its function, structure, identity, and feedback. In the literature, the main mathematical definitions of resilience are based on the dynamical systems theory, and specifically on attractors and basins of attraction [52]. In this study, the necessary conditions for system (1), after some perturbations, to return to its original state, are established in the following result:

Proposition 2.

Given the system of nonlinear differential Equation (1), the AM-free equilibrium point, $E_0=(m,0,0,0)$, is unstable. In addition, the coexistence equilibrium point, $E_7=(o^{*},x_1^{*},x_2^{*},x_3^{*})$, will be locally stable if $m\le \frac{r_i}{{\gamma }_i}$, for $i=1,2,3$.

Proof. 

To analyze the stability of the equilibrium points, we will use the indirect Lyapunov method. For this purpose, we return to the autonomous system of differential equations given by (1):

(11)$\begin{array}{ccc}\hfill o^{\prime }& =& f_0(o,x_1,x_2,x_3)=r\left(1-{\displaystyle \frac{o^2}{m^2}}\right)-o\sum _{i=1}^3{\beta }_i{x}_i,\hfill \end{array}$

(12)$\begin{array}{ccc}\hfill x_i^{\prime }& =& f_i(o,x_1,x_2,x_3)=r_i{x}_i\left(1-{\displaystyle \frac{m{x}_i}{k_{i}o}}\right)+{\gamma }_i(m-o)x_i,\mathrm{for}i=1,2,3.\hfill \end{array}$

The partial derivatives of (11), with respect to o and $x_i$, are given by:

$\frac{\partial f_0}{\partial o}=-\frac{2r}{m^2}o-\sum _{i=1}^3{\beta }_i{x}_i$

and

$\frac{\partial f_0}{\partial x_i}=-{\beta }_{i}o\mathrm{for}i=1,2,3.$

On the other hand, the partial derivatives of (12), with respect to o and $x_i$, are given by:

$\frac{\partial f_i}{\partial o}={\displaystyle \frac{m{r}_i{x}_i^2}{k_i{o}^2}}-{\gamma }_i{x}_i$

and

$\frac{\partial f_i}{\partial x_j}=\left\{\begin{array}{cc}{r}_i\left(1-\frac{2m}{k_{i}o}{x}_i\right)+{\gamma }_i(m-o)& \mathrm{if}i=j\\ 0& \mathrm{if}i\ne j\end{array}\right.$

for $i=1,2,3$.

Thus, the Jacobian matrix of system (1) is shown in Appendix A.1.

Evaluating the Jacobian matrix at the equilibrium point $E_0$ we obtain

(13)$J{\mid }_{E_0}=\left(\begin{array}{cccc}-\frac{2r}{m}& -{\beta }_{1}m& -{\beta }_{2}m& -{\beta }_{3}m\\ 0& r_1& 0& 0\\ 0& 0& r_2& 0\\ 0& 0& 0& r_3\end{array}\right).$

Note that the matrix in (13) is upper triangular, so its eigenvalues are the elements of the diagonal. Since this matrix has eigenvalues ${\lambda }_1=-\frac{2r}{m},{\lambda }_2=r_1,{\lambda }_3=r_2$ and ${\lambda }_4=r_3$, three of which are positive, then, the equilibrium point $E_0$ is unstable and this proves the first statement.

The Jacobian matrix evaluated at the equilibrium point $E_7$ is given by

(14)$J{\mid }_{E_7}=\left(\begin{array}{cccc}-S& -{\beta }_1{o}^{*}& -{\beta }_2{o}^{*}n& -{\beta }_3{o}^{*}\\ S_1& -T_1& 0& 0\\ S_2& 0& -T_2& 0\\ S_3& 0& 0& -T_3\end{array}\right),$

where

(15)$\begin{array}{ccc}\hfill S& =& \frac{2r}{m^2}{o}^{*}+\sum _{i=1}^3{\beta }_i{\displaystyle \frac{k_i{o}^{*}}{m{r}_i}}{T}_i,\hfill \end{array}$

(16)$\begin{array}{ccc}\hfill T_i& =& r_i+{\gamma }_i(m-o^{*})\mathrm{for}i=1,2,3,\hfill \end{array}$

(17)$\begin{array}{ccc}\hfill S_i& =& {\displaystyle \frac{k_i{T}_i}{m{r}_i}}\left(T_i-o^{*}{\gamma }_i\right)\mathrm{for}i=1,2,3.\hfill \end{array}$

Computing the determinant $\mid J-\lambda I\mid$ for the matrix (14), we have

(18)$D=\left|\begin{array}{cccc}-S-\lambda & -{\beta }_1{o}^{*}& -{\beta }_2{o}^{*}& -{\beta }_3{o}^{*}\\ S_1& -T_1-\lambda & 0& 0\\ S_2& 0& -T_2-\lambda & 0\\ S_3& 0& 0& -T_3-\lambda \end{array}\right|$

Consequently, the characteristic equation of the matrix (18) is given by

(19)$\begin{array}{c}\hfill {\lambda }^4+a{\lambda }^3+b{\lambda }^2+c\lambda +d=0,\end{array}$

and their coefficients are denoted in Appendix A.2.

The first observations are $o^{*}<m$ for being the root of Equation (6) with that property, from this we have that $T_i>0$ for each $i=1,2,3$, and $S>0$ for how it is defined. Now, it will be proven that, if for each $i=1,2,3$, $m\le \frac{r_i}{{\gamma }_i}$, then $S_i>0$. Suppose that for each $i=1,2,3$ it is satisfied that $m\le \frac{r_i}{{\gamma }_i}$, then as $o^{*}<m$, you have to $o^{*}<\frac{r_i}{{\gamma }_i}$, but this implies that ${\gamma }_i{o}^{*}<r_i$, i.e., $2(r_i-o^{*}{\gamma }_i)>0$, so that $r_i+{\gamma }_i(\frac{r_i}{{\gamma }_i}-o^{*})-o^{*}{\gamma }_i>0$, i.e., $T_i-o^{*}{\gamma }_i>0$, which implies that $S_i>0$ for $i=1,2,3$.

With this, we have proved that the coefficients of Equation (19) satisfy $a,b,c,d>0$. Moreover, analogous to the previous statement we can easily prove that $abc-c^2-a^{2}d>0$; therefore, Expression (19) is a Hurwitz polynomial. Now, we will apply the Routh–Hurwitz criterion to the characteristic Equation (19) to demonstrate the stability of the equilibrium point $E_7$. The Hurwitz matrices with the coefficients of Equation (19) are determined by:

$H_1=\left(a\right),H_2=\left(\begin{array}{cc}a& 1\\ c& b\end{array}\right),H_3=\left(\begin{array}{ccc}a& 1& 0\\ c& b& a\\ 0& d& c\end{array}\right)\mathrm{and}{H}_4=\left(\begin{array}{cccc}a& 1& 0& 0\\ c& b& a& 1\\ 0& d& c& b\\ 0& 0& 0& d\end{array}\right).$

Hence, with $a,b,c$, and d given in the form of Appendix A.2, the determinants of the matrices: $H_1$, $H_2$, $H_3$, and $H_4$ are greater than zero. Therefore, if $m\le \frac{r_i}{{\gamma }_i}$, for $i=1,2,3$, all eigenvalues of the characteristic Equation (19) have the negative real part and, consequently, the equilibrium point $E_7$ is locally stable. □

Proposition 2 establishes that the stability of system (1) can be studied through the equilibrium solution $E_7$. This solution allows inferring that the average utilization rate (${\gamma }_i$), oxygen saturation (m), and average reproduction rate ($r_i$) are the most important parameters in the variation of the system dynamics. From an ecological point of view, if $m\le \frac{r_i}{{\gamma }_i}$, for $i=1,2,3$, then the coexistence equilibrium is stable, implying that regardless of the initial conditions (abundances) in which the ecosystem is found, it will be stable, as long as this relationship is satisfied, in addition to being able to obtain information on long-term water quality through AM population fluctuations.

4. Results

In this section, the influence of oxygenation/aeration rate (r) on dissolved oxygen concentration (o) and different population densities of aquatic macroinvertebrates (state variables $x_1,x_2$, and $x_3$) are studied, as well as the effects of the coexistence equilibrium point, $E_7$, on water quality. Figure 2 shows the behaviors of the state variables $o,x_1,x_2$, and $x_3$ for different oxygenation/aeration rates. In general, when r is less than $7\%$, the AM abundances decrease (red band in b–d of Figure 2). Moreover, it takes more time for these abundances to reach their maximum and be constant. If r is between 20 and 80%, the maximum is reached in less time and takes less time to be constant. It should be noted that when the oxygenation rate is greater than 3 mgL⁻¹, there is no significant difference in the behavior of the population densities (blue band in b–d of Figure 2).

Figure 3 shows the population fluctuations of AMs under the equilibrium conditions established in Proposition 2. Moreover, Figure 3a indicates that, for all time t, the most pollution-tolerant AM families (Class III) are the most abundant, so it can be inferred that the water quality is poor. Through the behavior of the AM populations shown in Figure 3b, it can be inferred that for $0<t<30$, the water quality is bad, but for $t>30$, the water quality is regular because the most abundant AM families are Class II, which means that there is a transition in two different time intervals, at the same time that under those initial conditions, the ecosystem changes when the oxygenation/aeration rate increases. Likewise, Figure 3c shows that for $t<19$ the water quality is poor; however, for $9<t<68$ the water quality is fair, and for $t>68$, the water quality is good, this is inferred because the most abundant AM families are Class I (non-tolerant to pollution).

Validation

Rosas-Acevedo et al. [53] evaluated the water quality of the Coyuca Lagoon, Mexico, through environmental stress indicator insects. They made monthly visits to 10 sampling stations to collect insects between April 2011 and April 2012, collecting 2907 organisms, of which, 1343 were AM insects that were indicators of environmental stress. AMs are distributed as follows: 82 Class I organisms (Ephemeroptera, Trichoptera, Plecoptera, and Megaloptera), 1091 Class II organisms (Notonectidae, Diptera, Odonata, Hemiptera, Coleoptera, Chironomida), and 188 Class III organisms (Diptera), in other words, $x_1<x_3<x_2$. In addition, they reported that the minimum and maximum values observed for DO were 1.5 and 4 mg/L, respectively. From the above, they concluded that the dominant organisms were Class II, followed by Class III. While for the observed DO values, Class I organisms tended to disappear. They concluded that there was moderate contamination with a tendency to increase.

The logistic population growth model was fitted to the observed data using the least squares method, assuming an asymptotic trend response to the carrying capacity of the indicator organisms. The dashed curves, shown in Figure 4, are the logistic growth curves fitted to the observed data, using Equation (20) with ${\widehat{k}}_1\approx 0$, ${\widehat{r}}_1=-0.1758$, ${\widehat{k}}_2\approx 98$, ${\widehat{r}}_2=0.3387$, ${\widehat{k}}_3\approx 17$, and ${\widehat{r}}_3=0.1585$ for Classes I, II, and III, respectively.

(20)$y\left(t\right)=\frac{\widehat{k}*P_0}{P_0+(\widehat{k}-P_0)*{exp}^{-\widehat{r}t}}$

Note that $x_1<x_3<x_2$ allows inferring that there is moderate contamination, with a tendency to increase. The logistic growth curves do not consider the influence of DO. Likewise, the simulated population growth curves that consider the influence of the DO (proposed in this work) are shown in Figure 4. Observe that the population fluctuations in the observed vs. simulated data have similar behaviors but at different time intervals. For example, in the observed data, the abundance of Class I AMs is higher than that of Class III for $t<5$. Whereas for $t>5$, this relationship is reversed. On the other hand, in the simulated data, these same relationships hold for $t<15$ and $t\ge 15$, respectively (see Figure 4). It should be noted that the influence of oxygenation/aeration prevents the loss of diversity and abundance in short intervals of time.

5. Discussion

Oxygenation and aeration of a water body have proven to be useful for ecosystem restoration [54], for example, in Durowskie Lake, water oxygenation by wind aerators contributed to the increase in oxygen concentration, richness, and abundance of macroinvertebrates [55], which coincides with that reported by Kail et al. [56]. From Figure 2, the changes over time of the different state variables were simulated and it is evident that, the increase or decrease, the oxygenation/aeration rate, r, significantly affects the population dynamics of the organisms and dissolved oxygen concentration.

The qualitative analysis of the model determined analytically its equilibrium solutions. Proposition 1 established a set of equilibrium points belonging to the positive invariant region $\Omega$. The existence of the coexistence equilibrium solution, $E_7$, stands out for its importance in the evaluation of water quality. As a consequence and as a function of $E_7$, it was possible to classify the water quality through the dominant AM families. In the same direction, the equilibrium points described in part two of the same proposition classify the stress level of an ecosystem (water quality) as more strictly associated with the loss of diversity.

The stability of the equilibrium points demonstrated in Proposition 2 made it possible to identify those parameters that bring the system to stable or unstable behaviors, since it is essential to arrive at conclusions of value for researchers or managers since stability represents the ability to return to the state of equilibrium after temporary disturbances, which are often unpredictable. From Figure 3 it is observed that, regardless of the initial conditions of the state variables, the populations of organisms tend to stabilize when the conditions necessary for the equilibrium point, $E_7$, to be stable are satisfied. Moreover, this positively or negatively affects water quality when there are dominant classes of AMs.

Model validation, shown in Figure 4, evidenced the behavior of AM populations, with and without the influence of dissolved oxygen, through the rate of oxygenation/aeration. Highlighting the benefits of long-term oxygenation/aeration for the abundance and richness of bioindicator organisms. On the other hand, to strengthen the study of water quality, in subsequent works, taxa could be considered individually and take advantage of their functional responses or qualities of AMs as bioindicators, although this implies a greater number of equations and more complex analyses. Furthermore, since the conditions established here are general, they could be extrapolated to the study of different ecosystems with bioindicators with similar dynamics, adjusting the parameters to particular contexts.

6. Conclusions

Aquatic ecosystems involve many challenges for their restoration, understanding as systems, and analyses after being disturbed. The conceptualization of the model proposed here was formulated with the natural dynamics between two important environmental variables for the evaluation of water quality—bioindicator organisms and dissolved oxygen.

In this study, the changes over time of the different state variables were simulated and it is evident that the oxygenation/aeration rate, r, significantly affected the population dynamics of the organisms. In addition, knowing the current conditions of an ecosystem makes it possible to understand its current quality and, thus, establish a one-to-one correspondence between water quality and equilibrium points. On the other hand, the simulations showed that under equilibrium conditions (and regardless of initial conditions), the carrying capacity and oxygenation/aeration rate have the greatest influence over time in determining water quality.

Therefore, the study of water quality through aquatic macroinvertebrates and their close relationship with dissolved oxygen using the predictive model proposed here made it possible to show future scenarios that help guarantee water availability and its sustainable management.

Author Contributions

The authors contributed equally to this work. J.J.P.-P., J.C.H.-G., Y.E.M.-G., J.M.-R., and J.M.S. conceived, designed, and performed the numerical experiments; J.C.H.-G., J.M.-R., J.J.P.-P., Y.E.M.-G., and J.M.S. analyzed the data and wrote the paper. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

J.J.P.-P. thanks the Internationalization of Research BUAP, for the support received for the development of this work. In addition, to the Vicerrectoría de Investigación y Estudios de Posgrado (VIEP) BUAP and CONACYT for the complementary support through the Sistema Nacional de Investigadores (SNI).

Conflicts of Interest

The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

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

Figures and Tables

Figure 1. Conceptual diagram of the mathematical model. The state variables o, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] represent the DO concentrations in the water and the AM populations—intolerant, regularly tolerant, and tolerant to organic contamination, respectively. The processes are: (a) oxygenation/aeration, (b) DO consumption, (c) AM natality, (d) DO utilization, and (e) interspecific competition.

Figure 2. Influence of oxygenation/aeration rate r on AM population densities. (a) Variation of DO through time for different values of r. The population fluctuations of the different classes of AMs are shown in (b), (c), and (d) for [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], respectively. Note that, in all cases, if r is smaller, then AM reproduction is smaller and slower than when r is larger.

Figure 3. Behavior of AM population densities. In all cases; [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and the initial conditions are the same, except (a) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (poor water quality), (b) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (regular water quality), and (c) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (good water quality).

Figure 4. Population growth curves of observed vs. simulated data. The dotted lines do not consider the influence of the DO, while the continuous line curves do. Data from [53].

Table 1

Description and parameter values of system (1).

Parameter	Description	Value	Reference
r	oxygenation/aeration rate	0–14.5 mgL−1	Galic et al. [23]
m	oxygen saturation constant a	0 < m ≤ 14.5 mgL−1
βi	average respiration rate of AMs a	0–0.918	Galic et al. [23]
ri	average reproduction rate of AMs a	0–1.02	Galic et al. [23]
ki	average carrying capacity of AMs	varies
γi	average DO utilization rate	varies

^a—Assuming an oxygen saturation response for each of the rates (see Figure A1 in Appendix B).

Table 2

Relationships between the cardinalities of different AM classes, meaning, and classification of water quality.

Relation	Meaning	Water Quality
x3 < x2 < x1	Low pollution with a tendency to increase	good
x3 < x1 < x2	Regular pollution with a tendency to decrease	moderate
x2 < x3 < x1	Low pollution with a tendency to increase	good
x2 < x1 < x3	High pollution with a tendency to decrease	poor
x1 < x3 < x2	Regular pollution with a tendency to increase	moderate
x1 < x2 < x3	High pollution with a tendency to decrease	poor

Appendix A. Mathematical Calculation

Appendix A.1. Jacobian Matrix

The Jacobian of system (1) is calculated as follows:$$J=\left(\begin{array}{cccc}-{\displaystyle \frac{2r}{m^2}}o-\sum _{i=1}^3{\beta }_i{x}_i& -{\beta }_{1}o& -{\beta }_{2}o& -{\beta }_{3}o\\ {\displaystyle \frac{m{r}_1{x}_1^2}{k_1{o}^2}}-{\gamma }_1{x}_1& r_1\left(1-{\displaystyle \frac{2m}{k_{1}o}}{x}_1\right)+{\gamma }_1(m-o)& 0& 0\\ {\displaystyle \frac{m{r}_2{x}_2^2}{k_1{o}^2}}-{\gamma }_2{x}_2& 0& r_2\left(1-{\displaystyle \frac{2m}{k_{2}o}}{x}_2\right)+{\gamma }_2(m-o)& 0\\ {\displaystyle \frac{m{r}_1{x}_1^2}{k_1{o}^2}}-{\gamma }_1{x}_1& 0& 0& r_3\left(1-{\displaystyle \frac{2m}{k_{3}o}}{x}_3\right)+{\gamma }_3(m-o)\end{array}\right).$$

Appendix A.2. Coefficients

The coefficients of polynomial (19) are determined by:$$\begin{array}{ccc}\hfill a& =& S+T_1+T_2+T_3,\hfill \\ \hfill b& =& T_1{T}_2+T_3{T}_2+T_1{T}_3+S\left(T_1+T_2+T_3\right)+o^{*}\left({\beta }_1{S}_1+{\beta }_2{S}_2+{\beta }_3{S}_3\right),\hfill \\ \hfill c& =& T_1{T}_2{T}_3+S\left(T_1{T}_2+T_3{T}_2+T_1{T}_3\right)+o^{*}\left({\beta }_1{S}_1{T}_2+{\beta }_1{S}_1{T}_3+{\beta }_2{S}_2{T}_1+{\beta }_2{S}_2{T}_3+{\beta }_3{S}_3{T}_1+{\beta }_3{S}_3{T}_2\right),\hfill \\ \hfill d& =& S{T}_1{T}_2{T}_3+o^{*}\left({\beta }_1{S}_1{T}_2{T}_3+{\beta }_2{S}_2{T}_1{T}_3+{\beta }_3{S}_3{T}_1{T}_2\right).\hfill \end{array}$$

Appendix B. Parameters

Figure A1 shows AM respiration and reproduction rates when DO availability increases or decreases.

Figure A1. Impacts of hypoxia on reproduction and respiration of aquatic invertebrates. Data digitized at different oxygen concentrations, scaled to control conditions (normoxia) and the Michaelis–Menten model fit (±95% CI). Adapted from Galic et al. [23].

References

1. United Nations. Transforming our world: The 2030 agenda for sustainable development. General Assembley 70 Session; United Nations: New York, NY, USA, 2015.

2. Uddin, M.G.; Nash, S.; Olbert, A.I. A review of water quality index models and their use for assessing surface water quality. Ecol. Indic.; 2021; 122, pp. 3400-3411.

3. Almeida, C.; González, S.O.; Mallea, M.; González, P. A recreational water quality index using chemical, physical and microbiological parameters. Environ. Sci. Pollut. Res.; 2012; 19, pp. 3400-3411. [DOI: https://dx.doi.org/10.1007/s11356-012-0865-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22528988]

4. Farhat, N.; Kim, L.H.; Vrouwenvelder, J.S. Online characterization of bacterial processes in drinking water systems. Npj Clean Water; 2020; 3, pp. 1-7. [DOI: https://dx.doi.org/10.1038/s41545-020-0065-7]

5. Jørgensen, S.E.; Xu, F.L.; Salas, F.; Marques, J.C. Handbook of Ecological Indicators for Assessment of Ecosystem Health; CRC Press: Boca Raton, FL, USA, 2016; pp. 9-64.

6. Hernández-Mira, F.A.; Rosas-Acevedo, J.L.; Reyes-Umaña, M.; Violante-González, J.; Sigarreta-Almira, J.M.; Vakhania, N. Multimetric Index to Evaluate Water Quality in Lagoons: A Biological and Geomorphological Approach. Sustainability; 2021; 13, 4631. [DOI: https://dx.doi.org/10.3390/su13094631]

7. Bain, M.B.; Stevenson, N.J. Aquatic Habitat Assessment: Commom Methods; Asian Fisheries Society: Bethesda, MD, USA, 1999; 186. ISBN 1-888569-18-2

8. Ejigu, M.T. Overview of water quality modeling. Cogent Eng.; 2021; 8, 1891711. [DOI: https://dx.doi.org/10.1080/23311916.2021.1891711]

9. Zhen-Gang, J. Hydrodynamics and Water Quality: Modeling Rivers, Lakes, and Estuaries; John Wiley & Sons: Hoboken, NJ, USA, 2017.

10. Loucks, D.P.; Van Beek, E. Water Resource Systems Planning and Management: An Introduction to Methods, Models, and Applications; Springer: Cham, Switzerland, 2017.

11. Gonenc, I.E.; Wolflin, J.P. Coastal Lagoons: Ecosystem Processes and Modeling for Sustainable Use and Development; CRC Press: Boca Raton, FL, USA, 2004.

12. Jacoby, J.; Welch, E. Pollutant Effects in Freshwater: Applied Limnology; 3rd ed. CRC Press: Boca Raton, FL, USA, 2004.

13. Sánchez, E.; Colmenarejo, M.F.; Vicente, J.; Rubio, A.; García, M.G.; Travieso, L.; Borja, R. Use of the water quality index and dissolved oxygen deficit as simple indicators of watersheds pollution. Ecol. Indic.; 2007; 7, pp. 315-328. [DOI: https://dx.doi.org/10.1016/j.ecolind.2006.02.005]

14. Kannel, P.R.; Lee, S.; Lee, Y.S.; Kanel, S.R.; Khan, S.P. Application of water quality indices and dissolved oxygen as indicators for river water classification and urban impact assessment. Environ. Monit. Assess.; 2007; 132, pp. 93-110. [DOI: https://dx.doi.org/10.1007/s10661-006-9505-1]

15. Best, M.A.; Wither, A.W.; Coates, S. Dissolved oxygen as a physico-chemical supporting element in the Water Framework Directive. Mar. Pollut. Bull.; 2007; 55, pp. 53-64. [DOI: https://dx.doi.org/10.1016/j.marpolbul.2006.08.037]

16. Kannel, P.R.; Kanel, S.R.; Lee, S.; Lee, Y.S.; Gan, T.Y. A review of public domain water quality models for simulating dissolved oxygen in rivers and streams. Environ. Model. Assess.; 2011; 16, pp. 183-204. [DOI: https://dx.doi.org/10.1007/s10666-010-9235-1]

17. Kaller, M.D.; Kelso, W.E. Association of macroinvertebrate assemblages with dissolved oxygen concentration and wood surface area in selected subtropical streams of the southeastern USA. Aquat. Ecol.; 2007; 41, pp. 95-110. [DOI: https://dx.doi.org/10.1007/s10452-006-9046-2]

18. Wilhm, J.; McClintock, N. Dissolved oxygen concentration and diversity of benthic macroinvertebrates in an artificially destratified lake. Hydrobiologia; 1978; 57, pp. 163-166. [DOI: https://dx.doi.org/10.1007/BF00016460]

19. Wang, J.; Fu, Z.; Qiao, H.; Liu, F. Assessment of eutrophication and water quality in the estuarine area of Lake Wuli, Lake Taihu, China. Sci. Total Environ.; 2019; 650, pp. 1392-1402. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.09.137] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30308826]

20. Meier, H.E.M.; Eilola, K.; Almroth-Rosell, E.; Schimanke, S.; Kniebusch, M.; Höglund, A.; Pemberton, P.; Liu, Y.; Väli, G.; Saraiva, S. Disentangling the impact of nutrient load and climate changes on Baltic Sea hypoxia and eutrophication since 1850. Clim. Dyn.; 2019; 53, pp. 1145-1166. [DOI: https://dx.doi.org/10.1007/s00382-018-4296-y]

21. Zhang, J.; Wang, C.; Jiang, X.; Song, Z.; Xie, Z. Effects of human-induced eutrophication on macroinvertebrate spatiotemporal dynamics in Lake Dianchi, a large shallow plateau lake in China. Environ. Sci. Pollut. Res.; 2020; 27, pp. 1-15. [DOI: https://dx.doi.org/10.1007/s11356-020-07773-w]

22. Bazzanti, M.; Mastrantuono, L.; Solimini, A.G. Selecting macroinvertebrate taxa and metrics to assess eutrophication in different depth zones of Mediterranean lakes. Fundam. Appl. Limnol. Hydrobiol.; 2012; 180, pp. 133-143. [DOI: https://dx.doi.org/10.1127/1863-9135/2012/0200]

23. Galic, N.; Hawkins, T.; Forbes, V.E. Adverse impacts of hypoxia on aquatic invertebrates: A meta-analysis. Sci. Total Environ.; 2019; 652, pp. 736-743. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.10.225]

24. Etemi, F.Z.; Bytyçi, P.; Ismaili, M.; Fetoshi, O.; Ymeri, P.; Shala–Abazi, A.; Muja-Bajraktari, N.; Czikkely, M. The use of macroinvertebrate based biotic indices and diversity indices to evaluate the water quality of Lepenci river basin in Kosovo. J. Environ. Sci. Heal. Part A; 2020; 55, pp. 748-758. [DOI: https://dx.doi.org/10.1080/10934529.2020.1738172]

25. Slimani, N.; Sánchez-Fernández, D.; Guilbert, E.; Boumaiza, M.; Guareschi, S.; Thioulouse, J. Assessing potential surrogates of macroinvertebrate diversity in North-African Mediterranean aquatic ecosystems. Ecol. Indic.; 2019; 101, pp. 324-329. [DOI: https://dx.doi.org/10.1016/j.ecolind.2019.01.017]

26. Silva, D.R.O.; Herlihy, A.T.; Hughes, R.M.; Macedo, D.R.; Callisto, M. Assessing the extent and relative risk of aquatic stressors on stream macroinvertebrate assemblages in the neotropical savanna. Sci. Total Environ.; 2018; 633, pp. 179-188. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.03.127]

27. Croijmans, L.; de Jong, J.F.; Prins, H.H. Oxygen is a better predictor of macroinvertebrate richness than temperature—A systematic review. Environ. Res. Lett.; 2020; 38, pp. 1820-1832. [DOI: https://dx.doi.org/10.1088/1748-9326/ab9b42]

28. Su, P.; Wang, X.; Lin, Q.; Peng, J.; Song, J.; Fu, J.; Wang, S.; Cheng, D.; Bai, H.; Li, Q. Variability in macroinvertebrate community structure and its response to ecological factors of the Weihe River Basin, China. Ecol. Eng.; 2019; 140, 105595. [DOI: https://dx.doi.org/10.1016/j.ecoleng.2019.105595]

29. Mezgebu, A.; Lakew, A.; Lemma, B. Water quality assessment using benthic macroinvertebrates as bioindicators in streams and rivers around Sebeta, Ethiopia. Afr. J. Aquat. Sci.; 2019; 44, pp. 361-367. [DOI: https://dx.doi.org/10.2989/16085914.2019.1685450]

30. Liu, Z.; Fan, B.; Huang, Y.; Yu, P.; Li, Y.; Chen, M.; Cai, M.; Lv, W.; Jiang, Q.; Zhao, Y. Assessing the ecological health of the Chongming Dongtan Nature Reserve, China, using different benthic biotic indices. Mar. Pollut. Bull.; 2018; 146, pp. 76-84.

31. Liu, Z.; Chen, M.; Li, Y.; Huang, Y.; Fan, B.; Lv, W.; Yu, P.; Wu, D.; Zhao, Y. Different effects of reclamation methods on macrobenthos community structure in the Yangtze Estuary, China. Mar. Pollut. Bull.; 2018; 127, pp. 429-436.

32. Pineda-Pineda, J.J.; Rosas-Acevedo, J.L.; Sigarreta, J.M.; Hernández-Gómez, J.C.; Reyes-Umaña, M. Biotic Indices to Evaluate Water Quality: BMWP. Int. J. Environ. Ecol. Fam. Urban Stud. IJEEFUS; 2018; 8, pp. 23-36.

33. Zhou, X.D.; Xu, M.Z.; Lei, F.K.; Zhang, J.H.; Wang, Z.Y.; Luo, Y.Y. Responses of Macroinvertebrate Assemblages to Flow in the Qinghai-Tibet Plateau: Establishment and Application of a Multi-metric Habitat Suitability Model. Water Resour. Res.; 2022; 58, e2021WR030909. [DOI: https://dx.doi.org/10.1029/2021WR030909]

34. Pineda-Pineda, J.J.; Martínez-Martínez, C.T.; Méndez-Bermúdez, J.A.; Muñoz-Rojas, J.; Sigarreta, J.M. Application of Bipartite Networks to the Study of Water Quality. Sustainability; 2020; 12, 5143. [DOI: https://dx.doi.org/10.3390/su12125143]

35. Schleiter, I.M.; Borchardt, D.; Wagner, R.; Dapper, T.; Schmidt, K.D.; Schmidt, H.H.; Werner, H. Modelling water quality, bioindication and population dynamics in lotic ecosystems using neural networks. Ecol. Model.; 1999; 120, pp. 271-286. [DOI: https://dx.doi.org/10.1016/S0304-3800(99)00108-8]

36. Maier, H.R.; Dandy, G.C. The use of artificial neural networks for the prediction of water quality parameters. Water Resour. Res.; 1996; 32, pp. 1013-1022. [DOI: https://dx.doi.org/10.1029/96WR03529]

37. Villamarín, C.; Rieradevall, M.; Paul, M.J.; Barbour, M.T.; Prat, N. A tool to assess the ecological condition of tropical high Andean streams in Ecuador and Peru: The IMEERA index. Ecol. Indic.; 2013; 29, pp. 79-92. [DOI: https://dx.doi.org/10.1016/j.ecolind.2012.12.006]

38. Karaouzas, I.; Smeti, E.; Kalogianni, E.; Skoulikidis, N.T. Ecological status monitoring and assessment in Greek rivers: Do macroinvertebrate and diatom indices indicate same responses to anthropogenic pressures?. Ecol. Indic.; 2019; 101, pp. 126-132. [DOI: https://dx.doi.org/10.1016/j.ecolind.2019.01.011]

39. El Sayed, S.M.; Hegab, M.H.; Mola, H.R.; Ahmed, N.M.; Goher, M.E. An integrated water quality assessment of Damietta and Rosetta branches (Nile River, Egypt) using chemical and biological indices. Environ. Monit. Assess.; 2020; 192, pp. 1-16. [DOI: https://dx.doi.org/10.1007/s10661-020-8195-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32162005]

40. Streeter, H.W.; Phelps, E.B. A Study of the Pollution and Natural Purification of the Ohio River; Public Health Bulletin No 146 Public Health Service: Washington, DC, USA, 1925.

41. Wang, Q.; Li, S.; Jia, P.; Qi, C.; Ding, F. A review of surface water quality models. Sci. World J.; 2013; [DOI: https://dx.doi.org/10.1155/2013/231768]

42. da Silva Burigato Costa, C.M.; Leite, I.R.; Almeida, A.K.; de Almeida, I.K. Choosing an appropriate water quality model—A review. Environ. Monit. Assess.; 2021; 193, pp. 1-15. [DOI: https://dx.doi.org/10.1007/s10661-020-08786-1]

43. Park, R.A.; Clough, J.S.; Wellman, M.C. AQUATOX: Modeling environmental fate and ecological effects in aquatic ecosystems. Ecol. Model.; 2008; 213, pp. 1-15. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2008.01.015]

44. West Virginia Department of Environmental Protection. Available online: https://dep.wv.gov/WWE/getinvolved/sos/Documents/Benthic/VisualMacroGuide.pdf (accessed on 1 September 2022).

45. Researchgate. Available online: https://www.researchgate.net/profile/Pablo-Gutierrez-Fonseca/publication/295854904_Guia_fotografica_de_familias_de_macroinvertebrados_acuaticos_de_Puerto_Rico/links/56ce23e508aeb52500c36b4f/Guia-fotografica-de-familias-de-macroinvertebrados-acuaticos-de-Puerto-Rico.pdf (accessed on 1 September 2022).

46. Everaert, G.; De Neve, J.; Boets, P.; Dominguez-Granda, L.; Mereta, S.T.; Ambelu, A.; Hoang, T.H.; Goethals, P.L.M.; Thas, O. Comparison of the abiotic preferences of macroinvertebrates in tropical river basins. PLoS ONE; 2014; 9, e108898. [DOI: https://dx.doi.org/10.1371/journal.pone.0108898] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25279673]

47. Weiwei, L.; Youhui, H.; Zhiquan, L.; Yang, Y.; Bin, F.; Yunlong, Z. Application of macrobenthic diversity to estimate ecological health of artificial oyster reef in Yangtze Estuary, China. Mar. Pollut. Bull.; 2016; 103, pp. 137-143.

48. Weiwei, L.; Zhiquan, L.; Yang, Y.; Youhui, H.; Bin, F.; Qichen, J.; Yunlong, Z. Loss and self-restoration of macrobenthic diversity in reclamation habitats of estuarine islands in Yangtze Estuary, China. Mar. Pollut. Bull.; 2016; 103, pp. 128-136.

49. Pearl, R.; Reed, L.J. On the rate of growth of the population of the United States since 1790 and its mathematical representation. Proc. Natl. Acad. Sci. USA; 1920; 6, 275. [DOI: https://dx.doi.org/10.1073/pnas.6.6.275]

50. Hutchinson, G.E. An Introduction to Population Ecology; Yale University Press: New Haven, CT, USA, 1978.

51. Chapman, E.J.; Byron, C.J. The flexible application of carrying capacity in ecology. Glob. Ecol. Conserv.; 2018; 13, e00365.

52. Gunderson, L.H. Ecological resilience-in theory and application. Annu. Rev. Ecol. Syst.; 2000; 31, pp. 425-439. [DOI: https://dx.doi.org/10.1146/annurev.ecolsys.31.1.425]

53. Rosas-Acevedo, J.L.; Ávila-Pérez, H.; Sánchez-Infante, A.; Rosas-Acevedo, Y.; García-Ibañez, S.; Sampedro-Rosas, L.; Granados-Ramírez, J.G.; Juárez-López, A.L. índice BMWP, FBI y EPT para determinar la calidad del agua en la laguna de Coyuca de Benítez, Guerrero, México. Rev. Iberoam. Cienc.; 2014; 1, pp. 81-88.

54. Rosińska, J.; Kozak, A.; Dondajewska, R.; Kowalczewska-Madura, K.; Gołdyn, R. Water quality response to sustainable restoration measures–Case study of urban Swarzędzkie Lake. Ecol. Indic.; 2018; 84, pp. 437-449. [DOI: https://dx.doi.org/10.1016/j.ecolind.2017.09.009]

55. Gołdyn, R.; Podsiadłowski, S.; Dondajewska, R.; Kozak, A. The sustainable restoration of lakes—Towards the challenges of the Water Framework Directive. Ecohydrol. Hydrobiol.; 2014; 14, pp. 68-74. [DOI: https://dx.doi.org/10.1016/j.ecohyd.2013.12.001]

56. Kail, J.; Brabec, K.; Poppe, M.; Januschke, K. The effect of river restoration on fish, macroinvertebrates and aquatic macrophytes: A meta-analysis. Ecol. Indic.; 2015; 58, pp. 311-321. [DOI: https://dx.doi.org/10.1016/j.ecolind.2015.06.011]