548 American Academic Scientific Research Journal for Engineering, Technology, and Sciences ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 http://asrjetsjournal.org/ Finite Element Modal Analysis of Transient Water Flow in Aquifers Ricardo Pereira Guedes a *, Leonardo José Do Nascimento b , Paulo Marcelo Vieira Ribeiro c , Cicero Vitor Chaves Junior d , Jean Baptiste Joseph e a,b,c,d,e Dept. of Civil Engineering, Federal University of Pernambuco, Av. Prof. Moraes Rego,1235, Cidade Universitária, Recife, 50670-901, Pernambuco, Brazil a Email: ricardo.pguedes@ufpe.br, b Email leonardo.guimaraes@ufpe.br c Email: paulo.vribeiro@ufpe.br, d Email: cicero.vitor@ufpe.br e Email: jean.baptiste@ufpe.br Abstract In this paper a modal superposition method is applied for the numerical modeling of aquifers. The proximity of aquifers to populated regions requires special care in their management to avoid problems that affect the quantity and quality of the water they supply. To contribute to the management of this type of natural resource, we propose a numerical strategy based on modal analysis using the finite element method. This procedure assist water production scenarios, performing the mass balance where water extraction is done through wells, in aquifers that are subject to natural recharge. This mathematical procedure is based on the modal superposition for transient flow in porous media. To evaluate its efficiency, this strategy was compared with the classical finite element method. The advantage of the proposed method resides in the possibility of reusing the properties of the global matrix of the finite element method in transient problems, for different production conditions given by the distributed recharge and by the water extraction rate from the wells, solving the numerical problem with a more efficient use of computational resources. This strategy is useful in studies of uncertainty quantification, history matching and optimization of water production in aquifers, since these types of analysis are resource intensive for the very large number of numerical simulations required for these scenarios. Keywords: modal superposition method; finite element method; numerical simulations; aquifer. 1. Introduction An aquifer is defined by [1] as a geological formation that contains water and allows significant amounts of water to move through its porous structure, called porous groundwater reservoir. According to [2], groundwater is essential to maintain life on the planet: it is part of ecosystems and maintain the discharge of rivers, lakes, mangroves and swamps. ------------------------------------------------------------------------ * Corresponding author. http://asrjetsjournal.org/ American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 549 Even though this resource has been used for a long time, groundwater has not yet been used as it should, being often relegated to the background. In recent decades, many questions have been raised about the importance of groundwater from an environmental, social and economic point of view [2]. Many developing countries still do not have an adequate treatment of their aquifers. Even in countries with tropical climate, high rainfall and high supply of surface water in economically developed regions, it is highly recommended that water management governmental agencies give proper importance to the source of groundwater. To give an idea of the magnitude of this source of water supply, in a publication by [3], a survey of the world water balance was presented based on data from UNESCO (1978), and it was quantified that the water on planet Earth is 96.5% distributed in the oceans and 3.5% on the Earth's surface. These 3.5% of the water on the continents are divided into approximately 1% composed of saline groundwater or in saline lakes, leaving only approximately 2.5% of fresh water in the world. Of this fraction of freshwater in the world, 68.6% are in the arctic, such as icebergs and the polar crust, and 30.1% in aquifers and 1.3% on the surface of continents. In other words, aquifer water represents almost 96% of all water available for immediate consumption. This demonstrates the urgent need for an adequate management of this water resource. In order to understand the water cycle on our planet, much has been studied. This is because the water cycle or hydrological cycle comprises all the movement of water through the continent, oceans and the atmosphere, ranging from the evaporation of water to the formation of clouds, its precipitation, generating recharge for rivers, lakes and oceans. Plants and soil also perform an important role in water retention, a portion of which will return to the atmosphere through evaporation. Another portion of the water will infiltrate the soil and feed the groundwater reservoirs in the form of recharge, forming a cycle of vital importance to the supply of all forms of life. In order to put into practice a sustainable management of water resources, numerical tools are needed to assess the different exploitation scenarios of an aquifer, integrating its behavior to the hydrological cycle of a specific region. Natural recharge, for example, is a variable that is difficult to obtain experimentally, which can be obtained through numerical modeling with monitoring of the piezometric head of the aquifer. With information from the scenarios accompanied by real field data, it is possible to reduce the risk of damage, especially those caused in densely inhabited regions. One of the major challenges for the development of large cities is to guarantee the flow of rainwater in urban contexts [4]. Aquifers where volumes of water are extracted at rates greater than natural recharge can be damaged and may present problems such as saline intrusion (salinization of coastal aquifers) or surface subsidence and aquifer compaction due to excessive depressurization, resulting in damage to the region that benefits from the aquifer. Reference [5] pointed out that computational models for groundwater flow play an important role in the management of water resources, allowing decision-making based on flow control in wells and measurable aquifer properties. This article contributes with a computational model to accelerate the simulation of water flow in aquifers. In recent years, many advances have been made with the evolution of processors and mathematical algorithms, American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 550 allowing increasingly faster simulation, enabling studies of optimization and propagation of uncertainties. Furthermore, many problems that people did not dream of solving two decades ago are now routinely solved. [6] The porous media waterflow simulator proposed here has as input data the permeabilities of the geological formations that compose the aquifer as well as the recharge data and water extraction from the wells. In this numerical code in Matlab, the modal superposition method was used, based on discretization in finite elements of the problem domain that allows the analysis process of one aquifer where multiple scenarios are possible for recharge and exploitation of the aquifer be carried out in an agile and accurate way. A case study on a hypothetical aquifer is also presented in this paper in order to demonstrate the effectiveness of the proposed methodology. The method proposed here contributes to the ability to quickly perform several simulations for the same domain with different boundary conditions. According to [7], a real site study requires a series of alternative simulation runs with different boundary and initial conditions. In analyzing an aquifer, each scenario gives a better view of flow processes and shows data uncertainties and model limitations when model outputs are contrasted with field observations and measurements. THEORY 1.1. Flow in porous media Flow through porous media is widely studied in several engineering areas, such as groundwater hydrology, reservoir engineering, soil mechanics, chemical engineering, among others. In general, a porous medium can be defined as a solid matrix with empty spaces [1]. Porous materials can be soil, porous or fractured rocks, ceramics, fibrous aggregates, filter paper, sand filters, among others. To be classified as flow in porous media, there must be a solid matrix with interconnected voids so that the fluid can flow. When the interconnection between the pores is high, there is a high permeability of the geological formation. Aquifers are geological formations with high capacity to transmit and store water, that is, of high hydraulic conductivity (permeability) and porosity [8], respectively. Starting from the macroscopic equation of water mass conservation in a transient flow regime in a three- dimensional porous medium, already considering Darcy's law, we obtain (1) below: ' 0 ; , 1,2,3ij i j h h k q S i j x x t             (1) where kij is the hydraulic conductivity tensor, h is the hydraulic head, S is the specific storage and q' is the source/sink term, which may be due to recharge and water extraction from the wells [9]. In equation (1) and in the other equations in this paper we are using the index notation for the variables, where repeated indexes indicate summation. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 551 Hydraulic conductivity is related to the volumetric flow and the imposed hydraulic gradient, that depends on the porous medium and the viscosity of the percolating fluid, which in this case is the water. The specific storage (S) of an aquifer indicates the relationship between changes in the amount of water stored and the corresponding changes in hydraulic head [1]. 1.2. Numerical solutions to the physical problem To solve the physical problem represented by the equation in partial derivatives (1), based on the solution of a system of ordinary differential equations through modal analysis, the finite element method is initially used. The finite element method is a tool for approximating the solutions of governing equations of various physical phenomena [10]. Reference [11] emphasizes the importance of the finite element method for solving problems from engineering, mathematics, and physics. As well as the finite difference method, the finite element method can also be used to obtain the solution to the problem of water flow in an aquifer. This latter is the method adopted in this paper, without loss of generality, since the same methodology proposed here could also be developed based on the finite difference method. The development of the finite element method for transient flow in porous media is based on the Galerkin procedure, where the following approximate solution is proposed: 1 2 3 1 2 3 ˆ( , , , ) ( , , ) ( ) ; 1,2, ,J Jh x x x t N x x x h t J n  (2) where NJ are the shape functions of the finite element method and hJ the hydraulic head at the nodal points of the finite element mesh. n is the number of nodes in the finite element mesh. In the notation used here, the indices i and j refer to the coordinate axes x1, x2 and x3 and the indices I and J refer to the nodes of the finite element mesh, according to [9]. Following the classical procedure of the Finite Element Method, applying the Weighted Residual Method according to Galerkin's Weak Formulation to the conservation equation (1) and substituting the approximate solution given by equation (2), one arrives at [9]: ´J I J I J ij J I I i jR R B R dh N N SN N dR k h dR N qdB N q dR dt x x            (3) where R is the flow problem domain and B its boundary (regions b1 and b2) where the boundary condition is applied, as shown in Fig. 1. q is the outflow through the boundary, associated with the boundary condition of Neumann, ij i j h q k n x       , with the normal ni pointing out of the domain. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 552 R b1 b2 Figure 1: Flow domain in porous medium discretized by the Finite Element Method. The equation (3) described above can be rewritten in the matrix form [9]: [ ]{ } [ ]{ } { }+ =&S h K h Q (4) In this equation, [S] is the Specific Hydraulic Storage matrix, [K] is the conductance matrix, {Q} which is the vector related to recharge and/or water extraction. The vector {h} and its rate over time { h&} are the unknowns of the problem. The matrices in equation (4) are formed by the composition of the individual elements of the finite element mesh, according to their connectivities,     1 m e e S S   ,     1 m e e K K   ,     1 m e e Q Q   [6], where m is the number of elements in the discretized domain. From equation (3) we can obtain the equations below, for each individual element [9]. The conductance matrix [K] is represented by: e I Je IJ ij i j R N N K k dR x x      , (5) The specific storage array [S] is being represented by, e e IJ I J R S SN N dR  , (6) and finally, the recharge and sink vector {Q} can be represented by, ' e e e I I I B R Q N qdB N q dR    , (7) where this term of the equation will be constant over time From the specific storage matrix e IJS , alternatively, the concentration operation of the specific storage mass can be carried out through diagonalization [12]: ni American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 553 with 0e e e II IJ IJ J S S S for I J   (8) Being e IJS more convenient for solutions in transient problems in situations where it is necessary to invert the global matrices. From the finite element governing equation (4), a code can be generated that will be used in the comparison of the classical finite element method and the one proposed in the next section, with the application of modal superposition. In the classical finite element method, the solution applied for the time step can be approximated using finite differences according to equation (9) [9]:                   1 11 1 1 1 k k k k K S h K S h Q Q t t                         (9) where k and k+1 are the previous and next time steps, respectively, and  is the time weighting factor. According to [9],   1k Q  depends on the values of the flows q and q'. With equation (9) defined, it is possible to have the transient part of the finite element code defined. The solution of the system of non-homogeneous first order linear ordinary differential equations is represented by equation (10): 1 2 ( ) ( ) { } ( )n h t h t h h t ì üï ïï ïï ïï ïï ï= í ý ï ïï ïï ïï ïï ïî þ M (10) 1.3. Modal superposition As an alternative of solution for equation (4), the modal superposition method is introduced, which will be compared with the classical finite element method in equation (9) in terms of CPU time and accuracy. The modal superposition is a method traditionally used in dynamic analysis of structures to obtain the shape and natural frequency modes, the displacement, velocity and acceleration fields of the structural system under study. In particular, it is widely used in vibration analysis of structures [13]. According to [14], in recent decades modal analysis has become one of the main technologies for optimizing the dynamic characteristics of structures [14]. Roy R. Craig, J., & Kurdila, A. J [15] indicate the great use of this method by the advent of computers and the extensive use of Fast Fourier Transform (FFT) that allowed a new paradigm for the structural analysis [15]. He, J., & Fu, Z.-F. [14] comment that the first most significant proposal of modal analysis was in 1947 by C.C. Kennedy and C.D. Pancu (1947), which only became important after J.W. Cooley and J.W. Tukey (1965) developed an FFT algorithm in 1965. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 554 Due to the versatility of the modal superposition method in solving linear systems, it is possible to replicate the technique widely used to solve systems of second-order differential equations, found in systems of dynamic analysis of structures, for systems of first-order differential equations of other phenomena. An example is the equation system for flow in porous media, such as the one described here, where the unknows are the hydraulic head and its rate of change. The advantage of the modal superposition technique is to be able to solve linear systems of differential equations with multiple degrees of freedom and to decouple the equations in simple degrees of freedom, with the use of orthogonalization of the matrices and the change of base, performing a transformation in the vectors and matrices of the linear system. As demonstrated by [16], a base transformation is achieved through the principle of modal superposition applied to a system of multiple degrees of freedom with generalized (modal) coordinates in order to obtain an uncoupled system. In other words, the modal superposition method can transform coupled systems of differential equations into a system of independent or uncoupled equations, where each equation contains only the time variable [16]. The response of the multiple degrees of freedom system by modal superposition is then defined as the sum of the responses of individual modes [13]. 1.4. Mathematical representation The governing equation of the physics of transient flow in porous media is an equation in partial derivatives whose system of equations resulting from its discretization by the finite element method is inhomogeneous. For the solution of this non-homogeneous system of equations, it is necessary to compose two results: a general solution of a non-homogeneous system in any interval I, in addition to a particular solution of the system [17]. In other words, the solution of this inhomogeneous system of equation (4) is composed by = +c ph h h , where ch is the solution for the homogeneous system (general solution) and a ph is the particular solution of the inhomogeneous system, which gives the desired solution of the equation (5). 1.5. Homogeneous solution To calculate the solution of a homogeneous differential system of equations, it is necessary that the recharge/sink term is zero, that is, { } 0Q = , which leads to a solution proposed in equation (11) [17]: 1 2 1 2 n c c c n ch c h c h c h= + +L (11) the solution for each vector i ch can be represented in equation (12): American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 555 1 2 to 1..i i i i t c i n l l h e i n l l ì üï ïï ïï ïï ïï ïï ï= =í ý ï ïï ïï ïï ïï ïï ïî þ M (12) where λ and l are constants in the solution for each i ch . Equation (11) can be written in simplified form [17]: { } [ ] [ ] 1 2 to 1..i i n n i i t t c c i i i n l l h h e L e i n l l l ì üï ïï ïï ïï ïï ïï ï= = = =í ý ï ïï ïï ïï ïï ïï ïî þ å å M (13) Getting a simplification that leads to the formation of matrices and vectors of the system of equation. The matrix L represents the sum of the constants of l. The time derivative is simply presented as, [ ] [ ] [ ]{ } t ch e Lll=& (14) to facilitate the understanding of the subsequent substitutions we can represent equation (13) and equation (14) as, t ch e Ll= and t ch e Lll=& , respectively (15) Substituting now, equation (13) and equation (14), with the simplification represented in equation (15), in equation (4), we have [17]: [ ][ ] [ ][ ] 0t te S L e K Ll ll + = (16) From equation (16), through substitutions, it is possible to obtain an adequate mathematical formulation for the eigenvalues and eigenvectors calculations. To proceed, we need to divide (16) by tel , obtaining [17]:            ( ) 0 or S K L K L S L     (17) With the manipulation of equation (17), it can be compared to the equation below [18]: [ ]{ } { }A x x (18) Which is the classical representation of an eigenvalue and eigenvector problem, where λ is the eigenvalue and x is the eigenvector of the matrix [A]. For the flow formulation in porous media, here represented by equation (17), the matrix [A] in equation (18) is given according to equation (19): American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 556 1[ ] [ ] [ ]A S K  (19) and so, as represented by equation (20):        1 S K L L    (20) where [L] is the modal matrix corresponding to the matrix formed by the matrices whose solution to the equation (20) can be given by, ( ) with 1..i n i t c i i h L c e i n  (21) with n being the total number of equations generated in the system of equations, which coincides with the number of nodes of the finite element mesh, represented as below [17], 1 1 2 2 (1) (2) ( ) 1 1 1 1 (1) (2) ( ) 22 2 2 (1) (2) ( ) modes (eigenvectors) ( ) ( ) ( ) n n n t c tn c c tnc nn n n h t l l l c e h t c el l l h h t c el l l                                      (22) where the equation (22) represents the result of the system of equations for the homogeneous solution, which will be part of the general solution. 1.6. Nonhomogeneous solution In this step, the particular solution of the inhomogeneous system ( ph ) will be deduced. It is necessary to know the problem that is being proposed, in which the recharge and the extraction flow will be constant over the time intervals. If the equation for the particular solution is: [ ]{ } [ ]{ } { }p pS h K h Q+ =& (23) Using the undetermined coefficient method to arrive at the value of ( ph ) and knowing that the recharge vector {Q} will also be a constant vector in the considered time interval, the solution one comes up with is [17]: 1 2 ( ) ( ) constant vector over time. ( ) n p p p p h t h t h h t               (24) American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 557 The solution is a constant vector over time, meaning that its first derivative will be null [17]: = constant vector 0 p ph h  (25) We can substitute this result in equation (23): 1[ ]{0} [ ]{ } { } { } [ ] { }p pS K h Q h K Q-+ = \ = (26) Showing that the particular solution of the non-homogeneous equation can be written by equation (26). 1.7. General solution of the Differential Equation Finally, we have the general equation of the inhomogeneous system as the sum of the homogeneous and the particular solutions = +c ph h h , where ch is the solution for the homogeneous system (general solution) and the ph is the particular solution of the non-homogeneous system: 1.8. System decoupling The modal superposition method cannot be applied directly to coupled systems, according to equation (4). Thus, finding a coordinate system that does not exhibit any form of coupling is the essence of this procedure. Once equation (4) is uncoupled, the system of equations can be solved independently. The coordinates that allow the uncoupling of the system's equations are called principal bases, or normal bases [13]. To decouple a system of differential equations, it is necessary to introduce a set of alternative bases [13]: where equation (4) can be transformed into a set of n uncoupled equations, that is, their solutions can be independently determined. The normal or main basis {q} is defined by a transformation matrix [13] where [L] is the modal matrix (n x n) determined by solving an eigenvalue problem. ( ) 1[ ] { }i n i t i i h L c e K Q   (27) 1 2{ } { ( , , , )}nq q h h h (28) { } [ ]{ } { } [ ]{ }h L q h L q   (29) American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 558 In order to change a system of equations from the physical base to the normal base, it is necessary to multiply the modal matrices of the transformation in equation (4), according to equation (30): where, [ ]TL is the transpose of the modal matrix. The uncoupled system in the normal or main basis is obtained by substituting equation (28) in equation (29) [13]: which can be written as where [ ]S , [ ]K and [ ]Q are respectively the modal matrix of specific hydraulic storage, the modal matrix of conductance and the modal vector of recharge/extraction, given by [13]: At this point, we obtain the n uncoupled equations where: { } [ ][ ]T r r rS L S L= (37) { } [ ][ ]T r r rK L K L= (38) { } { }T r rQ L Q= (39) To transform to the physical bases, it is necessary to impose the initial conditions of the problem [13]: { (0)} [ ]{ (0)} { (0)} [ ]{ (0)}h L q h L q   (40) [ ] ([ ]{ } [ ]{ } { })TL S h K h Q+ =& (30) [ ] [ ][ ]{ } [ ] [ ][ ]{ } [ ] { }T T TL S L q L K L q L Q+ =& (31) [ ]{ } [ ]{ } { }S q K q Q+ =& (32) [ ] [ ] [ ][ ]TS L S L= (33) [ ] [ ] [ ][ ]TK L K L= (34) [ ] [ ] [ ]TQ L Q= (35) 1,2, ,r r r r rS q K q Q r n+ = =& L (36) American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 559 Multiplying equation (48) by [ ] [ ]T S we get: [ ] [ ]{ (0)} [ ] [ ][ ]{ (0)} [ ] [ ]{ (0)} [ ] [ ][ ]{ (0)} T T T T L S h L S L q L S h L S L q   (41) As [ ] [ ] [ ][ ]TS S= F F , one can thus obtain: [ ] [ ]{ (0)} [ ]{ (0)} [ ] [ ]{ (0)} [ ]{ (0)} T T L S h S q L S h S q   (42) In which the solution of the problem can be defined as 1 (0) { } [ ]{ (0)} 1, 2, , 1 (0) { } [ ]{ (0)} T r r T r r q L S h S r n q L S h S                    (43) 1.9. Uncoupled problem solution The solution to the modal superposition of the proposed problem is as described below: ( ) 1, 2, ,r n rr t r r rr r Q h l c e r n K                 (44) Where each term of the equations above is previously defined by equation (38) and equation (39) and for the value of the constant of integration cr it is obtained according to equation (45): { (0)} r r r r Q c q K æ ö÷ç= - ÷ç ÷÷çè ø (45) This completes the formulation that was implemented in a computer program in Matlab to solve the flow problem in porous media via modal superposition. In the next sections, applications of this program to a hypothetical case of flow in a confined aquifer will be presented. 2. Application 2.1. Underground water system As proposed in this paper, numerical codes were developed to compare the modal superposition method with the conventional finite element method. The quality and computational cost of the two numerical procedures were also compared. The computational was evaluated in terms of total CPU time and the results of the hydraulic head at specific points of the mesh. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 560 For the tests, it was necessary to simulate a hypothetical aquifer, considering aquifer configurations in the literature to guide the simulations, making the numerical experiment more realistic. In this section, the aquifer and the details of the simulation are described and the results of the simulations are presented. The CPU time and the discrepancies of results for meshes with different level of discretization were verified. 2.2. Simulated aquifer and meshes The aquifer has dimensions of 6000 by 6000 meters in length and width, composed of two reservoirs separated by a layer of low permeability (aquitard), as can be seen in Figure. 2. Its properties are shown in table 1. The upper receives a constant recharge of 81,0 ² m s m  . Although aquitard (material 2) has a lower permeability than aquifers (material 1), and is important, that it is not completely isolated, with the possibility of interaction between them. Figure 2: Aquifer material layers showing the different materials between aquitard and aquifers. Tabel 1: Permeability of materials. Material Permeability 1 Kx=1.8e-4 m/s Ky=1.8e-4 m/s Kz=1.8e-5 m/s 2 Kx=1e-15 m/s Ky=1e-15 m/s Kz=3.5e-9 m/s Discretization is a very important step in the numerical analysis of a given problem. The quality of the approximate solution obtained by the numerical method depends on the number of nodal points and elements that the problem domain is discretized. For the simulation, three levels of horizontal discretization of the domain were used. In order to verify the CPU time and the results for each type of discretization in both numerical methods presented here, as shown below: American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 561 Tabel 2: Types of Meshes. Mesh 1 Mesh 2 Mesh 3 Vertical discretization remained the same for all meshes, can be seen in Figure 2. The meshes have the same dimensions, only differentiating the number of elements and nodes, as can be seen in Table 3: Tabel 3: Node and element values for each mesh. Nodes Element. Mesh 1 6.724 28.800 Mesh 2 1.764 7.200 Mesh 3 3.844 16.200 2.3. Strategy adopted for the comparison of numerical solutions A well was considered at node 1, in all meshes with coordinates at the origin (x, y, z) = (0,0,0), as illustrated in Figure 3. Extraction rates of 30, 90 and 150 l/s were chosen so that their influence on each result can be studied. With the flows and recharge defined, monitoring points were proposed for the numerical simulation, the points are in Tabel 4. Figure 3: Well Location. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 562 Tabel 4: Monitoring points in the domain with the coordinates and the respective nodes for each mesh. Coordinates Chosen node data Point X Y Z Mesh 1 Mesh 2 Mesh 3 A1 3000 3000 90 4203 1103 2403 A2 3000 3000 60 2522 662 1442 A3 600 600 90 3531 927 2019 A4 600 600 60 1850 486 1058 A5 5400 5400 90 4875 1279 2787 A6 5400 5400 60 3194 838 1826 For each mesh and each proposed flow rate, the points in Tabel 4 were used to compare the results via conventional finite elements and by using the modal superposition technique. Each point has its fixed coordinate in the domain. As an example of the notation adopted for the results presented below, point A1 in mesh 1 with a flow rate of 30 l/s has the code 1-3A1 in the conventional finite element simulation and 1-3A1S for the simulation with modal superposition. Figure 4 shows the evolution of the values at point A4, where it is possible to observe the influence of the mesh discretization and the evolution of the simulations by the two methods. The values obtained in the two methods are very close, which was repeated for all points selected and listed in Tabel 4. Figure 4: Evolution of the hydraulic head at point A4 comparing the different discretizations of the meshes. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 563 Figure 5: Amplification of the rectangle in Figure 4 indicating the meshes used for each of the hydraulic at point A4. In Figure 6, for point A4, with meshes of different discretization and with different water extraction rates, the variation of hydraulic head over time can be observed. Once again, the consistency of the results is demonstrated, with very close solutions between the conventional finite element methods and the modal superposition. Figure 6: Time evolution of the hydraulic head at Point A4, for different discretization and imposed flow rates. An important variable to be observed is the total CPU time of the simulations. These values can be seen in Figure 7, which demonstrates the CPU time for the simulations in classical finite element method and another with modal superposition. The legend is referring to the type of mesh, (shown in Table 2), and the flow rate Mesh 1 Mesh 3 Mesh 2 American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 564 imposed on the well (30, 90 and 150 l/s) Figure 7: CPU time for simulations with the classical finite element method, and the modal superposition (with eigenvalues and superposition steps). As shown in Figure 7 the CPU time for both processes simulations used in this paper. In case of the modal superposition code process is performed in two steps:  The first step is the construction of matrices and the vector belonging to the domain ([S], [K] and {Q}) of the problem equation (4), and the calculation of the eigenvalue [L]. The eigenvalue matrix is an essential component of the modal superposition formulation described in section 2.5 and also for the system described in section 2.8. This step is represented in Figure 7 the figure legend by AV.  The second step is the calculation of the modal superposition system and the system decoupling procedure described in section 2.4. This step is represented in Figure 7 by SM. In order to compare the CPU time of the two numerical methods, the sum of the two steps described above was used as the total CPU time using the modal superposition, that is, the sum of the time spent for processing the eigenvalue matrix (AV) plus the processing time for calculating the modal superposition and that of the system decoupling (SM). It is observed that in meshes with greater discretization the processing time increases. This can be explained by the rise in the dimensions of the matrices of the algebraic system, requiring more CPU time to acquire the eigenvalues. For less discretized meshes, the CPU time was much smaller in the two processing steps of the modal superposition. The modal superposition method proved to be very useful for the type of problem proposed, where splitting the processing showed an advantage of the method, in which the part of the equation where the recharge is found can be modified in the simulation without the need to process the eigenvalue step, which is the one with the longest CPU time. In situations where the mesh does not need to have a greater discretization, the method saves a lot of CPU time. In cases of greater discretization (results of mesh 1, in figure 7), there is an increase in the American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 565 processing time of the superposition method in relation to the conventional finite element method. In the superposition method, the eigenvalue calculation step is crucial and the increase in discretization will result in an increase in the total CPU time. The method also allows a reduction of modes, reducing the modal matrix and consequently reducing the time in the eigenvalue processing step, thus resulting in a better performance of the code, that is, a reduction in the order of the problem is achieved with this step. The result of this order reduction can be shown with the results of Figure 8 and Figure9, which demonstrate the test performed to compare the behavior of the classical finite element method with that of modal superposition. In these figures, the discretization of mesh 1 was used, as shown in Tabel 2, and well flow rate of 30 l/s from the point of known coordinates A4 (Tabel 4). Figure 8: Comparison of CPU time. Figure 9: Comparison of the solution of the methods. Figure 8 presents the results for mesh 1 with a well flow of 30 l/s, with the discretization and flow already informed previously. It compares the processing time of the methods where in the first legend bar 1-3 is the processing time for the conventional finite element method and the other subsequent ones, with the legends bar 0 20 40 60 80 100 120 1-3 1-3 S 1-3 S 2 S ec o n d s CPU Time American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 566 1-3S and 1-3S2, are the processing time for the modal superposition method in which the normal processing of the method is represented and the other using the order reduction with 50% of the number of eigenvalues, respectively. It is shown that with modal superposition without order reduction, the processing time for the simulation used is the one that takes the longest time among those compared in Figure 8. But if the results of the first two bars are now compared with the third bar, which represents the superposition method with the order reduction of 50% of the eigenvalues, there is a reduction of almost 40% in processing time. It is important to note that the results found with this 50% reduction are shown in Figure 9, where the largest discrepancy found in relation to the conventional finite element solution was 0.18%. Showing that, even with a drastic order reduction of 50% of the system, the solution converges to the expected values, contributing to the reduction of the processing time and consequently with the decrease of the computational cost of the method. The codes were run on a computer with an Intel(R) Core (TM) i7-7700HQ CPU @ 2.80GHz, with 16.0GB of RAM. 3. Conclusion This paper was able to compare the flow solution in transient porous medium between conventional finite element and modal superposition methods, showing the versatility of the superposition method for simulations compatible with those demonstrated in the study carried out. The method has shown its potential for optimization studies and uncertainty analysis, where many simulations are needed, always looking for the best combination of parameters to maximize or minimize the design variables of a certain problem or to evaluate the impact of the variables that control the behavior of the aquifer during the water extraction process. In the previous sections, it has been shown that the method can be divided into two steps. Initially, it is necessary to set up the governing equation of the problem and consequently calculate the eigenvalue. The next step can be performed for any set of flows and recharge imposed on the aquifer and represents less than 50% of the CPU time for the analyzed cases. This method also allows order reduction, decreasing the number of eigenvalues needed to solve the problem, reducing the CPU time by 50% without loss of precision in the results, as demonstrated by the comparison with the classical finite element method. Acknowledgements This research work was carried out with the support of the Foundation to Support Science and Technology of the State of Pernambuco - Brazil (FACEPE). The first author thanks the Graduate Program in Civil Engineering (PPGEC) at the Federal University of Pernambuco (UFPE), which provided the necessary knowledge for the development of this research. References [1]. Bear, J. (1972). Dynamics of Fluids in Porous Media. Dover Publications. American Academic Scientific Research Journal for Engineering, Technology, and Sciences(ASRJETS)(2023) Volume 90, No 1, pp 548-567 567 [2]. Castilho-Barbosa, I. N., Carvalho, A. M., Miguel, G. F., & Carneiro, C. D. (2020). Conhecer para conservar aquíferos: como torná-los visíveis? Terrae didatica, 16, 1-12. [3]. Maidment, D. R. (1993). Handbook of Hydrology. McGraw-Hill. [4]. Guedes, R. P., de Araújo, M. P., & de Andrade, A. G. (2021). Necessidade do gerenciamento dos recursos hídricos em grandes cidades como Recife. Architecton - Revista de Arquitetura e Urbanismo, 06, 107-117. [5]. S. F. R. Khadri e C. Pande, “Ground water flow modeling for calibrating steady state using MODFLOW software: a case study of Mahesh River basin, India,” Model. Earth Syst. Environ., pp. 2- 17, 2016. [6]. W. H. Schilders, H. A. van der Vorst e J. Rommes, Model Order Reduction: Theory, Research Aspects and Applications, Berlin: Springer-Verlag Berlin Heidelberg, 2008. [7]. H.-j. G. Diersch, “Interactive , graphics-based finite element simulation of groundwater contamination processes,” Advances in Engineering Software, vol. 15, pp. 1-13, 1992. [8]. BEAR, J., & VERRUIJT, A. (1987). Modeling Groundwater Flow and Pollution. Springer Netherlands. [9]. HUYAKORN, P. S., & PINDER, G. F. (1983). Computational Methods in Subsurface Flow. ACADEMIC PRESS. [10]. Zienkiewicz, O. C., & Morgan, K. (1983). Finite Elements & Approximation . [11]. Logan, D. L. (2016). A First Course in the Finite Element Method - Sixth Edition. Cengage Learning. [12]. Zienkiewicz, O., & Taylor, R. (2000). The Finite Element Method (5 ª Ed. ed., Vol. 1 ). Butterworth Heinemann. [13]. Tedesco, J. W., McDougal, W. G., & Ross, C. A. (1998). Structure Dynamic : Theory and Applications. Prentice Hall. [14]. He, J., & Fu, Z.-F. (2001). Modal Analysis. Butterworth-Heinemann. [15]. Roy R. Craig, J., & Kurdila, A. J. (2006). Fundamentals of Structural Dynamics. John Wiley & Sons. Inc. [16]. Paz, M., & Kim, Y. H. (2019). Structural Dynamics. doi:https://doi.org/10.1007/978-3-319-94743-3 [17]. Zill, D. G. (2018). A First Course in Differential Equation with Modeling Applications (11th ed.). Boston: Cengage Learning. [18]. Greengerg, M. D. (1998). Advanced Engenieering Mathematics (2nd ed.). Prentice Hall.