A Framework Based on Finite Element Method (FEM) for Modelling and Assessing the Affection of the Local Thermal Weather Factors on the Performance of Anaerobic Lagoons for the Natural Treatment of Swine Wastewater

Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Department of Process Engineering, Institute for Environmental Studies and Natural Resources (i-UNAT), University of Las Palmas de Gran Canaria (ULPGC), 35017 Las Palmas, Spain; alejandro.ramos@ulpgc.es (A.R.-M.); sebastianovidio.perez@ulpgc.es (S.O.P.-B.); carlos.mendieta@ulpgc.es (C.M.-P.); federico.leon@ulpgc.es (F.L.-Z.) * Correspondence: saulo.brito101@alu.ulpgc.es † These authors contributed equally to this work.


Introduction
Anaerobic digestion (AD) is an eco-friendly biological process which is universally used for the treatment of agricultural, industrial and municipal wastewater around the world [1][2][3][4]. Its utilization is increasingly widely, due to its capacity for producing methane, which can be used afterwards as a heat source or for electricity generation, taking part within the low-carbon energy technologies and circular bio-economy [5]. In this context, anaerobic lagoons (AL) are natural wastewater treatment systems with a long hydraulic retention time, suitable for small communities due to the low energy demand and the operating costs [6][7][8][9]. By applying this kind of technology, the mechanical equipment, used for mixing processes in conventional plants, are avoidable. In addition, AL offer a number of advantages, such as the establishment of concentration profiles along the reactor, a buffering capacity in cases of overloads and greater protection against acidification [10]. However, due to the fact that AD is strongly influenced by temperature, there is a close dependence between AL and weather conditions, so its implementation may be limited in cold or low solar radiation areas [8,11].
The application of mathematical models builds understanding for both microbialrelated dynamic and kinetic processes, reveals optimisation possibilities, which lastly improves the digester's performance [12,13]. The IWA Anaerobic Digestion Model No.1 (ADM1) [14], created in 2002 to stablish a common platform for the modelling of AD processes [15], has been widely applied in waste treatment processes, due to its high feasibility, considering the fact that most of the processes of AD are included within ADM1 [16]. However this model has merely been applied to completely mixed reactors. The approach of models based on the ADM1 for unstirred waste water treatment systems has been little studied. In these models, complexity is increased and the effect of boundary conditions is essential. Moreover, the mathematical complexity required by these models does not entail a significant issue, due to the increasing technological and computational development [12].
In the past twenty years many researches based on mathematical models for treatment processes in lagoons have been carried out. Fleming [17] created the first models applying computational fluid dynamic (CFD) for the prediction of the performance of full-scale incompletely mixed anaerobic digesters. Wu and Chen [8] developed a CFD model for AL which combines physical and biological processes, and includes both heat conduction and solar radiation by a thermal model. In this model, a single-phase incompressible Newtonian fluid is considered. Goodarzi, Sookhak Lari, and Mossaiby [18] determined the effect of ambient and inlet temperature variations on the hydraulic performance of a typical rectangular pond. In all these described models, the biological processes are depicted by a single equation depending on the concentrations of the influent and effluent. Brito-Espino, Ramos-Martín, Pérez-Báez, and Mendieta-Pino [19] defined advection, diffusion and reaction phenomena for wastewater treatment in anaerobic plug flow reactors by non-linear, second order, partial differential equations. ADM1 is implemented within this model, and both biochemical and physical-chemical reactions of ADM1 are calculated by a flowchart for sequential processes. In this method, temperature is not considered. Nevertheless, very few researches have been conducted to develop a comprehensive model which integrates fluid flow, heat transfer, and cells behaviour in AL.
The aim of this work is to set-up a theoretical framework for wastewater treatment in unstirred flow anaerobic lagoons, by a model which allows the integration of fluid flow, heat transfer and cells behaviour, for the purpose of describing processes occurring in AL. The implementation of the ADM1 into the model and the consideration of the influence of the local thermal weather, identified with the boundary conditions, allows the model to portray the processes taking place in reality more precisely than [19]. In order to do this, an improved two dimensional mathematical model, based on the coupling of a set of parabolic partial differential equations (PDEs) and related to the phenomena associated to AL, has been developed. In addition, Dirichlet, Neumann and Robin boundary conditions have been established on the differential equations. This model combines the parametrization of different processes within the lagoon and its environment with the finite element analysis. Finally, the parallelization of the resulting algorithm has been performed in the simulation, therefore allowing an improved computational efficiency than the resulting form sequential processes in [19]. Thanks to the help of FreeFemm++ and the parallel solver package, available for this software, the processing of each one of the variables related to AD processes and the simultaneous exchange of the data has been feasible. Having said this, we conclude that the novelty of this study resides in the following aspects. Firstly, in the implementation of the ADM1 and the heat transfer phenomenon in a mathematical model which describes a unstirred fluid flow, in order to predict the spatial distribution of the different variables that take part in the processes within the AL. Furthermore, secondly, in the optimisation and designing of the algorithm, by parallel method, providing an accurate forecast of the real behaviour of the process, as is shown in the ADM1. In the simulation, two different scenarios have been chosen as examples; the first corresponds to a conventional AL which is subjected to the ambient temperature, and the second includes heat sources, induced by solar assisted [20] or through the biofuel recovery in the anaerobic process [21,22].

Overview
Pollutant is removed in AL through combination of physical, biochemical and physicalchemical phenomena. Advection, diffusion and heat transfer are the most common physical processes in these systems ( Figure 1). Both the organic matter and the suspended microorganism within lagoons are subjected to the mechanical transport with the bulk flow of the water (advection). At the same time, they tend to spread out and diffuse from higher to lower concentration as time varies (diffusion). The energy transfer in the system, due to a temperature gradient (heat transfer), is performed by conduction and convection processes. Atmospheric factors associated with the borders of the model on the Earth's surface include, beside the two previous, radiation process. Digestion process is carried by anaerobic microorganism's activity, bacteria an archaea, through a number of sequential and parallel reactions. The biochemical reactions consist of irreversible five-stage processes; disintegration, hydrolysis, acidogenesis, acetogenesis and methanogenesis reactions. Physical-chemical reactions are those reversible processes where cells are not involved. They are, firstly, the ion association/dissociation, and gas-liquid transfer [14].  Considering a system where the lagoon and its environment are included, different zones can be identified ( Figure 1). Depending on the local parameters-thermal conductivity, specific heat, density, and where biological processes take place-they are considered different zones. Boundary conditions are located on the borders. Undisturbed ground temperature T UGT is a ground thermal property situated at a depth where the ground temperature is approximately invariable, depth value depends on climatic conditions and is different in various regions of the Earth [23,24].

Governing Equations (Strong Formulations)
In this research, the mathematical model proposed is based on the two-dimensional advection-diffusion-reaction, Stokes, and heat transfer equations. This is accompanied by a series of boundary conditions. On the other hand, the IWA Anaerobic Digestion Model 1 (ADM1) has been implemented in the model.
The description of the model has been expressed in terms of primitive variables, mass, velocity, pressure, and temperature. In these equations it has been assumed that velocity and temperature field are in steady state conditions.
Governing equations and boundary conditions are summarized below.
where (φ) is a scalar field that represent concentrations of both substrates and cells of each of the biochemical reactions included in the anaerobic processes, u = (u 1 , u 2 ) is given by Equations (2), f (φ) is the source function, which is positive f (φ) > 0 for growth and production or negative f (φ) < 0 for decay and consumption, biomass and metabolites, respectively, Γ D and Γ N are Dirichlet and Neumann boundary conditions, respectively. This term is developed in Equation (7) , F(x, y) is a generation function, which is zero (0) in this case.

Stokes Equation
Stokes equation, together with the ADR has been used to describe the flow. It is usually used for fluid with slowly motion and with high viscosity [28,29]. In this research, a constant density and incompressible Newtonian fluid flow has been considered. Strong formulation and Dirichlet Γ D and Γ N boundary conditions are as follow.

The Energy Equation-Temperature Distribution
The energy equation is based on the conservation of energy and the Fourier heat conduction laws [30]. The internal energy balance equations, under a steady-state Eulerian description can be expressed as a function of temperature [30,31] where Γ D corresponds to Dirichlet condition. It is applied to the UGT ( Figure 1); Γ N is the Neumann condition. It describes the value of the gradient of the dependent field variable, normal to the boundary. Its calculation is based on Fourier Law; Γ R is the Robin condition. It describes the Earth's surface heating and cooling Γ R . This implied the use of Stefan-Boltzmann's law and Newtons' law of cooling to model the heat exchange, related to radiation and convection processes, respectively [32,33]. Stefan-Boltzmann constant is σ = 5.68 · 10 −8 (W · K −4 · m −2 ) [34]; T a,ext temperature of the externally surrounding surface; ε s is the Earths' surface emissivity, where 0 ≤ ε s ≤ 1 ; T sky is the sky radioactive temperature. These are used to estimate the radiative heat exchange with the Earth's atmosphere [35]. T sky and ε sky is used to estimate the radiative heat exchange with the Earth's atmosphere [35].
Here N are tenths cloud cover, and T dp (K) is the dew-point temperature to which it must be cooled to become saturated. It is obtained by a correlation found in [33] (6) and In this work, heat sources from biochemical reactions have not been considered.

Kinetic Equations
ADM1 is used for the description f (φ) (Equation (1)). Biochemical rate coefficients and kinetic rate equations are represented in the Tables S1 and S2 within the Supplementary Materials Section. First order kinetic was considered for the hydrolysis, acidogenesis, acetogenesis and methanogenesis. The following equations based on common kinetic expressions describe anaerobic treatment processes: f (S i ) and f (X i ) are the changes in substrates and cells concentration over time. These equations are based on the monod-type reaction kinetics [13,36]. In this model, it has been considered free ammonia and pH inhibitions, in addition to the butyrate and valerate competition [19].
The influence of temperature has been obtained by the Cardinal Temperature Model 1 (Appendix A) proposed by [37]

Solution Procedure
The finite element method, numerical technique based on the generation of a finite element geometric model, is used for the solution of the partial differential equations including in the problem.
In this methodology, the major steps include 1. The approach of the weak forms from the governing equations. The solutions are assumed to belong to Hilbert space, considering this space as an infinite dimensional function space with functions of specific properties that can be suitably managed in the same way as ordinary vectors in a vector space. They are represented in Table 1. 2. Discretization of the domains, both physical with more or less regular triangulation and related to time. In Figure 2, the discretization of the different sub-domains, nodes and triangle, is showed. 3. Selection of the shape functions, essential to provide an approximation of the solution within an element. These relate the coordinates of every point of a finite element with the positions of its nodes, 4. Formulation of the system of equations. 5. Solving systems of equations. The free software FreeFem++ has been used to solve them. It is a PDE solver with its own high-level programming language and accurate syntax for mathematical formulation. Freefem++ have high diversity of triangular finite elements (linear and quadratic, Lagrangian elements, discontinuous P 2 , etc.) to solve PDE in two (2D) and three (3D) dimensions.

Model Weak Equations Hilbert Spaces
Stokes

Calculation
The partial differential equation solver FreeFem++ was used to implement the algorithm for the calculation. Due to it advantages, open access software with a powerful generated mesh and a large collection package to visualize approximate solutions, makes Freefem++ an ideal tool to solve complex partial differential equations [38]. Parallel calculation by parallel computing on clusters of personal computer has been achievable with a Message Passing Interface (MPI) within Freefem++.

Model's Considerations
In order to solve the formulated problem, fitted for a specific case, it has been necessary to stablish geometric conditions, physical properties, initial conditions and boundary values. A summary of these characteristics is represented in Figure 2, Tables 2-4 and in the Supplementary Materials Section.
The system has been divided in different domains and subdomains. Ω 1 , Ω 2 and Ω 3 are included within Ω g and refer to the immediate ground around of lagoon; whereas Ω 4 , Ω 5 and Ω 6 in Ω r and concern the lagoon. AD occurs in Ω 5 considering Ω 4 and Ω 6 as transition zones (Figure 2).
In this case, the proposed anaerobic lagoon is located in temperate zones and is subjected to the environmental thermal conditions considering that there are no thermal loads on the sides of the domains, so the ground heat flow is transmitted vertically. Table 2. General parameters considered. Q represents hydraulic flow in the inlet and outlet pipe. S i and X i represent substrate and cell concentrations in the inlet pipe.   Figure 3 shows the temperature distribution in the proposed system under steadystate conditions for some examples of wastewater treatment plants whose information is included in Table 4. The lagoon contour has been illustrated in the first graphic. It is of interest to observe how the thermal behaviour within the lagoon depends on the boundary conditions, but also the hydraulic flow that is subject to the boundary values in the inlet and outlet pipe.  Table 4. Table 4. Specific weather parameters considered from four wastewater treatment plants (WWTP). It is provided the Universal Transverse Mercator (UTM) coordinates. T am are the annual means temperatures and T mm the monthly means. P1 and P3 are located in the coastal zones, while P2 and P4 are located in the mid-altitude zones.  Figure 4 represents the variation on concentrations, in steady state, happening in some of the processes taking place within the lagoon. These simulations describe, in the case of P1 (Table 4), the variation of sugar, propionate, acetic acid and their corresponding bacterial biomass concentration, along of the pond's length and depth, according to the boundary conditions as shown in Table 2.  As shown in these simulations, concentrations decreased throughout the pond length as a result of dispersion and biodegradation. The transition zones (Figure 2) have been considered as low microbial activity, so the net growth of cells is observable from the x-axis value equal to 2. As a result, propionate and acetic acid's source are located in this zone (Ω 5 ). The cells' growth is affected much more by the concentration of the substrates than the temperature's effects because temperature variations in this region differ very little from one point to another.

Organic Matter Removal and Behaviour of the Microbial Community
With respect to acidogenesis and methanogenesis, cells effectiveness on substrates removal is greater than in the acetogenesis due to the kinetic parameters. For the propionate µ max and K s are 0.49 (d −1 ) and 1.145 (mg(COD) · m −3 ), respectively, (see Table 3). Thus, cells growth value and, therefore, the substrate removal is lower than in the previous two cases (see Equation (7)). The resulting value of substrate concentration in the outlet, after the wastewater treatment process, is between 600 and 500 (mg(COD)/L). In the case of acetic acid, this same concentration, next to the pond outlet, is smaller, due to the accumulation of organic matter that has not been reached by the microbial community. Figure 5 shows charts representing substrate and propionic acid bacteria's concentrations for P1 (Table 4) along the axis AA. Cases 1 and 2, with different concentrations of microorganisms in the inlet pipe of the lagoon, 110 (mg(COD)·L −1 ) and 115 (mg(COD)·L −1 ), respectively, are compared. In both cases, the graphics share a similar trend, a downward slope which, al last, connects at the middle point of the axis. There is no net growth within the microbial population. The slope above mentioned is reduced, cells growth offsets the diffusion process in Ω 5 (Figure 2). Nevertheless, substrate removal in Case 2 is benefited by the highest concentration of cells at the inlet pipe 150 (mg(COD)·L −1 ). The decrease in values from the last section, for both graphs, occurs as a result of the dispersion phenomenon since microbial activity in Ω 6 has not been considered.  The lagoon is subjected to ambient temperature (see case P4 in Table 4 and Figure 3). E-2; It is included a bed heat source at the bottom of the lagoon, between 2 and 3 coordinates of the x axis with a temperature of 35 • C. E-3; In this occasion, that same heat source is located between coordinates 4 and 8 of the x axis. The distribution of temperature is showed in Figures 3 and 6. As expected, the removal efficiency is improved by the rising temperature of the heat source, as is observed in the cases 2 and 3. However, this graphic also shows that organic matter is eliminated more efficiently in case 2 than in 3 in a percentage of 10 %. Consequently a minor residual concentration in the outlet pipe is achieved.  By analysing the table, it can be said that the propionate concentration at the outlet of the lagoon, once acetogenic bacteria have removed a great part of substrate in the system, is among 280 and 135 mg·L −1 . By placing a heat source, strategically, at the bottom of the lagoon (E-2 and E-3) it is possible to reduce substrate concentrations at its outlet.

Conclusions
In this paper, we have proposed and assessed a methodology for anaerobic cells performance for wastewater treatment, in AL, under the influence of the temperature. It has been studied in terms of biomass and substrate concentrations. The model couples a series of PDEs, related to the phenomena associated to AL (ADRE, ADM1, Stokes and heat transfer), to each other.
Diffusion for horizontal and vertical directions, the movement of the bulk of the concentrations in accordance with a gradient, external temperature interactions, biochemical and physical-chemical reactions, and a set of boundary values were considered in this study.
This model builds understanding for microbial community's behaviour along the lagoon as a function of the temperature. Applying heat load in different points of the system, it has been possible to establish correlations through the graphics, as well as the comparison between diverse scenarios according to their corresponding boundary values. The results give us the possibility to obtain effective designs adapted to each circumstance, avoiding energy loss.
This methodology allows the optimization of unstirred flow systems, taking into account that the advantages of these systems make them more suitable for specific applications. The model can be used in the prediction of the effluent quality and in the design of AL to achieve better performances.
In view of the results, it can be concluded that this methodology has significant potential as a tool for both the design of AL, and the interactive learning of the microbial ecology in plug flow systems.

Abbreviations
The following abbreviations are used in this manuscript: