Next Article in Journal
Quantum Integral Inequalities in the Setting of Majorization Theory and Applications
Previous Article in Journal
Automatic Parking Path Optimization Based on Immune Moth Flame Algorithm for Intelligent Vehicles
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Continuous and Discrete Dynamical Models of Total Nitrogen Transformation in a Constructed Wetland: Sensitivity and Bifurcation Analysis

1
Department of Mathematics, Faculty of Science and Mathematics, Universitas Diponegoro, Jl. Prof. Jacub Rais, Semarang 50275, Indonesia
2
Department of Physics, Faculty of Science and Mathematics, Universitas Diponegoro, Jl. Prof. Jacub Rais, Semarang 50275, Indonesia
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(9), 1924; https://doi.org/10.3390/sym14091924
Submission received: 10 August 2022 / Revised: 8 September 2022 / Accepted: 9 September 2022 / Published: 14 September 2022

Abstract

:
In this research, we study a dynamical system of total nitrogen transformation in a mangrove-filled constructed wetland. The system’s variables are the mangrove biomass concentration and total nitrogen concentration in wastewater and in soil solution. We investigate the system’s dynamics by examining the local stability of the equilibriums, simulating the phase portrait and solutions and providing time-dependent parameter sensitivity analyses. The analysis shows that the level of garbage acts as the parameter for when mangrove biomass will disappear. Both the graphs of the system solutions and the sensitivity function in the case of biomass concentration and total nitrogen concentration in soil solution versus time show symmetrical features at specific time intervals. According to the sensitivity index when reaching equilibrium, the level of garbage is the most sensitive parameter to the system. In addition, we explore the model’s discrete form by investigating the conditions for the equilibrium’s local stability and presenting bifurcation diagrams for each parameter. The symmetrical aspects are visible in the visualization of the bifurcation diagram and the solutions’ chaotic behavior.

1. Introduction

Wastewater as an environmental problem has been on the rise as a result of the global effects of urbanization, climate change, and population increase [1]. Most wastewater contains nitrogen compounds, and if the wastewater is overly polluted with nitrogen, the water will smell and be unfit for human consumption [2]. Traditional biological wastewater treatment methods still have issues that prevent them from fully meeting the new environmental standards, particularly in terms of removing nitrogen compounds, which cause eutrophication in the water and seriously jeopardizes the stability of aquatic ecosystems [3]. Meanwhile, wastewater treatment via constructed wetlands as sustainable engineered systems to reduce water pollution in many different forms is spreading throughout the world during recent years [4,5]. The constructed wetlands are viewed as a substitute for traditional treatment due to its lower energy needs, ease of use and maintenance, and superior effectiveness by having microbes consume pollution produced by industrial and agricultural activity [6,7].
Constructed wetlands are a type of nature-based solution that may effectively and methodically imitate nature and its inherent traits and functions, and they can be used for a variety of land uses, including urban, agricultural, animal, and forest land uses [1]. Constructed wetlands incorporate a variety of wetland ecosystem functions via physical and ecological features, including biomass generation, habitat provision, cooling, and water purification [8,9]. The process of purifying water involves eliminating nutrients such as nitrogen, which is mostly accomplished by bacterial activity and plant uptake [1].
One of the primary nitrogen removal mechanisms in constructed wetlands is thought to involve the processes of denitrification [10]. The addition of carbon sources such as agricultural biomass materials can increase the removal of total nitrogen and significantly enhance the quality of constructed wetlands’ effluent in order to speed up the denitrification process [11]. Agricultural biomass materials are good optional carbon media because of their superior ability to supply carbon and low cost [12]. Mangroves are one of the plants that are frequently employed in constructed wetlands to generate high-quality agricultural biomass. They are self-sustaining and have unique responses to a variety of environmental stresses [13,14]. Thus, planting mangroves is one of the natural methods that may be used in constructed wetlands to improve nitrogen removal from wastewater.
One way to observe the dynamics of nitrogen removal in a constructed wetland is by using mathematical modeling such as a dynamical system. Bunwong et al. [2] proposed a mathematical model that describes the transformation of total nitrogen in a mangrove-filled constructed wetland. The model is in the form of a system of differential equations with three compartments: total nitrogen concentration in wastewater, total nitrogen concentration in soil solution, and mangrove biomass concentration. The model assumed that the growth rate of mangrove biomass production will follow the Michaelis–Menten kinetics, which illustrate the idea of limiting nutrients, i.e., there is a causal relationship between its exhaustion and the termination of growth [15], as evidenced by some empirical studies [16,17,18]. The model also assumed that the conversion of biomass is represented by a linear yield function following the suggestion in [19,20]. The analysis of the model was performed by scaling, that is all variables were scaled by the constant input rate of total nitrogen concentration. The authors also examined the existence of Hopf bifurcation for the yield function’s parameter.
The model of Bunwong et al. [2] was used and developed in the following studies. In [21], a nonlinear mathematical model of the rhizosphere microbial degradation in wetland plants was built based on impulsive state feedback control to examine the dynamics of pollutant concentration removal in constructed wetlands. The authors in [22] proposed a mathematical model of rhizosphere microbial degradation with impulsive diffusion on nutrients, which involves organic concentrations inside and outside the rhizosphere. In [23], the growth of wetland plants was examined in relation to rhizosphere dispersal and impulse input by employing the nutrient–wetland plant interaction in reducing the pollutant concentration. In order to explain the complex degradation and movement caused by introducing helpful microorganisms into the plant rhizosphere, the authors in [24] develop a three-dimensional diffusion model of the competitive rhizosphere microbes with impulsive feedback control.
In this research, we use continuous and discrete models to investigate the dynamics of total nitrogen transformation in a mangrove-filled constructed wetland. The continuous model employs the model of Bunwong et al. [2], but with a constant yield. We performed sensitivity analysis on the model’s parameters to determine how they affect the variables. By examining the sensitivity index after reaching steady-state circumstances, we may determine which parameter is the most sensitive to the system. The discrete model is derived by discretizing the continuous model. The bifurcation diagram of each parameter is presented, and the chaotic behavior is observed.

2. Dynamics of Total Nitrogen Transformation

We followed the model developed by [2] to describe the dynamics of total nitrogen transformation in a constructed wetland. Suppose T ( t ) is the concentration of mangrove biomass, W ( t ) is the concentration of total nitrogen in the wastewater, and  S ( t ) is the concentration of total nitrogen in the soil solution. The dynamical model of total nitrogen transformation in the wetland is given below:
d T ( t ) d t = β S ( t ) T ( t ) k + S ( t ) σ T ( t ) d W ( t ) d t = Q ( γ + α ) W ( t ) d S ( t ) d t = γ W ( t ) 1 a β S ( t ) T ( t ) k + S ( t ) ϕ S ( t )
where all the initial values are positive: T ( 0 ) > 0 , W ( 0 ) > 0 , S ( 0 ) > 0 . The variables are all non-negative, and the parameters are all positive. The model’s diagram is given in Figure 1. The arrows represent the movement of matter through the system, and each arrow also displays the parameterizations of these flows’ rates.
The explanation of Figure 1 is as follows. The mangroves take the total nitrogen from the source and convert it into biomass. It is supposed that a filter prevents debris from falling into the pond. The concentration of the total nitrogen input rate is a constant, denoted as Q. The loss of total nitrogen in the wastewater due to runoff or evaporation is α W . The total nitrogen exchange rate between the wastewater and soil solution is γ W . The growth rate of the mangroves follows the Michaelis–Menten kinetics β S k + S . The ratio of nutrients consumed to the biomass generated by the mangroves is represented by the yield constant a. In [2], the yield coefficient was assumed to be a linear function of S. However, in this paper, we studied the model with a constant yield. Through leaching or denitrification, the total nitrogen is lost in soil solutions at the rate ϕ S . The litter fall rate is σ T .

2.1. Equilibrium Points and Local Stability Analysis

We denote the equilibrium points of the system (1) by E * = ( T * , W * , S * ) . From  d W d t = 0 , we obtain the equilibrium W * = Q γ + α . Meanwhile, d T d t = 0 produces the equilibriums S * = σ k β σ (with β > σ ) or T * = 0 . From this case, we have three cases.
If d T d t = 0 happens when T * = 0 and S * σ k β σ , then
0 = d S d t = γ W ϕ S S * = γ ϕ W * = γ Q ϕ ( γ + α ) .
Thus, we have
E 1 * = 0 , Q γ + α , γ Q ϕ ( γ + α ) .
If d T d t = 0 happens when T * = 0 and S * = σ k β σ simultaneously, then
0 = d S d t = γ W ϕ S γ ϕ W * = S * or γ Q ϕ ( γ + α ) = σ k β σ .
In this case, we obtain
E 2 * = 0 , Q γ + α , σ k β σ ,
with γ Q σ ( γ + α ) = ϕ k β σ . This means that one of the parameters depends on the others.
If d T d t = 0 happens when T * 0 and S * = σ k β σ , then
0 = d S d t = γ W 1 a β S T k + S ϕ S T * = a γ Q σ ( γ + α ) ϕ k β σ .
Thus, we obtain
E 3 * = a γ Q σ ( γ + α ) ϕ k β σ , Q γ + α , σ k β σ ,
with the positivity conditions for the equilibrium: γ Q σ ( γ + α ) > ϕ k β σ and β > σ .
To analyze the local stability of each equilibrium, first, we calculate the Jacobian matrix of the system (1) at any point E = ( T , W , S ) , which is given below:
J ( E ) = β S k + S σ 0 β k T ( k + S ) 2 0 ( γ + α ) 0 β S a ( k + S ) γ β k T a ( k + S ) 2 + ϕ
The stability of each equilibrium is given by the following theorems.
Theorem 1.
E 1 * is stable if σ > σ ¯ , where σ ¯ = β γ Q k ϕ ( γ + α ) + γ Q .
Proof. 
Let σ ¯ = β γ Q k ϕ ( γ + α ) + γ Q . The Jacobian matrix with respect to E 1 * is
J ( E 1 * ) = σ ¯ σ 0 0 0 ( γ + α ) 0 σ ¯ a γ ϕ ,
with the characteristic equation
σ ¯ σ λ ( γ + α ) λ ϕ λ = 0 .
Thus, we have λ 1 = σ ¯ σ , λ 2 = ( γ + α ) , and  λ 3 = ϕ . It is clear that λ 2 and λ 3 are less than zero. Meanwhile, λ 1 will be negative if σ > σ ¯ . This proves the theorem.    □
In other words, Theorem 1 says that the system can be stable with a zero concentration of mangrove biomass ( T * = 0 ) when the level of garbage is too high ( σ > σ ¯ ).
Theorem 2.
E 2 * is not stable.
Proof. 
The Jacobian matrix with respect to E 2 * is
J ( E 2 * ) = 0 0 0 0 ( γ + α ) 0 σ a γ ϕ ,
with the characteristic equation
λ ( ( γ + α ) λ ) ( ϕ λ ) = 0 .
This means there is a zero eigenvalue, which implies that the system is not stable.    □
Theorem 3.
E 3 * is stable.
Proof. 
The Jacobian matrix with respect to E 3 * is
J ( E 3 * ) = 0 0 a A 0 ( γ + α ) 0 σ a γ ( A + ϕ ) ,
where A = ( β σ ) 2 k β γ Q σ ( γ + α ) ϕ k β σ . The characteristic equation of J ( E 3 * ) is
λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0 ,
where
a 1 = γ + α + ϕ + A ,                        
a 2 = ( γ + α ) ( A + ϕ ) + σ A ,
a 3 = σ A ( γ + α ) .                                    
By the Routh–Hurwitz criterion, the equilibrium E 3 * is stable if a 0 , a 1 , a 2 > 0 and a 1 a 2 > a 0 . Note that all coefficients of the characteristic equation are positive. It can be seen that a 1 a 2 = ( γ + α ) ( A + ϕ ) + σ A ( γ + α + ϕ + A ) = ( γ + α ) ( A + ϕ ) ( γ + α + ϕ + A ) + σ A ( ϕ + A ) + σ A ( γ + α ) > σ A ( γ + α ) = a 0 . Therefore, by the Routh–Hurwitz criterion, E 3 * is stable.    □

2.2. Numerical Solution

Numerical simulations of the solution of the system (1) were performed using the parameter values as presented in Table 1. The phase portraits of the system (1) with numerous different initial values are presented in Figure 2. The solution-pairs ( T ( t ) , W ( t ) , S ( t ) ) behave in two ways: converging to the equilibriums E 1 * and E 3 * . The graphical solutions that converge to E 1 * and E 3 * are given in Figure 3. The solutions will converge to E 1 * (the concentration of mangrove biomass is zero) if the level of garbage is too high ( σ > σ ¯ ); otherwise, the solutions will converge to E 3 * (positive equilibrium). Based on Table 1, then we have σ ¯ = 0.2143 . In the simulations of Figure 3a,b, we used σ = 0.25 . In Figure 3, we can also conclude that, if the exchange rate of total nitrogen between wastewater and soil solution is bigger than the loss of total nitrogen in the soil solution ( γ > ϕ ), then the concentration of total nitrogen in the wastewater will converge to a value that is lower than the concentration of the total nitrogen in the soil solution. On the other hand, if  γ < ϕ (here, we used γ = 0.2 and ϕ = 0.3 ), then the concentration of total nitrogen in the wastewater will converge to a value that is higher than the concentration of the total nitrogen in the soil solution. Another interesting result is the appearance of symmetrical characteristics between the mangrove biomass concentration and total nitrogen concentration in the soil solution at specific intervals where the two complement each other.

2.3. Sensitivity Analysis

A sensitivity analysis of a dynamical system can provide rich qualitative behavior by observing the impacts of the parameters on the variables over time, and we can determine which parameter is the most sensitive to the system. This technique is widely used in the literature [25,26,27,28,29].
Let the vector of variables be X = ( T , W , S ) and the vector of parameters be P = ( β , k , σ , Q , γ , α , a , ϕ ) , and the right side in the model (1) is F = ( F , G , H ) , where F = β S T k + S σ T , G = Q ( γ + α ) W , and  H = γ W 1 a β S T k + S ϕ S . Define a sensitivity function V = X / P . Next, we have
d V d t = d d t X P = P d X d t = F X X P + F P = J ( X ) V + F P .
Note that J ( X ) is a 3 × 3 Jacobian matrix given in (5), V is a 3 × 8 matrix, and  F / P is a 3 × 8 matrix.
We denote T β = v β T , , T ϕ = v ϕ T , W β = v β W , , W ϕ = v ϕ W , and  S β = v β S , , S ϕ = v ϕ S . Thus, we have
V = v β T v k T v σ T v Q T v γ T v α T v a T v ϕ T v β W v k W v σ W v Q W v γ W v α W v a W v ϕ W v β S v k S v σ S v Q S v γ S v α S v a S v ϕ S .
For the case of F P , we have
F P = S T k + S β S T ( k + S ) 2 T 0 0 0 0 0 0 0 0 1 W W 0 0 1 a S T k + S 1 a β S T ( k + S ) 2 0 0 W 0 1 a 2 β S T k + S S .
Each pair of respective entries in the matrix Equation (9) forms a differential equation. Therefore, from (9), we have a system of twenty four differential equations (the details of the equations are provided in Appendix A). We simulated the sensitivity by solving (9) numerically using ode45 in Matlab with initial values V 0 = 0 3 × 8 ( 3 × 8 zeroes matrix). The simulation was performed with the parameter values in Table 1. The results are presented in Figure 4.
In the model (1), the differential equation of total nitrogen concentration in wastewater (W) is a standalone or independent equation. Therefore, the only parameters that affect W are the parameters that appear in the equation, that is Q, γ , and α . In other words, W is only sensitive with the change of Q, γ , and α , not for the other parameters. This fact was confirmed by the simulation of the parameters’ sensitivity in Figure 4. In Figure 4d, the parameter Q affects W positively, while in Figure 4e,f, both parameters γ and α affect W negatively. This is obvious because of the differential equation form of W.
For the case of the parameters of the Michaelis–Menten function, in Figure 4a, we can see that the parameter β affects T positively and S negatively, and on the other hand, the parameter k affects T negatively and S positively, as shown in Figure 4b. For the case of the other parameters, there are interesting findings regarding the effects of the changes of the parameters on the variables T and S. In Figure 4c, the parameter of the garbage level σ affects T negatively, which is obvious, and S positively. In Figure 4d to Figure 4h, all the parameters Q, γ , α , a, and ϕ affect S positively and then negatively or otherwise from the starting time until some time point, but then, the sensitivity converges to zero. In other words, in the long term, those parameters will not affect the total nitrogen concentration in the soil solution at all. For the rest of the cases, the parameters Q, γ , and a affect T positively, as shown in Figure 4d,e,g, while the parameters α and ϕ both affect T negatively, as shown in Figure 4f,h.
The majority of the graphical sensitivity functions display the symmetrical properties of the parameter sensitivity in the case of the mangrove biomass concentration and total nitrogen concentration in the soil solution. Another interesting result from the sensitivity analysis is that we can observe which parameter has the greatest effect on the system. In other words, we can know which parameter is most sensitive. In Figure 5, we present the sensitivity index of the parameters after arriving at the steady-state conditions. From the figure, we can see that the level of garbage parameter σ is the most sensitive parameter.

3. Discrete Form Model

Discrete dynamical systems are commonly utilized to simulate the dynamics of numerous natural phenomena. Occasionally, the discrete model is generated by discretizing the continuous model, because discrete models are easier to manage. One advantage of the discrete model is that we can easily display a bifurcation diagram and observe chaos, although the analysis may not be easy. Several studies have been conducted on discrete dynamical systems and the bifurcation of the system’s parameters [30,31,32,33,34].
The discrete form of the model (1) is given by
T t + 1 = T t + β S t T t k + S t σ T t W t + 1 = W t + Q ( γ + α ) W t S t + 1 = S t + γ W t 1 a β S t T t k + S t ϕ S t ,
with discrete time t = 0 , 1 , 2 , .
The equilibriums of the discrete system (10) are the same as the equilibriums of the system (1), that is E 1 * as in (2), E 2 * as in (3), and E 3 * as in (4). Meanwhile, the Jacobian matrix of the system (10) at any point E = ( T , W , S ) is given by
J ( E ) = 1 + β S k + S σ 0 β k T [ k + S ] 2 0 1 ( γ + α ) 0 β S a [ k + S ] γ 1 β k T a [ k + S ] 2 ϕ
The local stability of equilibriums E 1 * and E 2 * is given in the following theorems.
Theorem 4.
E 1 * is stable if σ ¯ < σ < 2 + σ ¯ , γ + α < 2 , and ϕ < 2 , where σ ¯ = β γ Q k ϕ ( γ + α ) + γ Q .
Proof. 
The Jacobian matrix with respect to E 1 * is
J ( E 1 * ) = 1 + σ ¯ σ 0 0 0 1 ( γ + α ) 0 σ ¯ a γ 1 ϕ ,
with the characteristic equation
( 1 σ ¯ σ λ ) ( 1 ( γ + α ) λ ) ( 1 ϕ λ ) = 0 .
The eigenvalues are λ 1 = 1 + σ ¯ σ , λ 2 = 1 ( γ + α ) , and λ 3 = 1 ϕ . The equilibrium E 1 * is stable if all the eigenvalues satisfy | λ 1 , 2 , 3 | < 1 . We have
λ 1 < 1 if σ > σ ¯ , λ 2 < 1 , λ 3 < 1 ,
and
λ 1 > 1 if σ < 2 + σ ¯ , λ 2 > 1 if γ + α < 2 , λ 3 > 1 if ϕ < 2 .
Thus, we deduce that E 1 * is stable if σ ¯ < σ < 2 + σ ¯ , γ + α < 2 , and ϕ < 2 . □
Note that, since σ ¯ = β γ Q k ϕ ( γ + α ) + γ Q , then the condition σ ¯ < σ < 2 + σ ¯ can be changed into the condition for the other parameter. For example, in terms of parameter β , we have the condition
( σ 2 ) k ϕ ( γ + α ) + γ Q γ Q < β < σ k ϕ ( γ + α ) + γ Q γ Q .
Since β > 0 , then the left-hand side of (12) can be replaced by zero if the value of σ is less than 2. The case of the other parameters can be performed in a similar way.
Theorem 5.
E 2 * is not stable.
Proof. 
The Jacobian matrix with respect to E 2 * is
J ( E 2 * ) = 1 0 0 0 1 ( γ + α ) 0 σ a γ 1 ϕ ,
with the characteristic equation:
( 1 λ ) ( 1 ( γ + α ) λ ) ( 1 ϕ λ ) = 0 .
Thus, there is an eigenvalue that satisfies | λ | = 1 . Therefore, E 2 * is not stable. □

3.1. The Local Stability of Non-Zero Equilibrium

Now, we analyze the local stability of equilibrium E 3 * . The Jacobian matrix with respect to E 3 * is
J ( E 3 * ) = 1 0 a A 0 1 ( γ + α ) 0 σ a γ 1 ( A + ϕ ) ,
where A = ( β σ ) 2 k β γ Q σ ( γ + α ) ϕ k β σ .
The characteristic equation of the Jacobian matrix (the details of the calculation are given in Appendix B) is
λ 3 + b 1 λ 2 + b 2 λ + b 3 = 0 ,
where b 1 = a 1 3 , b 2 = a 2 2 a 1 + 3 , b 3 = a 3 + a 1 a 2 1 , where a 1 , a 2 , a 3 are defined in (6)–(8).
According to the Jury criterion [35], E 3 * is stable if the coefficients of the characteristic polynomial (13) satisfy
1 + b 1 + b 2 + b 3 > 0 , 1 b 1 + b 2 b 3 > 0 , 1 b 2 + b 1 b 3 b 3 2 > 0 , b 2 > 3 .
In the previous section, it was shown that all a 1 , a 2 , and a 3 are positive. The left-hand sides of the Jury criterion (14) are evaluated in terms of a 1 , a 2 , and a 3 :
1 + b 1 + b 2 + b 3 = a 3 ,
1 b 1 + b 2 b 3 = 8 4 a 1 + 2 a 2 a 3 ,
1 b 2 + b 1 b 3 b 3 2 = 3 2 a 3 a 2 2 a 3 2 + a 1 a 2 + 2 a 2 a 3 a 1 a 3 ,
b 2 = a 2 2 a 1 + 3 .
It is clear that Equation (15) is positive. Equation (16) will be positive if 4 a 1 2 a 2 + a 3 < 8 . Equation (17) will be positive if 2 a 3 + a 2 2 + a 3 2 a 1 a 2 2 a 2 a 3 + a 1 a 3 < 3 . Lastly, Equation (18) will be greater than 3 if 2 a 1 < a 2 .

3.2. Numerical Simulations

We performed some numerical simulations to depict the behavior of the equilibrium of the system (10) when the parameters vary. This behavior can be represented by a bifurcation diagram. All the simulations were performed using the parameter values in Table 1. In the following figures, we present the bifurcation diagrams of all parameters of the model (10): β , k, σ , Q, γ , α , a, and ϕ , respectively.
First, we depict the bifurcation diagram of parameter β in Figure 6a. We can see that the equilibrium L * has zero value for the value of β in the range 0 until some number that is bigger than 0.2, but less than 0.3. This value is none other than σ k ϕ ( γ + α ) + γ Q γ Q as in (12). Thus, we deduce that the system converges to E 1 * when 0 < β < σ k ϕ ( γ + α ) + γ Q γ Q . When β > σ k ϕ ( γ + α ) + γ Q γ Q , the system starts to converge to equilibrium E 3 * . However, as we can see in the figure, the system becomes unstable when the value of β is around 0.5. Perhaps it is possible to obtain this value analytically from the stability condition in Section 3.1, but it may be difficult to deal with. With the same reasoning, in the rest of this paper, we only present the numerical simulations. Back to Figure 6a, after the equilibrium E 3 * is not stable, the system produces periodic orbits, and then, at some value of β , the system becomes chaotic. The behavior of the solution of the system (10) when it becomes chaotic is depicted in Figure 6b. This simulation used β = 0.8 . From this result, we can conclude that the system may be unstable when the parameter β is quite big.
Next, we present the bifurcation diagram of parameter k in Figure 7a. The figure shows that the system becomes unstable and leads to chaos when the value of k is too small (near zero). The simulation of the chaotic behavior of the system when k = 0.15 is given in Figure 7b.
The bifurcation of parameter σ is presented in Figure 8a. We can see that the system converges to equilibrium E 1 * when σ is bigger than some value between 0.2 and 0.25. This value is none other than σ ¯ . When σ < σ ¯ , the system converges to equilibrium E 3 * . However, when σ is quite small, we can observe that the system becomes unstable and leads to chaos. The chaotic behavior of the solution when σ = 0.03 is shown in Figure 8b. By this result, we can deduce that too small a value of σ may make the system become unstable.
From the equilibrium of the total nitrogen concentration in the wastewater W * = Q γ + α , it is clear that a higher Q makes a higher W * . This is evidenced by Figure 9a. From the figure, we can also observe that the system converges to equilibrium E 1 * when Q is too low, and a higher Q makes the system converge to equilibrium E 3 * . However, at some point, too high a Q will cause the system to become unstable and even lead to chaos. The behavior of the system when it becomes chaotic is shown in Figure 9b. The simulation used Q = 5.3 .
For parameters γ and α , it is obvious that the increase in either of these parameters’ values will make the equilibrium W * decrease. However, when these parameter values are too high, as we can observe in Figure 10a and Figure 11a, they cause not only W * to possibly be unstable, but also the system to possibly become unstable. The unstable behavior of this system is simulated in Figure 10a and Figure 11a for γ = 1.83 and α = 1.72 , respectively. We can see that the system swings, with bigger fluctuations as time passes.
For the case of parameter a, as shown in Figure 12a,b, we do not see any interesting result when the bifurcation diagram is plotted. This happens because parameter a does not appear in the Jacobian matrices; thus, it does not contribute to the stability condition. The only clear result is that a higher a results in a higher equilibrium T * .
The bifurcation diagram of parameter ϕ is presented in Figure 13a. Most of the figure does not show any interesting dynamics. However, when we depict the behavior of the solution with ϕ = 1.96 , as given in Figure 13b, we obtain a quite interesting dynamics of the solutions T ( t ) and S ( t ) . They oscillate first before reaching the equilibrium.

4. Discussion and Conclusions

Even with an abundance of nutrients, plants in nature are limited in their ability to develop. The use of a Michaelis–Menten function in the model of Bunwong et al. [2] for the growth of mangrove biomass is capable of representing the limitation. According to the function, as the total nitrogen concentration in the soil solution increases, the mangrove biomass growth will converge to the upper bound. However, we demonstrated that the growth of the mangrove biomass, as well as the total nitrogen concentration in the soil solution become unstable and may cause chaos in the discrete model with a specific high level of β or low level of k.
The level of garbage in the constructed wetland makes a big contribution to determining the stability of the total nitrogen transformation system. If the level of garbage becomes too high, the mangrove biomass will disappear. On the other hand, when the level of garbage is too low, from the discrete model, we showed that the system may be unstable and can lead to chaos. The level of garbage affects the mangrove biomass and the total nitrogen concentration in the soil solution over time, although it does not affect the total nitrogen concentration in the wastewater. This parameter also is the most sensitive parameter of the model, as shown by the sensitivity index of all the parameters after reaching the steady-state conditions. Based on these results, the level of garbage in a constructed wetland is the most essential aspect that must be carefully maintained so that the nitrogen removal in the wastewater remains stable.
We can see symmetrical features in both the continuous and discrete models for the mangrove biomass concentration and total nitrogen concentration in the soil solution. In the case of the continuous model, the symmetry emerges at some interval times in the graph of the model’s solution and in the graph of the parameter sensitivity function. Meanwhile, in the case of the discrete model, the symmetry appears in the bifurcation diagrams of some parameters and in the chaotic behavior of the solution.

Author Contributions

Conceptualization, S. and M.F.A.; methodology, M.F.A.; formal analysis, M.F.A.; investigation, S. and M.F.A.; resources, S.; writing—original draft preparation, M.F.A.; writing—review and editing, S., M.F.A., S.K. and D.P.S.; visualization, M.F.A.; supervision, S. and D.P.S.; project administration, S.K.; funding acquisition, S. All authors have read and agreed to the published version of the manuscript.

Funding

The authors are grateful to DRPM KEMENDIKBUD-RISTEK RI for funding the research, through the PDUPT Grant 2022, with Contract No. 187-21/UN7.6.1/PP/2022.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The details of the system of differential equations (9) are provided below.
d v β T d t = β S k + S σ v β T + β k T ( k + S ) 2 v β S + S T k + S , v β T ( 0 ) = 0 ;    
d v k T d t = β S k + S σ v k T + β k T ( k + S ) 2 v k S β S T ( k + S ) 2 , v k T ( 0 ) = 0 ;    
d v σ T d t = β S k + S σ v σ T + β k T ( k + S ) 2 v σ S T , v σ T ( 0 ) = 0 ;  
d v Q T d t = β S k + S σ v Q T + β k T ( k + S ) 2 v Q S , v Q T ( 0 ) = 0 ;
d v γ T d t = β S k + S σ v γ T + β k T ( k + S ) 2 v γ S , v γ T ( 0 ) = 0 ;  
d v α T d t = β S k + S σ v α T + β k T ( k + S ) 2 v α S , v α T ( 0 ) = 0 ;  
d v a T d t = β S k + S σ v a T + β k T ( k + S ) 2 v a S , v a T ( 0 ) = 0 ;  
d v ϕ T d t = β S k + S σ v ϕ T + β k T ( k + S ) 2 v ϕ S , v ϕ T ( 0 ) = 0 ;  
d v β W d t = ( γ + α ) v β W , v β W ( 0 ) = 0 ;        
d v k W d t = ( γ + α ) v k W , v k W ( 0 ) = 0 ;        
d v σ W d t = ( γ + α ) v σ W , v σ W ( 0 ) = 0 ;        
d v Q W d t = ( γ + α ) v Q W + 1 , v Q W ( 0 ) = 0 ;    
d v γ W d t = ( γ + α ) v γ W W , v γ W ( 0 ) = 0 ;  
d v α W d t = ( γ + α ) v α W W , v α W ( 0 ) = 0 ;  
d v a W d t = ( γ + α ) v a W , v a W ( 0 ) = 0 ;  
d v ϕ W d t = ( γ + α ) v ϕ W , v ϕ W ( 0 ) = 0 ;  
d v β S d t = β S a ( k + S ) v β T + γ v β W β k T a ( k + S ) 2 + ϕ v β S 1 a S T k + S , v β S ( 0 ) = 0 ;
d v k S d t = β S a ( k + S ) v k T + γ v k W β k T a ( k + S ) 2 + ϕ v k S + 1 a β S T ( k + S ) 2 , v k S ( 0 ) = 0 ;
d v σ S d t = β S a ( k + S ) v σ T + γ v σ W β k T a ( k + S ) 2 + ϕ v σ S , v σ S ( 0 ) = 0 ;    
d v Q S d t = β S a ( k + S ) v Q T + γ v Q W β k T a ( k + S ) 2 + ϕ v Q S , v Q S ( 0 ) = 0 ;    
d v γ S d t = β S a ( k + S ) v γ T + γ v γ W β k T a ( k + S ) 2 + ϕ v γ S + W , v γ S ( 0 ) = 0 ;    
d v α S d t = β S a ( k + S ) v α T + γ v α W β k T a ( k + S ) 2 + ϕ v α S , v α S ( 0 ) = 0 ;
d v a S d t = β S a ( k + S ) v a T + γ v a W β k T a ( k + S ) 2 + ϕ v a S + 1 a 2 β S T k + S , v a S ( 0 ) = 0 ;        
d v ϕ S d t = β S a ( k + S ) v ϕ T + γ v ϕ W β k T a ( k + S ) 2 + ϕ v ϕ S S , v ϕ S ( 0 ) = 0 .

Appendix B

The detailed calculation to obtain Equation (13) is given below. Let A = ( β σ ) 2 k β γ Q σ ( γ + α ) ϕ k β σ .
0 = d e t 1 0 a A 0 1 ( γ + α ) 0 σ a γ 1 ( A + ϕ ) λ 0 0 0 λ 0 0 0 λ = 1 ( γ + α ) λ ( 1 λ ) ( 1 ( A + ϕ ) λ ) + σ A = 1 ( γ + α ) λ λ 2 + ( A + ϕ 2 ) λ + 1 + σ A ( A + ϕ ) .
Multiplying both sides with 1 gives
0 = λ + ( γ + α ) 1 λ 2 + ( A + ϕ 2 ) λ + 1 + σ A ( A + ϕ ) = λ 3 + γ + α + ϕ + A 3 λ 2 + ( γ + α ) ( A + ϕ ) + σ A 2 ( γ + α + ϕ + A ) + 3 λ + σ A ( γ + α ) + γ + α + ϕ + A ( γ + α ) ( A + ϕ ) + σ A 1 .
In (6)–(8), we denote a 1 = γ + α + ϕ + A , a 2 = ( γ + α ) ( A + ϕ ) + σ A , and a 3 = σ A ( γ + α ) . Thus, (A25) can be simplified as
λ 3 + ( a 1 3 ) λ 2 + ( a 2 2 a 1 + 3 ) λ + a 3 + a 1 a 2 1 = 0 .
By denoting b 1 = a 1 3 , b 2 = a 2 2 a 1 + 3 , b 3 = a 3 + a 1 a 2 1 , Equation (A26) becomes
λ 3 + b 1 λ 2 + b 2 λ + b 3 = 0 .

References

  1. Choi, H.; Geronimo, F.K.; Jeon, M.; Kim, L.H. Evaluation of bacterial community in constructed wetlands treating different sources of wastewater. Ecol. Eng. 2022, 182, 106703. [Google Scholar] [CrossRef]
  2. Bunwong, K.; Sae-jie, W.; Lenbury, Y. Modelling nitrogen dynamics of a constructed wetland: Nutrient removal process with variable yield. Nonlinear Anal. Theory Methods Appl. 2009, 71, e1538–e1546. [Google Scholar] [CrossRef]
  3. Wu, H.; Zhang, J.; Ngo, H.H.; Guo, W.; Hu, Z.; Liang, S.; Fan, J.; Liu, H. A review on the sustainability of constructed wetlands for wastewater treatment: Design and operation. Bioresour. Technol. 2015, 175, 594–601. [Google Scholar] [CrossRef] [PubMed]
  4. Senzia, M.A.; Mashauri, D.A.; Mayo, A.W. Suitability of constructed wetlands and waste stabilisation ponds in wastewater treatment: Nitrogen transformation and removal. Phys. Chem. Earth Parts A/B/C 2003, 28, 1117–1124. [Google Scholar] [CrossRef]
  5. Han, Z.; Dong, J.; Shen, Z.; Mou, R.; Zhou, Y.; Chen, X.; Fu, X.; Yang, C. Nitrogen removal of anaerobically digested swine wastewater by pilot-scale tidal flow constructed wetland based on in-situ biological regeneration of zeolite. Chemosphere 2019, 217, 364–373. [Google Scholar] [CrossRef]
  6. Strigul, N.S.; Kravchenko, L.V. Mathematical modeling of PGPR inoculation into the rhizosphere. Environ. Model. Softw. 2006, 21, 1158–1171. [Google Scholar] [CrossRef]
  7. Jia, L.; Gou, E.; Liu, H.; Lu, S.; Wu, S.; Wu, H. Exploring utilization of recycled agricultural biomass in constructed wetlands: Characterization of the driving force for high-rate nitrogen removal. Environ. Sci. Technol. 2019, 53, 1258–1268. [Google Scholar] [CrossRef]
  8. Choi, H.; Geronimo, F.K.F.; Jeon, M.; Kim, L.H. Investigation of the Factors Affecting the Treatment Performance of a Stormwater Horizontal Subsurface Flow Constructed Wetland Treating Road and Parking lot Runoff. Water 2021, 13, 1242. [Google Scholar] [CrossRef]
  9. Rousseau, D.P.L.; Louage, F.; Wang, Q.; Zhang, R. Constructed Wetlands for Urban Wastewater Treatment: An Overview. In Encyclopedia of Inland Waters, 2nd ed.; Mehner, T., Tockner, K., Eds.; Elsevier: Amsterdam, The Netherlands, 2022; pp. 272–284. [Google Scholar]
  10. Chen, Y.; Zhang, J.; Guo, Z.; Li, M.; Wu, H. Optimizing agricultural biomass application to enhance nitrogen removal in vertical flow constructed wetlands for treating low-carbon wastewater. Environ. Res. 2022, 209, 112867. [Google Scholar] [CrossRef]
  11. Zhang, Y.; Li, M.; Dong, L.; Han, C.; Li, M.; Wu, H. Effects of biochar dosage on treatment performance, enzyme activity and microbial community in aerated constructed wetlands for treating low C/N domestic sewage. Environ. Technol. Innovat. 2021, 24, 101919. [Google Scholar] [CrossRef]
  12. Jia, L.; Li, C.; Zhang, Y.; Chen, Y.; Li, M.; Wu, S.; Wu, H. Microbial community responses to agricultural biomass addition in aerated constructed wetlands treating low carbon wastewater. J. Environ. Manag. 2020, 270, 110912. [Google Scholar] [CrossRef]
  13. Wu, Y.; Chung, A.; Tam, N.F.Y.; Pi, N.; Wong, M.H. Constructed mangrove wetland as secondary treatment system for municipal wastewater. Ecol. Eng. 2008, 34, 137–146. [Google Scholar] [CrossRef]
  14. Leung, J.Y.S.; Cai, Q.; Tam, N.F.Y. Comparing subsurface flow constructed wetlands with mangrove plants and freshwater wetland plants for removing nutrients and toxic pollutants. Ecol. Eng. 2016, 95, 129–137. [Google Scholar] [CrossRef]
  15. Lobry, J.R.; Flandrois, J.P.; Carret, G.; Pave, A. Monod’s bacterial growth model revisited. Bull. Math. Biol. 1992, 54, 117–122. [Google Scholar] [CrossRef]
  16. Chiu, C.Y.; Lee, S.C.; Juang, H.T.; Hur, M.T.; Hwang, Y.H. Nitrogen nutritional status and fate of applied N in Mangrove soils. Bot. Bull. Acad. Sin. 1996, 37, 191–196. [Google Scholar]
  17. Pisman, T.I.; Pechurkin, N.S.; Mariasova, T.S.; Somova, L.A.; Sarangova, A.B. A mathematical model of “plants-microorganisms” interaction on complete mineral medium and under nitrogen limitation. Adv. Space Res. 1999, 24, 383–387. [Google Scholar] [CrossRef]
  18. Walker, R.L.; Burns, I.G.; Moorby, J. Responses of plant growth rate to nitrogen supply: A comparison of relative addition and N interruption treatments. J. Exp. Bot. 2001, 52, 309–317. [Google Scholar] [CrossRef]
  19. Pilyugin, S.S.; Waltman, P. Multiple limit cycles in the chemostat with variable yield. Math. Biosci. 2003, 182, 151–166. [Google Scholar] [CrossRef]
  20. Zhu, L.; Huang, X. Multiple limit cycles in a continuous culture vessel with variable yield. Nonlinear Anal. 2006, 64, 887–894. [Google Scholar] [CrossRef]
  21. Zhao, Z.; Pang, L.; Zhao, Z.; Luo, C. Impulsive State Feedback Control of the Rhizosphere Microbial Degradation in the Wetland Plant. Discret. Dyn. Nat. Soc. 2015, 2015, 612354. [Google Scholar] [CrossRef]
  22. Zhao, Z.; Song, Y.; Pang, L. Mathematical modeling of rhizosphere microbial degradation with impulsive diffusion on nutrient. Adv. Differ. Equ. 2016, 2016, 24. [Google Scholar] [CrossRef] [Green Version]
  23. Zhao, Z.; Li, Q.; Chen, L. Effect of rhizosphere dispersal and impulsive input on the growth of wetland plant. Math. Comput. Simul. 2018, 152, 69–80. [Google Scholar] [CrossRef]
  24. Zhao, Z.; Yin, Q.; Li, Q.; Wu, X. Mathematical model for diffusion of the rhizosphere microbial degradation with impulsive feedback control. J. Biol. Dyn. 2020, 14, 566–577. [Google Scholar] [CrossRef]
  25. Suandi, D.; Ningrum, I.P.; Alifah, A.N.; Izzah, N.; Reza, M.P.; Muwahidah, I.K. Mathematical Modeling and Sensitivity Analysis of the Existence of Male Calico Cats Population Based on Cross Breeding of All Coat Colour Types. Commun. Biomath. Sci. 2019, 2, 96–104. [Google Scholar] [CrossRef]
  26. Mapfumo, K.Z.; Pagan’a, J.C.; Juma, V.O.; Kavallaris, N.I.; Madzvamuse, A. A Model for the Proliferation–Quiescence Transition in Human Cells. Mathematics 2022, 10, 2426. [Google Scholar] [CrossRef]
  27. Rentzeperis, F.; Wallace, D. Local and global sensitivity analysis of spheroid and xenograft models of the acid-mediated development of tumor malignancy. Appl. Math. Model. 2022, 109, 629–650. [Google Scholar] [CrossRef]
  28. Alsahafi, S.; Woodcock, S. Exploring HIV Dynamics and an Optimal Control Strategy. Mathematics 2022, 10, 749. [Google Scholar] [CrossRef]
  29. Reinharz, V.; Churkin, A.; Dahari, H.; Barash, D. Advances in Parameter Estimation and Learning from Data for Mathematical Models of Hepatitis C Viral Kinetics. Mathematics 2022, 10, 2136. [Google Scholar] [CrossRef]
  30. Ansori, M.F.; Sumarti, N.; Sidarto, K.A.; Guandi, I. Analyzing a macroprudential instrument during the covid-19 pandemic using border collision bifurcation. Rev. Electron. Comun. Y Trab. ASEPUMA Rect. 2021, 22, 113–125. [Google Scholar] [CrossRef]
  31. Du, X.; Han, X.; Lei, C. Behavior Analysis of a Class of Discrete-Time Dynamical System with Capture Rate. Mathematics 2022, 10, 2410. [Google Scholar] [CrossRef]
  32. Saeed, T.; Djeddi, K.; Guirao, J.L.G.; Alsulami, H.H.; Alhodaly, M.S. A Discrete Dynamics Approach to a Tumor System. Mathematics 2022, 10, 1774. [Google Scholar] [CrossRef]
  33. Liu, B.; Wu, R. Bifurcation and Patterns Analysis for a Spatiotemporal Discrete Gierer-Meinhardt System. Mathematics 2022, 10, 243. [Google Scholar] [CrossRef]
  34. He, Z.Y.; Abbes, A.; Jahanshahi, H.; Alotaibi, N.D.; Wang, Y. Fractional-Order Discrete-Time SIR Epidemic Model with Vaccination: Chaos and Complexity. Mathematics 2022, 10, 165. [Google Scholar] [CrossRef]
  35. Gandolfo, G. Economic Dynamics: Methods and Models, 2nd ed.; Elsevier Science Publisher BV: Amsterdam, The Netherlands, 1985. [Google Scholar]
Figure 1. Diagram of total nitrogen transformation.
Figure 1. Diagram of total nitrogen transformation.
Symmetry 14 01924 g001
Figure 2. Phase portrait of the system (1) with some initial conditions. Panel (a) is for σ > σ ¯ , and Panel (b) is for σ < σ ¯ .
Figure 2. Phase portrait of the system (1) with some initial conditions. Panel (a) is for σ > σ ¯ , and Panel (b) is for σ < σ ¯ .
Symmetry 14 01924 g002
Figure 3. Various solutions of the system (1). (a) The solutions converge to E 1 * with σ > σ ¯ and γ > ϕ . (b) The solutions converge to E 1 * with σ > σ ¯ and γ < ϕ . (c) The solutions converge to E 3 * with σ < σ ¯ and γ > ϕ . (d) The solutions converge to E 3 * with σ < σ ¯ and γ < ϕ .
Figure 3. Various solutions of the system (1). (a) The solutions converge to E 1 * with σ > σ ¯ and γ > ϕ . (b) The solutions converge to E 1 * with σ > σ ¯ and γ < ϕ . (c) The solutions converge to E 3 * with σ < σ ¯ and γ > ϕ . (d) The solutions converge to E 3 * with σ < σ ¯ and γ < ϕ .
Symmetry 14 01924 g003
Figure 4. Sensitivity of the parameters to the variables. Panels (ah) are, respectively, for the cases of β , k, σ , Q, γ , α , a, and ϕ .
Figure 4. Sensitivity of the parameters to the variables. Panels (ah) are, respectively, for the cases of β , k, σ , Q, γ , α , a, and ϕ .
Symmetry 14 01924 g004
Figure 5. Sensitivity index of the parameters to the variables when arriving at the equilibrium.
Figure 5. Sensitivity index of the parameters to the variables when arriving at the equilibrium.
Symmetry 14 01924 g005
Figure 6. (a) Bifurcation diagram of parameter β . (b) Chaotic behavior of the solution of the system (10) when β = 0.8 .
Figure 6. (a) Bifurcation diagram of parameter β . (b) Chaotic behavior of the solution of the system (10) when β = 0.8 .
Symmetry 14 01924 g006
Figure 7. (a) Bifurcation diagram of parameter k. (b) Chaotic behavior of the solution of the system (10) when k = 0.15 .
Figure 7. (a) Bifurcation diagram of parameter k. (b) Chaotic behavior of the solution of the system (10) when k = 0.15 .
Symmetry 14 01924 g007
Figure 8. (a) Bifurcation diagram of parameter σ . (b) Chaotic behavior of the solution of the system (10) when σ = 0.03 .
Figure 8. (a) Bifurcation diagram of parameter σ . (b) Chaotic behavior of the solution of the system (10) when σ = 0.03 .
Symmetry 14 01924 g008
Figure 9. (a) Bifurcation diagram of parameter Q. (b) Chaotic behavior of the solution of the system (10) when Q = 5.3 .
Figure 9. (a) Bifurcation diagram of parameter Q. (b) Chaotic behavior of the solution of the system (10) when Q = 5.3 .
Symmetry 14 01924 g009
Figure 10. (a) Bifurcation diagram of parameter γ . (b) The solution’s behavior of the system (10) when γ = 1.83 .
Figure 10. (a) Bifurcation diagram of parameter γ . (b) The solution’s behavior of the system (10) when γ = 1.83 .
Symmetry 14 01924 g010
Figure 11. (a) Bifurcation diagram of parameter α . (b) The solution’s behavior of the system (10) when α = 1.72 .
Figure 11. (a) Bifurcation diagram of parameter α . (b) The solution’s behavior of the system (10) when α = 1.72 .
Symmetry 14 01924 g011
Figure 12. (a) Bifurcation diagram of parameter a. (b) The solution’s behavior of the system (10) when a = 0.2 .
Figure 12. (a) Bifurcation diagram of parameter a. (b) The solution’s behavior of the system (10) when a = 0.2 .
Symmetry 14 01924 g012
Figure 13. (a) Bifurcation diagram of parameter ϕ . (b) The solution’s behavior of the system (10) when ϕ = 1.96 .
Figure 13. (a) Bifurcation diagram of parameter ϕ . (b) The solution’s behavior of the system (10) when ϕ = 1.96 .
Symmetry 14 01924 g013
Table 1. Parameter values used in the simulations.
Table 1. Parameter values used in the simulations.
NotationDescriptionValueUnit
β Maximum possible value of the plant growth rate at infinite total nitrogen concentration in the soil solution0.25 time 1
kSemi-saturation1 [ M ]
σ Level of garbage0.1 time 1
QInput level of total nitrogen concentration1 [ M ] × time 1
γ Exchange rate of total nitrogen between wastewater and soil solution0.3 time 1
α Rate of total nitrogen loss in wastewater through runoff or evaporation0.2 time 1
aConversion of nutrients consumed from biomass produced0.5Dimensionless
ϕ Rate of total nitrogen loss in the soil solution by leaching or denitrification0.1 time 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sunarsih; Ansori, M.F.; Khabibah, S.; Sasongko, D.P. Continuous and Discrete Dynamical Models of Total Nitrogen Transformation in a Constructed Wetland: Sensitivity and Bifurcation Analysis. Symmetry 2022, 14, 1924. https://doi.org/10.3390/sym14091924

AMA Style

Sunarsih, Ansori MF, Khabibah S, Sasongko DP. Continuous and Discrete Dynamical Models of Total Nitrogen Transformation in a Constructed Wetland: Sensitivity and Bifurcation Analysis. Symmetry. 2022; 14(9):1924. https://doi.org/10.3390/sym14091924

Chicago/Turabian Style

Sunarsih, Moch. Fandi Ansori, Siti Khabibah, and Dwi Purwantoro Sasongko. 2022. "Continuous and Discrete Dynamical Models of Total Nitrogen Transformation in a Constructed Wetland: Sensitivity and Bifurcation Analysis" Symmetry 14, no. 9: 1924. https://doi.org/10.3390/sym14091924

APA Style

Sunarsih, Ansori, M. F., Khabibah, S., & Sasongko, D. P. (2022). Continuous and Discrete Dynamical Models of Total Nitrogen Transformation in a Constructed Wetland: Sensitivity and Bifurcation Analysis. Symmetry, 14(9), 1924. https://doi.org/10.3390/sym14091924

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop