Fuzzy Logic Theory and Spatiotemporal Modeling of the Fungus Phakopsora pachyrhizi Based on Differential Equations

Abstract

Brazil has been one of the largest soybean producers in recent years. The soybean is a legume commonly found in family meals. Among the diseases affecting the grains, Asian soybean rust is one of the most concerning. The fungus causing the disease is spread by the wind, making it difficult to control. Although it has been researched since its first records, not much data are available regarding the macro propagation behavior of spores. Therefore, this research aimed to model its dispersion based on a partial differential equation, the diffusion–advection equation, used by researchers to model the behavior of any pollutant. The terms of this equation were developed from real data, processed by fuzzy logic, and the simulation results were compared with disease records throughout a harvest. By using this approach to model the spatiotemporal dynamics of this fungus, it was possible to simulate its spread satisfactorily. Additionally, its results were used as input variables for a fuzzy system that estimates the susceptibility of a given location to disease development.

1. Introduction

Soy is a grain used in the preparation of food in many families. It is cultivated in all Brazilian regions, and thus, it has become one of the country’s most exported products in recent decades. Asian soybean rust is one of the diseases that most worry producers. If left unchecked, it can cause severe economic damage to legume crops.

Asian soybean rust is a plant disease that has been present in Brazilian agribusiness for 20 years [1]. It causes early grain ripening, reducing crop productivity significantly. It is caused by the fungus Phakopsora pachyrhizi (PP), which has soybeans as its principal host. The wind carries its spores, and volunteer soybeans on the sides of highways guarantee the disease’s continuity from one harvest to the next [2,3]. This research was dedicated to simulating the spatiotemporal dynamics of this fungus using a new approach.

Transport models have evolved for airborne diseases with a spore dispersal component based on atmospheric emissions models that calculate the dispersal of suspended aerosols emitted by a point source. These spores are treated as polluting particles, except for the issue of spore viability. Thus, long-distance dispersion, like any physical–dynamic problem in which both advection and diffusion are dominant, can be modeled using an Eulerian diffusion–advection approach [4]. Eulerian diffusion–advection models assume that atmospheric diffusion is analogous to molecular diffusion. Furthermore, they comply with Fick’s Law, which determines that the diffusion rate is proportional to the concentration gradient of the diffusing material [5].

A diffusion–advection equation was constructed to estimate the fungus concentration, and researchers used a differential equation to model the behavior of pollutants. From this equation, Missio (2008) [6] modeled the behavior of foot-and-mouth disease in cattle; Wolmuth (2009) and Krindges and Meyer (2011) [7,8] simulated the evolution of pollutants in Lago Manso, in Mato Grosso, in two and three dimensions; and Meyer and Diniz (2007) [9] modeled the dispersion of pollutants in an air–water system. The equation involves a diffusive term, which refers to the dispersion capacity of the pollutant; a transport term, represented by the vector field that transfers the pollutant from one location in the system to another; the decay, which can be interpreted as a mortality rate; and the term source, which guarantees the production of this pollutant.

The domain chosen for the construction of numerical approximations was the state of Mato Grosso, as it is the largest soybean producer in the country. The model was fed with data from the National Institute of Meteorology—INMET, and the results were compared with records of disease occurrences made available by the Antirust Consortium.

The text initially presents the problem’s classical and variational modeling and the development of spatial, temporal, and domain discretizations. It also explains how the terms diffusive, source, transport, and decay were constructed, constituting the equation. In sequence, it explains how the parameters used were defined and presents the results of the numerical simulations. Finally, the impressions left regarding the solutions obtained are recorded.

2. Materials and Methods

To develop this research, a partial differential equation was created that simulated the spread of the fungus PP in space and time during the 2018/2019 harvest over the state of Mato Grosso (Brazil). Mato Grosso is the largest soybean producer in the country. In the 2022/2023 harvest alone, it produced more than the entire South region (the second-largest grain-producing region in the country). In the 2023/2024 harvest, it was responsible for over a quarter of all production [10].

The diffusion–advection equation is frequently used to simulate pollutant dissemination (or concentration). Examples include a concentration of animals infected by a virus, river contaminants, or a specific poison’s dispersion in the air.

The equation was used to model the dispersion of pollutants in an air–water system. Initially, the air transport of the pollutant and its deposit on the surface of the water body were modeled with a two-dimensional vertical domain. Next, the aquatic environment was included for polluting substances that penetrate water, with a two-dimensional horizontal domain [9].

A system of partial differential equations of the susceptible, infected, and recovered (SIR) type was constructed, and the three equations were adjustments to the diffusion–advection equation. Thus, to simulate the behavior of foot-and-mouth disease in cattle, the first equation modeled the population susceptible to the disease; the second, the infected population; and the third, the recovered population [6].

The evolution of pollutants from agro-industrial activity in the Lago Manso region in Mato Grosso was also modeled using the diffusion–advection equation from two different approaches. In the first, the domain was two-dimensional, as it characterized the surface of a body of water. With the algorithmic tool built from this research, it was possible to generate scenarios using wind variations and reducing the flow of agrochemical products [7]. In the second, the domain covered was three-dimensional, including the lake’s depth. In addition to producing material and method for qualitative evaluation, with attention to this problem, the objective was to build a generic model that could be applied to any situation that fell into the same category as the second research [8].

For all these reasons, this equation was chosen to model the spatiotemporal dynamics of the fungus that causes Asian soybean rust. Next, the classic model of the diffusion–advection equation is presented, and a brief explanation of the specific terms is provided. Some concepts concerning variational formulation, discretization of equations, and domain are also addressed and contextualized to address the problem of this research.

2.1. Classical Model of the Equation

Calling $C(t,x,y)$ the concentration of the pollutant in the medium under study at the point $(x,y)$ in the plane and at the instant t, for the diffusion–advection equation for C, we obtain the following equation

$${\displaystyle \frac{\partial C}{\partial t}}+\left\{\mathrm{div}\left(\mathrm{flow}\right)\right\}+\left\{\mathrm{decay}\right\}=\left\{\mathrm{source}\right\}$$

(1)

More specifically, the pollutant flow (modeled in the equation by div(flow)) can be divided into two parts: diffusive flow (effective diffusion), which contains the effects of molecular and turbulent diffusion, and advective flow (transport), which is responsible for the macroscopic movements of the medium [7]:

$$\overrightarrow{J}=\overrightarrow{J_1}+\overrightarrow{J_2}$$

with $\overrightarrow{J}$ representing the flow, $\overrightarrow{J_1}$ the diffusive flow, and $\overrightarrow{J_2}$ the advective flow.

Making the separation between diffusive flow and advective flow, the Equation (1) can be rewritten as follows [7,11,12,13,14]:

$${\displaystyle \frac{\partial C}{\partial t}}=\left\{\mathrm{diffusion}\right\}-\left\{\mathrm{transport}\right\}-\left\{\mathrm{decay}\right\}+\left\{\mathrm{source}\right\}.$$

(2)

As the diffusive flow ($\overrightarrow{J_1}$) fulfills Fick’s law [15], it is possible to state that matter tends to spread, moving from places of greater concentration to areas of lower concentration. It is proportional to the pollutant density gradient ($\overrightarrow{J_1}=-\alpha \nabla C$). The advective flow ($\overrightarrow{J_2}$) is the movement caused by external agents, such as the medium’s velocity field. In this case, it is given by the product of velocity ($\overrightarrow{v}$) by concentration ($\overrightarrow{J_2}=\overrightarrow{v}C$) [7,16].

The system’s loss of mass, determined by the changes the substance undergoes over time due to the most diverse factors, determines the decay (or degradation). In many situations, we assume this loss is linearly proportional to the concentration. The means by which the substance is introduced (or removed) from the system are sources (or sinks) [7,16].

It is important to note that the gradient, divergence, and Laplacian operators in the following descriptions only relate to the spatial variables, not the temporal variable. Thus, for the classical modeling of these phenomena in a space of arbitrary size, we obtain the following equations [15,17,18]:

$$\begin{array}{ccc}\hfill \left\{\mathrm{diffusion}\right\}& =& \nabla ·\left(-\alpha \nabla C\right)\hfill \\ \hfill \left\{\mathrm{transport}\right\}& =& \nabla ·\left(\overrightarrow{v}C\right)\hfill \\ \hfill \left\{\mathrm{decay}\right\}& =& \sigma C\hfill \end{array}$$

Thus, the equation that models the phenomenon of dispersion for a given substance in a two-dimensional domain is given by

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial C}{\partial t}}=\nabla ·\left(\alpha \nabla C\right)-\nabla ·\left(\overrightarrow{v}C\right)-\sigma C+f,& (x,y)\in \mathrm{\Omega }\subset {\mathbb{R}}^2,& t\in (0,T]\subset \mathbb{R},\hfill \end{array}$$

being

• $\alpha =\alpha (t,x,y)$ the effective diffusion coefficient in the medium;
• $\overrightarrow{v}=\overrightarrow{v}(t,x,y)=(v_1(t,x,y),v_2(t,x,y))$ the velocity field in the middle;
• $\sigma =\sigma (t,x,y)$ the total decay coefficient in the medium;
• $f=f(t,x,y)$ is the source term.

The initial condition is given by $C(0,x,y)=C_0(x,y),(x,y)\in \mathrm{\Omega }\subset {\mathbb{R}}^2$. The following boundary conditions are considered, whose graphical representation are recorded in Figure 1:

1.

Non-homogeneous von Neumann-type for pollutant entry (in red)

$\alpha {\left.{\displaystyle \frac{\partial C}{\partial \eta }}\right|}_{{\mathrm{\Gamma }}_q}=\theta (t,x,y)$, being $\theta$ a given function and $q=1$;

2.

Homogeneous von Neumann-type for frontier without loss of pollutants (in blue)

$\alpha {\left.{\displaystyle \frac{\partial C}{\partial \eta }}\right|}_{{\mathrm{\Gamma }}_q}=0$, being $q=2$.

Thus, the classical model of the diffusion–advection equation involves a diffusive term, an advective transport, a decay coefficient, and a source term. The diffusive term consists of the divergence of the product between the diffusion coefficient and the concentration gradient. This coefficient can be constant or vary depending on time, space, environmental conditions, etc. Advective transport is given by the divergence of the scalar product between a vector field and the concentration. Decay is the product of the degradation coefficient and concentration. This coefficient can also be constant or variable. The term source, finally, characterizes the introduction of new objects and the concentration we want to model.

Finally, the equation for which the numerical solution is obtained is given by

$$\left\{\begin{array}{c}\frac{\partial C}{\partial t}=\nabla ·\left(\alpha \nabla C\right)-\nabla ·\left(\overrightarrow{v}C\right)-\sigma C+f,(x,y)\in \mathrm{\Omega }\subset {\mathbb{R}}^2,t\in (0,T]\subset \mathbb{R}\hfill \\ \nabla ·\overrightarrow{v}=0\hfill \\ \alpha {\left.{\displaystyle \frac{\partial C}{\partial \eta }}\right|}_{{\mathrm{\Gamma }}_1}=\theta (t,x,y)\hfill \\ \alpha {\left.{\displaystyle \frac{\partial C}{\partial \eta }}\right|}_{{\mathrm{\Gamma }}_2}=0\hfill \\ C(0,x,y)=C_0(x,y),\forall (x,y)\in \mathrm{\Omega }\hfill \end{array}\right.$$

(3)

It is possible to affirm that it is a second-order linear differential equation. Thus, for Equation (3) to be satisfied, it is necessary to guarantee the continuity of the solution function C and its partial derivatives up to the second-order level. However, this is only sometimes possible. For these cases, it is necessary to construct a variational formulation, also called a weak formulation of the differential equation.

2.2. Variational Formulation

Although differential equations mathematically model several problems, many need analytical solutions. In these cases, the variational formulation of the equation is used, among other techniques, for obtaining numerical solutions. This reformulation makes it possible to obtain approximate values of the equation at pre-determined points in its domain, enabling a general understanding of the phenomenon studied.

The variational formulation rewrites the classical formulation of the differential equation (or system of equations), so its solutions must be sought in a space with fewer differentiability restrictions. These solutions are called weak solutions [6]. For linear differential equations, it is possible to state that the solution to an initial value problem exists in any interval that contains the starting point and in which the coefficients of the differential equation are continuous. If any of the coefficients is discontinuous at specific points, the solution may be discontinuous or even cease to exist at those points in the interval [19].

Thus, to satisfy the differential equation in its classical model (3), the solution function C must be, obligatorily, at least of class $C^2$; in addition to being continuous throughout the domain, its first and second derivatives must also be continuous. In practice, however, the function that models the source term or initial conditions is only sometimes continuous concerning time or space. In other words, environmental problems generally do not have the regularity conditions that the classical formulation requires [7,20].

Approximation methods, whose instruments do not have the regularity requirements of the robust formulation, are an alternative to these problems [7]. Thus, the natural option is to use the variational formulation (also called weak formulation of the problem), which “weakens” the hypotheses of the solution’s regularity, “widening” the conditions necessary to find or construct the sought solution [14].

In variational methods, the objective is to obtain an approximate solution of a differential equation within a convenient space of functions of finite dimension. In this subspace, the approximate solution is an exact projection of the solution. It can be written as a linear combination of a finite number of known functions, which form the basis of this subspace. Thus, with the resolution of the system of algebraic equations resulting from the discretization process, the coefficients of this linear expansion are sought to satisfy the differential equation. Thus, the problem must be rewritten in its weak form [21]. However, before obtaining the variational formulation (or weak formulation) of equation (3), it is necessary to understand some important definitions and state Green’s theorem, which is presented below.

Assuming $u\in \mathrm{\Omega }$ and $v\in \mathrm{\Omega }$ functions, with $\mathrm{\Omega }\subset {\mathbb{R}}^2$. If

• u and its first and second derivatives are continuous;
• v and its first derivatives are continuous.

Green’s theorem guarantees that the following equality holds [22]:

$${\int }_{\mathrm{\Omega }}(-\Delta u)vdx={\int }_{\mathrm{\Omega }}\nabla u·\nabla vdx-{\int }_{\partial \mathrm{\Omega }}(\eta ·\nabla u)vd\partial \mathrm{\Omega }.$$

where $\eta$ is the unit vector out of $\partial \mathrm{\Omega }$.

Let $L^2\left(\mathrm{\Omega }\right)$ be the space of functions $u\in \mathrm{\Omega }$, which are Lebesgue integrable, over $\mathrm{\Omega }$, with $\mathrm{\Omega }\subset {\mathbb{R}}^2$. We have

$$L^2\left(\mathrm{\Omega }\right)=\left\{u:\mathrm{\Omega }\to \mathbb{R},\mathrm{with}{\int }_{\mathrm{\Omega }}{\left|u(x,y)\right|}^{2}d\mu <\infty \right\}.$$

It considers $u\in L^2\left(\mathrm{\Omega }\right)$, so u has weak derivatives of all orders. It is not possible to say that these weak derivatives also belong to $L^2\left(\mathrm{\Omega }\right)$. Thus, it is necessary to define the space $H^m\left(\mathrm{\Omega }\right)$, so that if the function $u(x,y)$ belongs to $H^m\left(\mathrm{\Omega }\right)$, then $u(x,y)$ and its derivatives ${\displaystyle \frac{{\partial }^{i}u}{\partial x^i}}$ and ${\displaystyle \frac{{\partial }^{i}u}{\partial y^i}}$, with $1\le i\le m$, belong to $L^2\left(\mathrm{\Omega }\right)$ [14]. For the weak formulation of the Equation (3), it is enough to guarantee that the first partial derivatives belong to $L^2\left(\mathrm{\Omega }\right)$. Thus, the Sobolev space is defined:

$$H^1\left(\mathrm{\Omega }\right)=\left\{u(x,y)\in L^2\left(\mathrm{\Omega }\right):{\displaystyle \frac{\partial u}{\partial x}}\in L^2\left(\mathrm{\Omega }\right)\mathrm{and}{\displaystyle \frac{\partial u}{\partial y}}\in L^2\left(\mathrm{\Omega }\right)\right\}.$$

Internal products and norms are defined for the spaces $L^2\left(\mathrm{\Omega }\right)$ and $H^1\left(\mathrm{\Omega }\right)$ to facilitate the search for test functions. From these, an approximate solution to the problem is constructed. For the development of this variational formulation, the following are established as the inner product and norm of $L^2\left(\mathrm{\Omega }\right)$, respectively:

$${(u,v)}_{L^2\left(\mathrm{\Omega }\right)}={(u,v)}_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }}uv\mathrm{d}\mu \mathrm{and} {\parallel u\parallel }_{L^2\left(\mathrm{\Omega }\right)}^2={\parallel u\parallel }_{\mathrm{\Omega }}^2={(u,u)}_{\mathrm{\Omega }}$$

over the domain $\mathrm{\Omega }$, with $\mu$ as the Lebesgue integral measure. Furthermore,

$${\langle u,v\rangle }_{L^2\left(\mathrm{\Gamma }\right)}={\langle u,v\rangle }_{\mathrm{\Gamma }}={\int }_{\mathrm{\Gamma }}uv\mathrm{d}\partial \mathrm{\Omega }$$

the inner product over part of the frontier of $\mathrm{\Omega }$ being $\mathrm{\Gamma }\subset \partial \mathrm{\Omega }$.

In $H^1\left(\mathrm{\Omega }\right)\subset L^2\left(\mathrm{\Omega }\right)$, we have as inner product and norm, respectively,

$${(u,v)}_{H^1\left(\mathrm{\Omega }\right)}={(u,v)}_{\mathrm{\Omega }}+{(\nabla u\parallel \nabla v)}_{\mathrm{\Omega }} \mathrm{and} {\parallel u\parallel }_{H^1\left(\mathrm{\Omega }\right)}^2={\parallel u\parallel }_{\mathrm{\Omega }}^2+{\parallel \nabla u\parallel }_{\mathrm{\Omega }}^2$$

with ${(\nabla u\parallel \nabla v)}_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }}{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial x}}$d$\mu +{\int }_{\mathrm{\Omega }}{\displaystyle \frac{\partial u}{\partial y}}{\displaystyle \frac{\partial v}{\partial y}}$d$\mu$, $u\in H^1\left(\mathrm{\Omega }\right)$, and $v\in H^1\left(\mathrm{\Omega }\right)$.

That said, working with the variational formulation consists of considering the derivatives of the classical model in the weak sense and taking the inner product of each term of (3) using the test function $w\in H^1\left(\mathrm{\Omega }\right)$ [6,7,14].

Multiplying each term of (3) by $w\in H^1\left(\mathrm{\Omega }\right)$, we obtain

$${\displaystyle \frac{\partial C}{\partial t}}w=\nabla ·\left(\alpha \nabla C\right)w-\nabla ·\left(\overrightarrow{v}C\right)w-\left(\sigma C\right)w+fw.$$

(4)

Integrating (4) concerning spatial variables,

$${\int }_{\mathrm{\Omega }}{\displaystyle \frac{\partial C}{\partial t}}wd\mu ={\int }_{\mathrm{\Omega }}\nabla ·\left(\alpha \nabla C\right)wd\mu -{\int }_{\mathrm{\Omega }}\nabla ·\left(\overrightarrow{v}C\right)wd\mu -{\int }_{\mathrm{\Omega }}\left(\sigma C\right)wd\mu +{\int }_{\mathrm{\Omega }}fwd\mu$$

(5)

for all $w\in H^1\left(\mathrm{\Omega }\right)$.

Just as in Diniz (2003) [12], which considered $\nabla ·\overrightarrow{v}=0$ to approximate a “well-behaved” field in the direction of the air flow, in this research, we also adopted $\nabla ·\overrightarrow{v}=0$. Therefore, using $\nabla ·\left(\overrightarrow{v}C\right)=\left(\nabla ·\overrightarrow{v}\right)·C+\overrightarrow{v}·\nabla C$, we have $\nabla ·\left(\overrightarrow{v}C\right)=\overrightarrow{v}·\nabla C$. Thus, substituting in (5) and isolating the source term, we have

$${\int }_{\mathrm{\Omega }}{\displaystyle \frac{\partial C}{\partial t}}wd\mu -\alpha {\int }_{\mathrm{\Omega }}\nabla ·\left(\nabla C\right)wd\mu +{\int }_{\mathrm{\Omega }}(\overrightarrow{v}·\nabla C)wd\mu +{\int }_{\mathrm{\Omega }}\left(\sigma C\right)wd\mu ={\int }_{\mathrm{\Omega }}fwd\mu$$

(6)

for all $w\in H^1\left(\mathrm{\Omega }\right)$.

Using Green’s theorem for the diffusive term, we have

$$\alpha {\int }_{\mathrm{\Omega }}\nabla ·(\nabla C)wd\mu =-\alpha {\int }_{\mathrm{\Omega }}(\nabla C·\nabla w)d\mu +\alpha {\int }_{\partial \mathrm{\Omega }}{\displaystyle \frac{\partial C}{\partial \eta }}wd\partial \mathrm{\Omega }.$$

(7)

Knowing that $\partial \mathrm{\Omega }={\bigcup }_{q=1}^2{\mathrm{\Gamma }}_q$ and given the boundary conditions, we have

$$\alpha {\int }_{\partial \mathrm{\Omega }}{\displaystyle \frac{\partial C}{\partial \eta }}wd\partial \mathrm{\Omega }={\int }_{{\mathrm{\Gamma }}_1}\theta wd\partial \mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_2}0wd\partial \mathrm{\Omega }.$$

(8)

In this way, replacing (8) in (7), the diffusive term becomes

$$\alpha {\int }_{\mathrm{\Omega }}\nabla ·(\nabla C)wd\mu =-\alpha {\int }_{\mathrm{\Omega }}(\nabla C·\nabla w)d\mu +\left({\int }_{{\mathrm{\Gamma }}_1}\theta wd\partial \mathrm{\Omega }\right).$$

(9)

Using (9) in (6), we have

$${\int }_{\mathrm{\Omega }}{\displaystyle \frac{\partial C}{\partial t}}wd\mu +\alpha {\int }_{\mathrm{\Omega }}(\nabla C·\nabla w)d\mu +{\int }_{\mathrm{\Omega }}(\overrightarrow{v}·\nabla C)wd\mu +{\int }_{\mathrm{\Omega }}\left(\sigma C\right)wd\mu ={\int }_{\mathrm{\Omega }}fwd\mu +{\int }_{{\mathrm{\Gamma }}_1}\theta wd\partial \mathrm{\Omega }$$

(10)

for all $w\in H^1\left(\mathrm{\Omega }\right)$.

Rewriting the equation using the inner product notations introduced previously, we have

$${\left({\displaystyle \frac{\partial C}{\partial t}},w\right)}_{\mathrm{\Omega }}+\alpha {\left(\nabla C\parallel \nabla w\right)}_{\mathrm{\Omega }}+{\left((\overrightarrow{v}·\nabla C),w\right)}_{\mathrm{\Omega }}+{\left(\left(\sigma C\right),w\right)}_{\mathrm{\Omega }}={\left(f,w\right)}_{\mathrm{\Omega }}+{\langle \theta ,w\rangle }_{{\mathrm{\Gamma }}_1}$$

(11)

for all $w\in H^1\left(\mathrm{\Omega }\right)$. The Equation (11) is, therefore, the variational formulation of the Equation (3).

Thus, a solution $C\in \mathcal{V}$ that satisfies (11) is sought, with $\mathcal{V}$ being the space of functions that belong to $L^2\left(\mathrm{\Omega }\right)$ and whose spatial (in x and y) and temporal derivatives also belong to $L^2\left(\mathrm{\Omega }\right)$:

$$\mathcal{V}=\left\{C\in L^2[(0,T),H^1\left(\mathrm{\Omega }\right)]:{\displaystyle \frac{\partial C}{\partial t}}\in L^2\left(\mathrm{\Omega }\right),\forall t\in [0,T]\right\}.$$

Based on the demonstrations developed in his doctoral thesis, Inforzato (2008) [13] proved the existence of a unique solution to his equation. Making the necessary adjustments guarantees the existence and uniqueness of the solution to the Equation (11). Once the variational formulation (11) has been constructed, some suitable numerical method must approximate the solution. The Galerkin finite element method [11] uses the variational formulation of the equations, allowing the discretization of spatial variables and obtaining numerical approximations of the solution.

2.2.1. Spatial Discretization

Let ${\mathcal{V}}^{\left(h\right)}$ be a finite dimensional subspace of $H^1\left(\mathrm{\Omega }\right)$ generated by the basis $\mathcal{B}=\{{\phi }_1,{\phi }_2,\dots ,{\phi }_{n_h}\}$. The approximation of the solution of (11), or its projection onto the subspace ${\mathcal{V}}^{\left(h\right)}$ can be written as a linear combination of the elements of the base $\mathcal{B}$:

$$C_h=\sum _{j=1}^{n_h}{c}_j\left(t\right){\phi }_j(x,y).$$

(12)

These functions (${\phi }_i$), called base functions, allow the construction of an approximation of the solution to the variational problem. In this case, the problem domain is a finite spatial discretization on which the basis $\mathcal{B}$ of functions of the finite subspace ${\mathcal{V}}^{\left(h\right)}$ is produced [16].

Replacing $C_h$ in the variational formulation (11), we have

$${\left({\displaystyle \frac{\partial C_h}{\partial t}},w\right)}_{\mathrm{\Omega }}+\alpha {\left(\nabla C_h\parallel \nabla w\right)}_{\mathrm{\Omega }}+{\left((\overrightarrow{v}·\nabla C_h),w\right)}_{\mathrm{\Omega }}+{\left(\left(\sigma C_h\right),w\right)}_{\mathrm{\Omega }}={\left(f,w\right)}_{\mathrm{\Omega }}+{\langle \theta ,w\rangle }_{{\mathrm{\Gamma }}_1}$$

(13)

for all $w\in {\mathcal{V}}^{\left(h\right)}$. By using (12), we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial C_h}{\partial t}}=\sum _{j=1}^{n_h}{\displaystyle \frac{d{c}_j}{dt}}{\phi }_j& \hspace{1em}\mathrm{and}\hspace{1em}& \nabla C_h=\sum _{j=1}^{n_h}{c}_j\nabla {\phi }_j.\hfill \end{array}$$

(14)

Using (12) and (14) in (13), we have

$${\left(\sum _{j=1}^{n_h}{\displaystyle \frac{d{c}_j}{dt}}{\phi }_j,w\right)}_{\mathrm{\Omega }}+\alpha {\left(\sum _{j=1}^{n_h}{c}_j\nabla {\phi }_j\parallel \nabla w\right)}_{\mathrm{\Omega }}+{\left(\left(\overrightarrow{v}·\sum _{j=1}^{n_h}{c}_j\nabla {\phi }_j\right),w\right)}_{\mathrm{\Omega }}+{\left(\left(\sigma \sum _{j=1}^{n_h}{c}_j{\phi }_j\right),w\right)}_{\mathrm{\Omega }}={\left(f,w\right)}_{\mathrm{\Omega }}+{\langle \theta ,w\rangle }_{{\mathrm{\Gamma }}_1}$$

(15)

for all $w\in {\mathcal{V}}^{\left(h\right)}$.

This procedure corresponds to the concept of separation of variables, and as the functions $c_j\left(t\right)$ do not depend on $(x,y)$, it is possible to remove them from the inner product [16]. Then,

$$\sum _{j=1}^{n_h}{\displaystyle \frac{d{c}_j}{dt}}{\left({\phi }_j,w\right)}_{\mathrm{\Omega }}+\alpha \sum _{j=1}^{n_h}{c}_j{\left(\nabla {\phi }_j\parallel \nabla w\right)}_{\mathrm{\Omega }}+\sum _{j=1}^{n_h}{c}_j{\left((\overrightarrow{v}·\nabla {\phi }_j),w\right)}_{\mathrm{\Omega }}+\sum _{j=1}^{n_h}{c}_j{\left(\left(\sigma {\phi }_j\right),w\right)}_{\mathrm{\Omega }}={\left(f,w\right)}_{\mathrm{\Omega }}+{\langle \theta ,w\rangle }_{{\mathrm{\Gamma }}_1}$$

(16)

for all $w\in {\mathcal{V}}^{\left(h\right)}$.

In the Galerkin Method, the test functions (also called weight functions) w are of the same class as the base functions [16]. In this way, as the Equation (16) is valid for all $w\in {\mathcal{V}}^{\left(h\right)}$, it is sufficient to calculate it for the elements in the base of ${\mathcal{V}}^{\left(h\right)}$. Therefore,

$$\sum _{j=1}^{n_h}{\displaystyle \frac{d{c}_j}{dt}}{\left({\phi }_j,{\phi }_i\right)}_{\mathrm{\Omega }}+\alpha \sum _{j=1}^{n_h}{c}_j{\left(\nabla {\phi }_j\parallel \nabla {\phi }_i\right)}_{\mathrm{\Omega }}+\sum _{j=1}^{n_h}{c}_j{\left((\overrightarrow{v}·\nabla {\phi }_j),{\phi }_i\right)}_{\mathrm{\Omega }}+\sum _{j=1}^{n_h}{c}_j{\left(\left(\sigma {\phi }_j\right),{\phi }_i\right)}_{\mathrm{\Omega }}={\left(f,{\phi }_i\right)}_{\mathrm{\Omega }}+{\langle \theta ,{\phi }_i\rangle }_{{\mathrm{\Gamma }}_1}$$

(17)

for all ${\phi }_i\in \mathcal{B}=\{{\phi }_1,{\phi }_2,\dots ,{\phi }_{n_h}\}$.

Organizing the Equation (17) possibility to write this sentence in the form of a system of first-order ordinary differential equations, we obtain

$$\sum _{j=1}^{n_h}{\displaystyle \frac{d{c}_j}{dt}}{\left({\phi }_j,{\phi }_i\right)}_{\mathrm{\Omega }}+\sum _{j=1}^{n_h}{c}_j\left[\alpha {\left(\nabla {\phi }_j\parallel \nabla {\phi }_i\right)}_{\mathrm{\Omega }}+{\left((\overrightarrow{v}·\nabla {\phi }_j),{\phi }_i\right)}_{\mathrm{\Omega }}+{\left(\left(\sigma {\phi }_j\right),{\phi }_i\right)}_{\mathrm{\Omega }}\right]={\left(f,{\phi }_i\right)}_{\mathrm{\Omega }}+{\langle \theta ,{\phi }_i\rangle }_{{\mathrm{\Gamma }}_1}$$

(18)

for all ${\phi }_i\in \mathcal{B}$.

Making $C_h={\left[c_1\left(t\right),c_2\left(t\right),\dots ,c_{n_h}\left(t\right)\right]}^T$ and $C_h^{\prime }={\left[{\displaystyle \frac{d{c}_1}{dt}},{\displaystyle \frac{d{c}_2}{dt}},\dots ,{\displaystyle \frac{d{c}_{n_h}}{dt}}\right]}^T$ from this expression, we have

$$M{C}_h^{\prime }\left(t\right)+W{C}_h\left(t\right)=H+N,$$

(19)

the matrices being $M={\left(m_{ij}\right)}_{n_h\times n_h}$, $W={\left(w_{ij}\right)}_{n_h\times n_h}$, $H={\left(h_i\right)}_{n_h\times 1}$, and $N={\left(n_i\right)}_{n_h\times 1}$ given by:

$$\begin{array}{c}{m}_{ij}={({\phi }_j,{\phi }_i)}_{\mathrm{\Omega }};\hfill \\ w_{ij}=\alpha {\left(\nabla {\phi }_j\parallel \nabla {\phi }_i\right)}_{\mathrm{\Omega }}+{\left((\overrightarrow{v}·\nabla {\phi }_j),{\phi }_i\right)}_{\mathrm{\Omega }}+{\left(\left(\sigma {\phi }_j\right),{\phi }_i\right)}_{\mathrm{\Omega }};\hfill \\ h_i={\left(f,{\phi }_i\right)}_{\mathrm{\Omega }};\hfill \\ n_i={\langle \theta ,{\phi }_i\rangle }_{{\mathrm{\Gamma }}_1}.\hfill \end{array}$$

(20)

The vectors $C_h^{\prime }\left(t\right)$ and $C_h\left(t\right)$ are obtained from the system resolution. Thus, we have a system of first-order ordinary differential equations for which solutions are obtained with relative theoretical ease.

2.2.2. Temporal Discretization

The temporal discretization of the equation occurs using the finite difference scheme known as Crank–Nicolson [9,11]. This method gives rise to an unconditionally stable scheme of order two and has been widely used. It consists of taking the average instant between two subsequent instants; that is, using the time-centered difference $t^{n+{\displaystyle \frac{1}{2}}}$, obtaining the following approximations:

$$\begin{array}{cccc}& {\left.{\displaystyle \frac{d{c}_j}{dt}}\right|}_{t^{n+{\displaystyle \frac{1}{2}}}}={\displaystyle \frac{c_j^{n+1}-c_j^n}{\Delta t}},& c_j^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{c_j^{n+1}+c_j^n}{2}},& \\ M^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{M^{n+1}+M^n}{2}},& W^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{W^{n+1}+W^n}{2}},& H^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{H^{n+1}+H^n}{2}},& N^{n+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{N^{n+1}+N^n}{2}}.\end{array}$$

Replacing t in (19) by $t^{n+{\displaystyle \frac{1}{2}}}$, we have

$$M^{n+{\displaystyle \frac{1}{2}}}\left[{\displaystyle \frac{d{C}_h^{n+\frac{1}{2}}}{dt}}\right]+W^{n+{\displaystyle \frac{1}{2}}}\left[C_h^{n+{\displaystyle \frac{1}{2}}}\right]=H^{n+{\displaystyle \frac{1}{2}}}+N^{n+{\displaystyle \frac{1}{2}}}.$$

Thus,

$$M^{n+{\displaystyle \frac{1}{2}}}\left[{\displaystyle \frac{C_h^{n+1}-C_h^n}{\Delta t}}\right]+W^{n+{\displaystyle \frac{1}{2}}}\left[{\displaystyle \frac{C_h^{n+1}+C_h^n}{2}}\right]=H^{n+{\displaystyle \frac{1}{2}}}+N^{n+{\displaystyle \frac{1}{2}}}.$$

Factoring the equation by taking out $C_h^{n+1}$ and $C_h^n$ and multiplying everything by $\Delta t$, the solution to the system of differential equations (19) can be approximated by the solution of the following system:

$$\left[M^{n+{\displaystyle \frac{1}{2}}}+{\displaystyle \frac{\Delta t}{2}}{W}^{n+{\displaystyle \frac{1}{2}}}\right]·C_h^{n+1}=\left[M^{n+{\displaystyle \frac{1}{2}}}-{\displaystyle \frac{\Delta t}{2}}{W}^{n+{\displaystyle \frac{1}{2}}}\right]·C_h^n+\Delta t\left(H^{n+{\displaystyle \frac{1}{2}}}+N^{n+{\displaystyle \frac{1}{2}}}\right).$$

(21)

The system (21) is solved iteratively, in time, starting from an initial condition. Furthermore, some inner products must be calculated to construct the M, W, H, and N matrices, as defined in (20). To this end, the $\phi$ functions used in the process are defined below.

2.2.3. Domain Discretization

Discretizing the domain involves replacing the analysis in a continuous region with the analysis in a finite set of points. A standard method used for this is finite elements. Thus, over the domain of the equation, a mesh is built that is formed by geometric elements, which, in the case of the two-dimensional domain, can be triangles and quadrilaterals, among others. From this mesh and the appropriate choice of functions that make up the basis of the space ${\mathcal{V}}^{\left(h\right)}$, it is possible to obtain spatial approximations for the solution via finite elements.

The finite element technique is suitable for composing approximations to solve boundary value problems. It consists of dividing the solution domain into a finite number of simple subdomains, also called simple elements, and taking convenient functions (test functions) in each. This way, an approximate partial differential equation solution can be obtained locally for each element. All these domains and their functions are strategically brought together, thus generating the total solution [7].

The method consists of appropriately choosing the $\phi$ functions so that, together with their derivatives, they are non-zero in a small portion of the domain. This reduces the domain’s integrals to integrals in part of the domain, where the functions and their derivatives are simultaneously different from zero, which simplifies the integration process [21].

The procedure boils down to the following:

1.

Create a mesh in $\mathrm{\Omega }$, defining which geometric figure is used for the finite elements;

2.

Determine the degree of the approximation polynomials (which act on each element).

Being appropriate to the problem, the base functions $\mathcal{B}=\{{\phi }_1,{\phi }_2,\dots ,{\phi }_{n_h}\}$ can be chosen among the types piecewise linear, satisfying the following condition:

$${\phi }_i(x_j,y_j)=\left\{\begin{array}{cc}\hfill 1,& \mathrm{se}i=j\hfill \\ \hfill 0,& \mathrm{se}i\ne j\hfill \end{array}\right.,$$

(22)

with $(x_i,y_i)$ being the coordinates of the vertices of the geometric figures, also called mesh nodes [20].

To solve the system (21), the discretization of the region $\mathrm{\Omega }$ occurs using first-order triangular finite elements. The base functions are defined to satisfy the condition (22). Denote, with ${\left\{k_e\right\}}_{e=1}^{nte},$ a finite family of $nte$ triangles $k_e$ that are two-by-two disjoint or have as intersection only one edge or a vertex, such that

$${\mathrm{\Omega }}_h=\bigcup _{e=1}^{nte}{k}_e.$$

Each mesh triangle is a simple subdomain, where the integrals are calculated locally. Set the triangle whose vertices are $(1,0)$, $(0,1)$, and $(0,0)$ as reference ($\widehat{k}$). This way, the functions $\phi$ restricted to this reference element are

$${\phi }_{\widehat{1}}(\xi ,\eta )=\xi ;$$

(23)

$${\phi }_{\widehat{2}}(\xi ,\eta )=\eta ;$$

(24)

$${\phi }_{\widehat{3}}(\xi ,\eta )=1-\xi -\eta .$$

(25)

The technique used to implement the finite element method requires the definition of an affine transformation $T_e:\widehat{k}\to k_e$ that relates the reference element $\widehat{k}$ to a generic element $k_e$. That is, the way (23), (24) and (25) were defined, $T_e$ must transform

$$\begin{array}{c}(1,0)\to (x_1,y_1);\\ (0,1)\to (x_2,y_2);\\ (0,0)\to (x_3,y_3).\end{array}$$

Thus, it is possible to define the following linear transformation:

$$T_e(\xi ,\eta )=\left(\begin{array}{cc}{x}_1-x_3& x_2-x_3\\ y_1-y_3& y_2-y_3\end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)+\left(\begin{array}{c}{x}_3\\ y_3\end{array}\right)$$

(26)

which is equivalent to

$$\begin{array}{c}x(\xi ,\eta )=x_3+(x_1-x_3)\xi +(x_2-x_3)\eta \mathrm{and}\\ y(\xi ,\eta )=y_3+(y_1-y_3)\xi +(y_2-y_3)\eta \end{array}$$

and whose Jacobian matrix is given by

$$J{T}_e(\xi ,\eta )=\left(\begin{array}{cc}{x}_1-x_3& x_2-x_3\\ y_1-y_3& y_2-y_3\end{array}\right).$$

Knowing that ${\phi }_{\widehat{i}}(\xi ,\eta )={\phi }_i(x(\xi ,\eta ),y(\xi ,\eta ))$, the derivatives of the functions global values in terms of local derivatives are given by:

$$\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}={\displaystyle \frac{\partial {\phi }_i}{\partial x}}{\displaystyle \frac{\partial x}{\partial \xi }}+{\displaystyle \frac{\partial {\phi }_i}{\partial y}}{\displaystyle \frac{\partial y}{\partial \xi }}={\displaystyle \frac{\partial {\phi }_i}{\partial x}}(x_1-x_3)+{\displaystyle \frac{\partial {\phi }_i}{\partial y}}(y_1-y_3) \mathrm{and}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}={\displaystyle \frac{\partial {\phi }_i}{\partial x}}{\displaystyle \frac{\partial x}{\partial \eta }}+{\displaystyle \frac{\partial {\phi }_i}{\partial y}}{\displaystyle \frac{\partial y}{\partial \eta }}={\displaystyle \frac{\partial {\phi }_i}{\partial x}}(x_2-x_3)+{\displaystyle \frac{\partial {\phi }_i}{\partial y}}(y_2-y_3).\end{array}$$

Matrix-wise, we have

$$\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)=\left(\begin{array}{cc}{x}_1-x_3& y_1-y_3\\ x_2-x_3& y_2-y_3\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\\ {\displaystyle \frac{\partial {\phi }_i}{\partial y}}\end{array}\right)={\left(J{T}_e\right)}^T\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\\ {\displaystyle \frac{\partial {\phi }_i}{\partial y}}\end{array}\right).$$

Calling $T={\left({\left(J{T}_e\right)}^T\right)}^{-1}$, we have

$$\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\\ {\displaystyle \frac{\partial {\phi }_i}{\partial y}}\end{array}\right)=T\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right).$$

Using this relationship, it is possible to write

$$\begin{array}{c}{\displaystyle \frac{\partial {\phi }_i}{\partial x}}=\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)=T(1,:)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)\\ {\displaystyle \frac{\partial {\phi }_i}{\partial y}}=\left(T_{21}\hspace{1em}{T}_{22}\right)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)=T(2,:)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)\end{array}$$

with $T(1,:)$ and $T(2,:)$ being the first and second lines of T, respectively.

With the derivatives of the global basis functions written as a linear combination of the local basis functions, it is possible to write the integrals over the global elements as a function of the integrals over the reference element. Calculations are divided into cases.

First case:

$$\begin{array}{c}\hfill {\left({\phi }_j,{\phi }_i\right)}_{k_e}={\int \int }_{k_e}{\phi }_j{\phi }_{i}dydx={\int \int }_{\widehat{k}}{\phi }_{\widehat{j}}{\phi }_{\widehat{i}}|detJ{T}_e|d\eta d\xi =|detJ{T}_e|{\int \int }_{\widehat{k}}{\phi }_{\widehat{j}}{\phi }_{\widehat{i}}d\eta d\xi =|detJ{T}_e|{({\phi }_{\widehat{j}},{\phi }_{\widehat{i}})}_{\widehat{k}}.\end{array}$$

Second case:

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{\partial {\phi }_j}{\partial x}},{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\right)}_{k_e}& ={\int \int }_{k_e}{\displaystyle \frac{\partial {\phi }_j}{\partial x}}{\displaystyle \frac{\partial {\phi }_i}{\partial x}}dydx\hfill \\ & ={\int \int }_{\widehat{k}}\left[\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{c}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}\end{array}\right)\right]\left[\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{c}\frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }\\ \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }\end{array}\right)\right]|detJ{T}_e|d\eta d\xi \hfill \\ & =|detJ{T}_e|{\int \int }_{\widehat{k}}\left[\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{c}\frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }\\ \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }\end{array}\right)\right]\left[\left(\begin{array}{cc}\frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }& \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }\end{array}\right)\left(\begin{array}{c}{T}_{11}\\ T_{12}\end{array}\right)\right]d\eta d\xi \hfill \\ & =|detJ{T}_e|{\int \int }_{\widehat{k}}\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{cc}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}& {\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\\ {\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}& {\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\end{array}\right)\left(\begin{array}{c}{T}_{11}\\ T_{12}\end{array}\right)d\eta d\xi \hfill \\ & =|detJ{T}_e|\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{cc}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\right)}_{\widehat{k}}& {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\right)}_{\widehat{k}}\\ {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\right)}_{\widehat{k}}& {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\right)}_{\widehat{k}}\end{array}\right)\left(\begin{array}{c}{T}_{11}\\ T_{12}\end{array}\right).\hfill \end{array}$$

Calling $Sm=\left(\begin{array}{cc}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\right)}_{\widehat{k}}& {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\right)}_{\widehat{k}}\\ {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }}\right)}_{\widehat{k}}& {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }}\right)}_{\widehat{k}}\end{array}\right)$, and repeating the process to ${\left({\displaystyle \frac{\partial {\phi }_j}{\partial x}},{\displaystyle \frac{\partial {\phi }_i}{\partial y}}\right)}_{k_e}$, ${\left({\displaystyle \frac{\partial {\phi }_j}{\partial y}},{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\right)}_{k_e}$ and ${\left({\displaystyle \frac{\partial {\phi }_j}{\partial y}},{\displaystyle \frac{\partial {\phi }_i}{\partial y}}\right)}_{k_e}$, we have

$$\begin{array}{c}{\left(\frac{\partial {\phi }_j}{\partial x},\frac{\partial {\phi }_i}{\partial x}\right)}_{k_e}=|detJ{T}_e|T(1,:)·Sm·T{(1,:)}^{\prime };\\ {\left(\frac{\partial {\phi }_j}{\partial x},\frac{\partial {\phi }_i}{\partial y}\right)}_{k_e}=|detJ{T}_e|T(1,:)·Sm·T{(2,:)}^{\prime };\\ {\left(\frac{\partial {\phi }_j}{\partial y},\frac{\partial {\phi }_i}{\partial x}\right)}_{k_e}=|detJ{T}_e|T(2,:)·Sm·T{(1,:)}^{\prime };\\ {\left(\frac{\partial {\phi }_j}{\partial y},\frac{\partial {\phi }_i}{\partial y}\right)}_{k_e}=|detJ{T}_e|T(2,:)·Sm·T{(2,:)}^{\prime };\end{array}$$

Using this notation, the term ${(\nabla {\phi }_j\parallel \nabla {\phi }_i)}_{k_e}$ is as follows:

$$\begin{array}{ccc}\hfill {(\nabla {\phi }_j\parallel \nabla {\phi }_i)}_{k_e}& =& {\left({\displaystyle \frac{\partial {\phi }_j}{\partial x}},{\displaystyle \frac{\partial {\phi }_i}{\partial x}}\right)}_{k_e}+{\left({\displaystyle \frac{\partial {\phi }_j}{\partial y}},{\displaystyle \frac{\partial {\phi }_i}{\partial y}}\right)}_{k_e}\hfill \\ & =& |detJ{T}_e|·[T(1,:)·Sm·T{(1,:)}^{\prime }+T(2,:)·Sm·T{(2,:)}^{\prime }].\hfill \end{array}$$

Third case:

$$\begin{array}{cc}\hfill {\left({\displaystyle \frac{\partial {\phi }_j}{\partial x}},{\phi }_i\right)}_{k_e}& ={\int \int }_{k_e}{\displaystyle \frac{\partial {\phi }_j}{\partial x}}{\phi }_{i}dydx\hfill \\ & ={\int \int }_{\widehat{k}}\left(T_{11}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}+T_{12}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}\right){\phi }_{\widehat{i}}|detJ{T}_e|d\eta d\xi \hfill \\ & =|detJ{T}_e|{\int \int }_{\widehat{k}}\left(T_{11}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}{\phi }_{\widehat{i}}+T_{12}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}{\phi }_{\widehat{i}}\right)d\eta d\xi \hfill \\ & =|detJ{T}_e|\left[T_{11}{\int \int }_{\widehat{k}}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}{\phi }_{\widehat{i}}d\eta d\xi +T_{12}{\int \int }_{\widehat{k}}{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}{\phi }_{\widehat{i}}d\eta d\xi \right]\hfill \\ & =|detJ{T}_e|\left[T_{11}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}+T_{12}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\right]\hfill \\ & =|detJ{T}_e|\left(T_{11}\hspace{1em}{T}_{12}\right)\left(\begin{array}{c}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\\ {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\end{array}\right)=|detJ{T}_e|T(1,:)\left(\begin{array}{c}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\\ \\ {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\end{array}\right).\hfill \end{array}$$

Repeating the process to ${\left({\displaystyle \frac{\partial {\phi }_j}{\partial y}},{\phi }_i\right)}_{k_e}$, we have

$${\left({\displaystyle \frac{\partial {\phi }_j}{\partial y}},{\phi }_i\right)}_{k_e}=|detJ{T}_e|T(2,:)\left(\begin{array}{c}{\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}\\ {\left({\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }},{\phi }_{\widehat{i}}\right)}_{\widehat{k}}.\end{array}\right)$$

When dealing with integrals over boundary elements (line segments), the same notation as for triangles is used, noting that they have different dimensions. The $\phi$ in the reference element $\widehat{k}$ are

$${\phi }_{\widehat{1}}\left(\xi \right)=1-\xi ;$$

(27)

$${\phi }_{\widehat{2}}\left(\xi \right)=\xi .$$

(28)

The linear transformation $T_e:\widehat{k}\to k_e$ that relates a reference element $\widehat{k}$ to a generic element $k_e$; that is, that relates

$$\begin{array}{c}\hfill 0\to (x_1,y_1) \mathrm{and}\\ \hfill 1\to (x_2,y_2),\end{array}$$

is defined by

$$T_e\left(\xi \right)=\left(\begin{array}{c}{x}_2-x_1\\ y_2-y_1\end{array}\right)\xi +\left(\begin{array}{c}{x}_1\\ y_1\end{array}\right)$$

which is equivalent to

$$\begin{array}{c}\hfill x\left(\xi \right)=x_1+(x_2-x_1)\xi ;\\ \hfill y\left(\xi \right)=y_1+(y_2-y_1)\xi .\end{array}$$

Thus, we have ${\phi }_{\widehat{i}}\left(\xi \right)=({\phi }_i(x\left(\xi \right),y\left(\xi \right))$ and

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial T_e}{\partial \xi }}={\displaystyle \frac{\partial x}{\partial \xi }}\overrightarrow{i}+{\displaystyle \frac{\partial y}{\partial \xi }}\overrightarrow{j}=(x_2-x_1)\overrightarrow{i}+(y_2-y_1)\overrightarrow{j}.\end{array}$$

Regarding the border, there are two cases.

First case:

$${({\phi }_j,{\phi }_i)}_{k_e}={\int }_{k_e}{\phi }_j{\phi }_{i}dx={\int }_{\widehat{k}}{\phi }_{\widehat{j}}\left(\xi \right){\phi }_{\widehat{i}}\left(\xi \right)∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥d\xi =∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥{\int }_{\widehat{k}}{\phi }_{\widehat{j}}{\phi }_{\widehat{i}}d\xi =∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥{({\phi }_{\widehat{j}},{\phi }_{\widehat{i}})}_{\widehat{k}}$$

Second case:

$${(1,{\phi }_i)}_{k_e}={\int }_{k_e}{\phi }_{i}dx={\int }_{\widehat{k}}{\phi }_{\widehat{i}}\left(\xi \right)∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥d\xi =∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥{\int }_{\widehat{k}}{\phi }_{\widehat{i}}d\xi =∥{\displaystyle \frac{\partial T_e}{\partial \xi }}∥{(1,{\phi }_{\widehat{i}})}_{\widehat{k}}$$

These inner products are used to determine the elements of the M, W, H, and N matrices, as defined in (20). These internal products over the elements and contour elements are combinations of integrals over the reference element (triangle and reference line segment). These integrals over the reference elements are calculated, and their results are stored in matrices nominated for stiffness submatrices.

Replacing ${\phi }_{\widehat{i}}$ and ${\phi }_{\widehat{j}}$ with (23), (24), and (25) in the integrals to be calculated, and considering all combinations between functions, we have

$$\begin{array}{cc}\hfill {\left({({\phi }_{\widehat{i}},{\phi }_{\widehat{j}})}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right)\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}{\xi }^2& \xi \eta & -\xi \left(\xi +\eta -1\right)\\ \xi \eta & {\eta }^2& -\eta \left(\xi +\eta -1\right)\\ -\xi \left(\xi +\eta -1\right)& -\eta \left(\xi +\eta -1\right)& {\left(\xi +\eta -1\right)}^2\end{array}\right)d\eta d\xi \hfill \\ & =\left(\begin{array}{ccc}{\displaystyle \frac{1}{12}}& {\displaystyle \frac{1}{24}}& {\displaystyle \frac{1}{24}}\\ {\displaystyle \frac{1}{24}}& {\displaystyle \frac{1}{12}}& {\displaystyle \frac{1}{24}}\\ {\displaystyle \frac{1}{24}}& {\displaystyle \frac{1}{24}}& {\displaystyle \frac{1}{12}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }{\displaystyle \frac{\partial }{\partial \xi }}\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \xi }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}1& 0& -1\\ 0& 0& 0\\ -1& 0& 1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}{\displaystyle \frac{1}{2}}& 0& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0\\ -{\displaystyle \frac{1}{2}}& 0& {\displaystyle \frac{1}{2}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }{\displaystyle \frac{\partial }{\partial \eta }}\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \xi }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}0& 0& 0\\ 1& 0& -1\\ -1& 0& 1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}0& 0& 0\\ {\displaystyle \frac{1}{2}}& 0& -{\displaystyle \frac{1}{2}}\\ -{\displaystyle \frac{1}{2}}& 0& {\displaystyle \frac{1}{2}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \xi }},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }{\displaystyle \frac{\partial }{\partial \xi }}\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \eta }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}0& 1& -1\\ 0& 0& 0\\ 0& -1& 1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}0& {\displaystyle \frac{1}{2}}& -{\displaystyle \frac{1}{2}}\\ 0& 0& 0\\ 0& -{\displaystyle \frac{1}{2}}& {\displaystyle \frac{1}{2}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\displaystyle \frac{\partial {\phi }_{\widehat{i}}}{\partial \eta }},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }{\displaystyle \frac{\partial }{\partial \eta }}\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \eta }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& -1\\ 0& -1& 1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}0& 0& 0\\ 0& {\displaystyle \frac{1}{2}}& -{\displaystyle \frac{1}{2}}\\ 0& -{\displaystyle \frac{1}{2}}& {\displaystyle \frac{1}{2}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\phi }_{\widehat{i}},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \xi }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \xi }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}\xi & 0& -\xi \\ \eta & 0& -\eta \\ 1-\eta -\xi & 0& \xi +\eta -1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}{\displaystyle \frac{1}{6}}& 0& -{\displaystyle \frac{1}{6}}\\ {\displaystyle \frac{1}{6}}& 0& -{\displaystyle \frac{1}{6}}\\ {\displaystyle \frac{1}{6}}& 0& -{\displaystyle \frac{1}{6}}\end{array}\right);\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\left({\left({\phi }_{\widehat{i}},{\displaystyle \frac{\partial {\phi }_{\widehat{j}}}{\partial \eta }}\right)}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{c}\xi \\ \eta \\ 1-\xi -\eta \end{array}\right){\displaystyle \frac{\partial }{\partial \eta }}\left(\xi \hspace{1em}\eta \hspace{1em}1-\xi -\eta \right)d\eta d\xi \hfill \\ & ={\int }_0^1{\int }_0^{1-\xi }\left(\begin{array}{ccc}0& \xi & -\xi \\ 0& \eta & -\eta \\ 0& 1-\eta -\xi & \xi +\eta -1\end{array}\right)d\eta d\xi =\left(\begin{array}{ccc}0& {\displaystyle \frac{1}{6}}& -{\displaystyle \frac{1}{6}}\\ 0& {\displaystyle \frac{1}{6}}& -{\displaystyle \frac{1}{6}}\\ 0& {\displaystyle \frac{1}{6}}& -{\displaystyle \frac{1}{6}}\end{array}\right).\hfill \end{array}$$

For the case of boundary elements, replacing ${\phi }_{\widehat{i}}$ and ${\phi }_{\widehat{j}}$ by (27) and (28) in the integrals to be calculated, and considering all combinations between functions, we have

$$\begin{array}{cc}\hfill {\left({({\phi }_{\widehat{i}},{\phi }_{\widehat{j}})}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}& ={\int }_0^1\left(\begin{array}{c}1-\xi \\ \xi \end{array}\right)\left(1-\xi \hspace{1em}\xi \right)d\xi ={\int }_0^1\left(\begin{array}{cc}{\left(\xi -1\right)}^2& -\xi \left(\xi -1\right)\\ -\xi \left(\xi -1\right)& {\xi }^2\end{array}\right)d\xi \hfill \\ & =\left(\begin{array}{cc}{\displaystyle \frac{1}{3}}& {\displaystyle \frac{1}{6}}\\ {\displaystyle \frac{1}{6}}& {\displaystyle \frac{1}{3}}\end{array}\right);\hfill \end{array}$$

i.e., similarly, for the second inner product, we have

$${\left({(1,{\phi }_{\widehat{j}})}_{\widehat{k}}\right)}_{\widehat{i}\widehat{j}}={\int }_0^{1}1\left(\begin{array}{c}1-\xi \\ \xi \end{array}\right)d\xi =\left(\begin{array}{c}{\displaystyle \frac{1}{2}}\\ {\displaystyle \frac{1}{2}}\end{array}\right)$$

Thus, when assembling the matrices M, W, H, and N, which make up the system (21), each element corresponds to combinations of internal products between the base functions and/or their derivatives, as defined in (20). These matrices are symmetric, and the number of nodes in the finite element mesh determines their sizes.

The research used the free software Gmsh to construct the mesh. Based on data made available by the Brazilian Institute of Geography and Statistics (IBGE), the domain chosen for the diffusion–advection equation was the territory limited by the state of Mato Grosso. The mesh has 5,417,224 triangles and 2,712,315 nodes, of which 7404 are border nodes and 2,704,911 are internal nodes.

3. Results

This research aimed to model the spatial and temporal spread of the fungus that causes Asian soybean rust. For this purpose, a diffusion–advection equation was constructed whose coefficients were obtained from estimates based on actual data. To this end, the disease occurrence records made available by the Antirust Consortium for the 2015/2016 to 2019/2020 harvests and meteorological data provided by the National Institute of Meteorology—INMET for the period between 2018 and 2019 were used as a reference.

Spores that cause airborne diseases can be treated as particulate pollutants, and their dispersal can also be modeled using an Eulerian diffusion–advection approach [4]. Therefore, we chose to use this equation to spatially and temporally describe the spread of the fungus PP. Furthermore, numerical simulations allowed the creation of more than a hundred maps, which reproduced how the pathogen spread over the studied domain throughout the large part of a crop year. This section presents the developments and results for the terms constituting this equation and maps representing the entire period.

3.1. Obtaining the Terms of the Equation

Equation (2) presents the four terms adjusted from one equation to another. Thus, the definition of the diffusion coefficient, the velocity field (transport), the mortality rate (decay), and the source term guarantees each diffusion–advection equation is a unique model for the most diverse phenomena. Some parameters in numerical simulations varied in time and space, as with transport and decay. The diffusion coefficient was based on research by Diniz (2003) [12], and the values signed to the source on research by Melching et al. (1979) and Alves et al. (2006) [23,24]. In the following sections, the constructions and determinations of the equation coefficients are explained in detail.

3.1.1. Diffusion

The diffusion coefficient can be understood as the pollutant’s dispersion capacity. For the numerical simulations, the parameter definition was based on Diniz’s research [12], which simulated the propagation of pollutants in an air–water system in six different scenarios.

Table 1. Diffusion Coefficient in Air [12].

Intensity	Value (km2/h)
Moderate	0.11
	0.15
	0.30
High	0.50

summarizes all values considered in the simulations.

The four values were tested in numerical simulations. However, the relationship between the advective component and the diffusion coefficient can cause severe instability problems for numerical approximation methods [12]. To guarantee this stability, the diffusive term adopted for this equation was based on the parameter considered to be of high intensity. It was kept constant throughout the domain $\alpha =0.60$ ${\mathrm{km}}^2$/h, equivalent to approximately 100 ${\mathrm{km}}^2$ per week.

3.1.2. Transport

The vector field that carries the pathogen from one place to another in the system is represented in the equation by transport. Its construction was based on wind direction and speed data, which are available at INMET automatic stations. Thus, a corresponding vector was obtained at all points of the grid throughout the studied period, enabling the spatiotemporal variation of the advective term. Figure 2 presents a section of the vector field at three different scales.

Additionally, vectors parallel to highways were added to the vector field at all mesh nodes within 1 km of a major road. The dimensions of these vectors are proportional to the distance from the highway (the closer they are, the greater the intensity of the vector). In this way, it was possible to simulate the influence of vehicle circulation on the vector field.

3.1.3. Decay

Decay represents the pathogen’s mortality rate. To estimate this term, a fuzzy system was built, using as input variables the favorability data (conditions favorable to the disease’s progress), the host map (locations where soybeans are planted), and the application of agricultural pesticides (Fungicide). Figure 3 presents the structure of the fuzzy controller.

For infection and sporulation to occur, environmental conditions must be favorable. Therefore, the temperature associated with the leaf wetness duration is crucial for maintaining spores [25,26]. This data determines the favorability of the fuzzy system created by Alvet et al. (2011) [27]. This variable constitutes one of the fuzzy sets that determines the decay. Figure 4 presents MB, which refers to values for favorable conditions: Very Low, Low (B), Medium (M), High (A), and Very High (MA).

Furthermore, the survival conditions of fungus spores improve as the distance to properties decreases. Thus, the entry Host refers to the distance the point is from a soybean plant. Points up to 5 km away are considered to be in the vicinity of producing properties. They will be classified between Close (P) and Far (D), as seen in Figure 5.

The variable Fungicide is related to the application (A) or non-application (NA) of chemical control in plantations or areas susceptible to soybean plants (Figure 6). In areas where the climate is highly favorable or depending on disease pressure in the region, preventive applications of fungicides guarantee disease control [28]. Furthermore, choosing fungicides and applying the technology are essential steps [29].

These three indicators influence the variable Sigma, representing the diffusion–advection equation’s decay coefficient ($\sigma$) (Figure 7). It assumes values between 0 and 0.4 spores/hour. The rule base that determines the relationship between these variables comprises six rules, as seen in

Table 2. Fuzzy Decay Rule Base.

1. IF favorability is MB AND hosts is P AND fungicide is NA THEN sigma is M. 2. IF favorability is B AND hosts is P AND fungicide is NA THEN sigma is M. 3. IF favorability is M AND hosts is P AND fungicide is NA THEN sigma is B. 4. IF favorability is A AND hosts is P AND fungicide is NA THEN sigma is B. 5. IF favorability is MA AND hosts is P AND fungicide is NA THEN sigma is B. 6. IF host is D OR fungicide is A THEN sigma is A.

.

The Mamdani [30] inference method was used to solve the system, with center of area defuzzification. As each day and node have an associated favorability value, the decay coefficient also has one. By keeping the components of the equation varying in time and space, we seek to bring the simulations as close as possible to reality.

3.1.4. Source

The term source can be interpreted as the places where the pathogen is “produced”. Soybean plants are the main hosts of the fungus. Thus, the legume must exist for there to be source.

Although every soybean plant is eliminated at the end of each harvest, the grains that fall on the roads during product transport germinate during the first rains and can be infected by the fungus, ensuring the existence of a green bridge between one harvest and another [3]. These plants who voluntarily sprout do not receive any treatment to prevent the pathogen’s spread. Therefore, sources were distributed on all the main highways in the state.

Figure 8 shows where Asian soybean rust was recorded in the 2015/2016, 2016/2017, 2017/2018, 2018/2019, and 2019/2020 harvests (represented in the figure by *). Most of them are located close to these highways, reinforcing the hypothesis of considering points close to highways as inoculum sources.

The edges of the highways are public responsibility and correspond to 20 m to the right and 20 m to the left of the setback lane [31]. Therefore, it was decided to increase this limit so that every node in the network located up to 50 m away from a main highway was considered a source for constructing the equation.

To determine the value to be adopted for this parameter (quantity of spores produced at these points), it was taken into account that, due to the injury to the soybean leaves, at least three spores (upper and lower) are produced after two weeks of vaccination. As there can be up to 50 lesions per ${\mathrm{cm}}^2$ of lesion, each ${\mathrm{cm}}^2$ produces approximately 75 spores per week [23,24]. Comparing these data with the average area of the triangles that make up the mesh, we arrive at $f=135\times {10}^6$ spores/(${\mathrm{km}}^2$/week) for each node considered a source.

3.2. Maps

Numerical simulations allowed the production of 182 maps that represented the evolution of the fungus PP throughout the 2018/2019 harvest in the state of Mato Grosso (Brazil). 31 October, the date of the first occurrence of the disease recorded in the harvest, was considered the zero instant. Concentrations were calculated until 30 April, when the soybean cycle ends in the state.

The machine used to run all tests and scenarios was a server containing 144 GB of RAM and two Intel Xeon Hexa-core processors, totaling six cores each. The algorithms created were improved to be parallelized. There were 1400 iterations, which took an average of 170 s to generate the stiffness matrices and calculate the solution in each time interval.

The constants, which, together with the variable parameters, constituted the terms (20) of these matrices, are summarized in

Table 3. Parameters used in the diffusion–advection equation.

Parameter	Value	Unit of Measurement
$\Delta t$	0.02	week
$\alpha$	100.8	${\mathrm{km}}^2$/week
$\theta$	$75\times {10}^5$	spore/(km/week)
f	$135\times {10}^6$	spore/(${\mathrm{km}}^2$/week)

. Knowing that lesions on infected plants produce approximately 75 spores/(${\mathrm{cm}}^2$/week), carrying out the necessary conversions, $\theta$ (boundary condition) was adopted as a constant.

In Figure 9, the dispersion of spores for the last day of each month is represented from 31 October to 31 March, and in Figure 10, spore concentrations (spore/km²) are presented for the points where the disease occurred in the 2018/2019 harvest. In the graphs of this figure, the abscissa axis shows the harvest days (from day 1 to day 182), and the ordinate shows the spore concentrations for each location (from 0 to 1.5 million spores/km²). The highest concentration recorded for these locations was 1,340,409 spores/km² in Lucas do Rio Verde, whose first recorded occurrence was 18 January.

The first signs of the disease take between 5 and 7 days to appear [2]. Expanding this interval a little, the concentration data for these locations were recorded for the interval between the ninth and fifth day before the occurrence was recorded [32].

The first signs of the disease take 5 and 7 days [2]. Expanding this interval a little, the concentration data for these locations were recorded between the ninth and fifth day before the occurrence was recorded [32].

Table 4. Concentration of spores in locations where occurrences were recorded.

City	Record	Average	Max
Sapezal	13 December 2018	118,538	291,443
Jaciara	18 December 2018	261,927	345,987
Tangará da Serra	26 December 2018	5664	13,479
Campo Novo do Parecis	27 December 2018	93,246	230,093
Querência	28 December 2018	645,929	844,275
Campo Verde	28 December 2018	605,664	1,023,302
Comodoro	3 January 2019	292,909	386,616
Gaúcha do Norte	8 January 2019	251,319	311,625
Ipiranga do Norte	10 January 2019	5606	70,725
Campos de Júlio	14 January 2019	198,609	283,814
Sorriso	15 January 2019	668,136	899,254
Lucas do Rio Verde	18 January 2019	1,253,994	1,340,409
Primavera do Leste	21 January 2019	287,791	534,270
Nova Mutum	22 January 2019	504,006	540,787
Feliz Natal	22 January 2019	1888	13,077
Itiquira	28 January 2019	1883	12,397
Rosário Oeste	22 February 2019	295,360	352,821
Sinop	15 March 2019	412,856	866,666

contains the dates on which the occurrence of the disease was recorded, the averages of the five selected data (from the fifth day before registration to the ninth), and the maximum value that the concentration reached in that location during the period studied. In Gaúcha do Norte, the highest concentration value recorded occurred within the probable infection range.

4. Discussion and Conclusions

For Asian soybean rust to develop, an interaction between three factors is necessary, the so-called tripod of the disease: host, pathogen, and favorable climatic conditions [33]. This research was used to model the spatiotemporal dynamics of the fungus that causes this disease.

Throughout the entire period studied, data from 182 days were obtained. These data were used as input variables for a fuzzy system that more accurately modeled the temporal and spatial development of Asian soybean rust for the state of Mato Grosso [32]. In this research, the controller was correct in his prediction in 81% of the records, and according to another analysis, the accuracy percentage was 92%.

Obtaining actual values experimentally to model the terms that make up the differential equation is often unfeasible. The construction of the terms of the equation used in this research considered several nuances that influence the phenomenon based on fuzzy logic. This efficient tool for dealing with subjectivities allowed actual data that are already commonly captured, such as temperature and relative humidity, to bring these coefficients closer to reality, which can be considered an advance.

The fuzzy system created for the decay, the vector field that reproduces the interference of traffic in the pathogen dispersion, and the implantation of sources along the highways are examples of how the different variables that act in the studied system were considered. The fuzzy logic that estimates the effect of meteorological variables on mortality rate could be improved through deeper investigation into the impact of agricultural pesticides on decomposition. Installing sources along highways satisfactorily modeled the evident problem of soybean plants along highways. The simulations showed that along the main roads, there is a higher concentration of spores.

During the development of this research, it was possible to realize that although there is a lot of research to understand Asian soybean rust, few of them intend to understand the macro behavior of the fungus. Although it is essential to study it, there are limitations regarding implementing the equation for a domain of this magnitude. Even though it is not simple, replicating this research in other soybean-producing states and neighboring countries that also produce the grain would allow a more comprehensive perception of how the fungus behaves.