A Petrov-Galerkin finite element method for 2D transient and steady state highly advective flows in porous media

A new Petrov-Galerkin finite element method for two-dimensional (2D) highly advective flows in porous media, which removes numerical oscillations and retains its precision compared to the conventional Galerkin finite element method, is presented. A new continuous weighting function for quadratic elements is proposed. Moreover, a numerical scheme is developed to ensure the weighting factors are accurately determined for 2D non-uniform flows and 2D distorted elements. Finally, a series of numerical examples are performed to demonstrate the capability of the approach. Comparison against existing methods in the simulation of a benchmark problem further verifies the robustness of the proposed method.


Introduction
Many engineering problems involve heat or solute transport by means of both diffusion (or conduction) and advection. In those cases where the latter is dominant, the flow is referred to as highly advective. Examples where this phenomenon applies are ground source energy systems (GSES) which utilise the thermal energy beneath the earth's surface, such as open-loop GSES.
In recent years, considerable effort has been placed into researching the GSES systems due to their attractiveness as renewable energy sources. In many cases, numerical methods have been employed to investigate their behaviour and performance. However, finite element modelling of highly advective flows can often be problematic as the extensively used Galerkin finite element method (GFEM) may result in spatial oscillations of the nodal solution. Several studies (e.g. [1][2][3]) have shown that the oscillations become more significant with increasing Péclet number, which is a function of the properties of the porous medium governing diffusion and the velocity of the fluid, as well as the finite element size. Recently, Cui et al. [4] have also demonstrated that the type of the element and the boundary condition employed affect the behaviour and magnitude of the oscillations. Although the oscillations can be eliminated by refining the finite element mesh, in problems such as those involving open-loop GSES, this approach results in an extremely large number of elements, thus becoming computationally expensive.
The deficiencies of the GFEM have led to the development of various upwind finite element methods to model highly advective flows. The three main methods are ( [5]): artificial diffusion, quadrature and Petrov-Galerkin (PG). Artificial (or balancing) diffusion involves the introduction of an extra term to the physical diffusion in the advectiondiffusion equation which can then be solved using the standard GFEM. The main advantage of this approach is that it is easy to implement. However, it may result in a reduction of accuracy during the transient stage, as it alters the physics of the problem ( [3]). The quadrature technique was developed by Hughes [6] who suggested moving the quadrature (or integration) points within the element for a more efficient upwind effect. Conversely, the basic principle behind the PG formulation is to modify the nodal weighting functions in order to weigh the contribution from the upstream element more heavily than that from the downstream one.
The first upwind formulation for steady state advection-diffusion problems was proposed by Christie et al. [7] who modified the weighting function for a one-dimensional (1D) linear element. Later, the same approach was extended to two-dimensional (2D) linear elements by Heinrich et al. [8], as well as 1D and 2D quadratic elements by Heinrich and Zienkiewicz [9]. Huyakorn [10] suggested using different expressions for the weighting functions for 1D and 2D linear elements. Upwinding of 1D cubic elements was explored by Christie and Mitchell [11]. Kelly et al. [12] extended the technique to simulate a steady state problem with a uniform inclined flow in a 2D mesh by adding an extra diffusion term in the direction of flow to the advection-diffusion equation. Later, Donea et al. [13] explored a different modification to the steady state governing equation by studying quadratic elements. It should be noted that modifying the weighting functions has the same effect as adding artificial diffusion, provided only one-dimensional steady state is considered. Hence, most of the abovementioned methods https://doi.org/10.1016/j.compgeo.2018.04.013 produce exact solutions for steady state advection-diffusion problems.
The problem of transient advection-diffusion was found to be more challenging. The first solution to the transient problem for linear 1D and 2D elements was proposed by Huyakorn and Nilkuha [14] and Ramakrishnan [15]. While Huyakorn and Nilkuha [14] used the same continuous weighting functions for linear elements as Huyakorn [10], the time derivative term was weighted using the standard Galerkin weighting functions which resulted in oscillating and over-diffused solutions. Other researchers attempted to overcome the problem by developing higher order weighting functions. Cubic expressions were adopted by Ramakrishnan [15], who observed a reduction in accuracy despite applying the modified weighting functions to all terms in the advection-diffusion equation. Later, Dick [16] proposed weighting functions which are a combination of continuous cubic and quadratic expressions for 1D and 2D linear elements, whereas Westerink and Shea [17] extended this approach to quartic weighting functions for 1D quadratic elements. Cardle [18] applied different weighting functions to the time derivative and the spatial derivative in the 1D advection-diffusion equation. Another approach to the transient problem was suggested by Yu and Heinrich [19,20] who developed continuous weighting functions which include time dependency (the bilinear time-space shape functions). The same method with different expressions for the time-space shape functions was adopted by Al-Khoury and Bonnier [21].
Brooks and Hughes [5] developed the streamline upwind/Petrov-Galerkin (SU/PG) method where discontinuous weighting functions, instead of the continuous functions mentioned above, were used. In this approach, multi-dimensional flows were also considered. However, the weighting functions as well as the adopted weighting factors for multidimensional cases were formulated provisionally, and Brooks and Hughes [5] suggested further research should be carried out to obtain a more rigorous approach for multi-dimensional cases. Based on the SU/ PG approach, further development and modification to this approach were proposed in Hughes et al. [22] and Tezduyar and Park [23] (Shock Capturing (SC) method), Hughes et al. [24] (Petrov-Galerkin Least Square (PGSL) method), and Tezduyar and Ganjoo [25] and Codina et al. [26] (time-space weighting functions). Compared to the original SU/PG method, improved results were observed when multi-dimensional highly advective flows were simulated. However, for the abovementioned methods, only results of problems involving flow through a mesh with regular (e.g. square or rectangular) elements were shown, while the performance in problems using irregular finite element meshes has not been demonstrated in the literature.
It is evident that some of the methods mentioned above (both in one-and multi-dimensional forms) were observed to produce overdiffused solutions. This reduction in accuracy was automatically assigned by some researchers in this field to all upwind techniques (e.g. [27]). While this is true for the artificial diffusion approach, Brooks and Hughes [5] argued that the Petrov-Galerkin finite element method (PGFEM) does not experience this problem, provided appropriate modified weighting functions are employed.
In the current study, the PGFEM with continuous weighting functions is adopted and further developed. While examining the continuous weighting functions proposed in the literature, the authors concluded that the use of the more complex higher order expressions (e.g. [15,17]) is not necessary, and that the time-space shape functions (e.g. [19,20]) are not compatible with the time marching scheme employed in the current paper. Additionally, it was also found that the PG weighting formulation for 2D 8-noded quadratic elements proposed by Heinrich and Zienkiewicz [9] cannot eliminate oscillations in either steady state or transient analyses. In fact, in none of the abovementioned studies using continuous weighting functions (e.g. [8][9][10][11]) was the PGFEM successfully applied to both transient and steady state problems.
As a result, a consistent PGFEM with continuous weighting functions for solving the time-dependent advection-diffusion equation is proposed here, using both linear and quadratic elements. This scheme adopts the weighting functions of Huyakorn [10] and Heinrich and Zienkiewicz [9] for 1D linear and 1D quadratic elements, respectively, while the approach of Heinrich et al. [8] and Huyakorn [10] is adopted for weighting functions of 2D linear elements. However, a new formulation for 2D quadratic elements is developed and presented in this paper as the continuous formulation proposed by Heinrich and Zienkiewicz [9] was found to be inadequate. Furthermore, the modified weighting functions are applied to all terms in the time-dependent advection-diffusion equation which was not the case in some of the existing approaches (e.g. [14]). Additionally, a new numerical scheme for determining the Péclet numbers for 2D elements is also proposed, to ensure that the method retains its accuracy in complex cases, such as those involving 2D non-uniform flows (where both the flow direction and velocity vary spatially) and those where 2D distorted elements may be present. It should be noted that none of the PGFEMs available in the literature has been shown to perform well for these types of problems.
For the purpose of demonstrating its capabilities, the newly proposed PGFEM in this paper has been implemented into the finite element software ICFEP ( [28]), which is the authors' bespoke computational platform and therefore allows access to source code, which is crucial for this type of development. A series of numerical studies are then performed, demonstrating the accuracy and stability of the new formulation for both transient and steady-state cases at different levels of complexity. Comparisons were also made between the proposed method and other methods using discontinuous weighting functions, such as the SU/PG and the SC methods, in the simulation of a multidimensional case widely documented in the literature (e.g. [22,23]). Finally, the paper demonstrates the excellent performance of the proposed PGFEM in simulating challenging multi-dimensional flow problems involving distorted elements where both flow velocity and flow direction vary across an element, which have not been shown anywhere in the available literature.

Coupled thermo-hydraulic formulation
Numerical modelling of convective flows in porous media requires the simulation of pore fluid flow coupled with heat transfer through both diffusion and advection, which is generally referred to as a coupled thermo-hydraulic (TH) problem. For incompressible pore fluid flow in a fully saturated porous medium, the continuity equation must be satisfied, which can be expressed as: where v x , v y and v z are the components of the velocity of the pore fluid in the x, y and z directions, respectively, ε v is the volumetric strain due to stress changes, Q f represents any pore fluid sources and/or sinks, and t is time. The seepage velocity, z T , is assumed to be governed by Darcy's law, which can be written as: where k [ ] f is the permeability matrix and ∇h is the gradient of the hydraulic head. In a coupled thermo-hydraulic problem, if it is assumed that: (a) the effects of temperature gradients and spatial variations in fluid density (e.g. buoyancy-driven flows) on pore fluid flow through a fully saturated porous medium are negligible and (b) the solid phase is rigid, Eq. (1) reduces to the equation of seepage, which can be expressed as: Assuming an instantaneous temperature equilibration between the fluid and the solid phases, the heat transfer in a porous medium is governed by the law of conservation of energy ( [29]): where Φ T is the heat content of the porous medium per unit volume, Q T is the heat flux per unit volume including heat conduction and heat advection, Q T are any heat sources/sinks and dV is the volume of the porous medium. If the solid phase is assumed to be rigid, Eq. (4) reduces to the transient heat advection-diffusion equation: is the conductivity matrix of the porous medium, and ρ and C p are the density and the specific heat capacity of the mixture, respectively, which, for a fully saturated porous medium, can be expressed in terms of its components as: where n is the porosity, and the subscripts f and s denote the fluid phase and the solid phase, respectively. By combining Eqs. (3) and (5), a coupled thermo-hydraulic FE formulation can be obtained for modelling the convective flow in a rigid porous medium. It is noted that the solution of Eq. (5) depends on the fluid velocity obtained from Eq. (1), which may vary with both time and space in an analysis.

Boundary conditions
In a coupled thermo-hydraulic FE analysis, the flow of water can be generated by either prescribing pore water pressures (i.e. a Dirichlet boundary condition (BC)) or specifying an inflow/outflow at the boundaries (i.e. a Neumann BC). Similarly, heat flow may be simulated by applying either a prescribed temperature BC (i.e. a Dirichlet BC) or a heat flux BC (i.e. a Neumann BC). It should be noted that at the boundary where the pore fluid enters or leaves, either a prescribed temperature or a coupled thermo-hydraulic boundary condition is necessary ( [4]), where the latter applies a heat flux, q T , equivalent to the energy associated with the fluid flowing across the boundary: where v { } f b represents the velocity of the pore fluid, T b is the current temperature at the boundary, and B denotes either a line, a surface or a volume over which the boundary condition is prescribed. It is noted that this non-linear boundary condition ensures that there is no diffusive heat flux across the boundary and only advective heat flux is allowed.

Péclet number
In the context of heat transfer, the dimensionless Péclet number, Pe, represents the ratio between the advective and diffusive heat fluxes. It is therefore a useful measure of the significance of the problematic advective heat transport. A Péclet number above one indicates that advection is the predominant heat flux, whereas a Péclet number below one signifies a diffusion-dominated heat flow. Mathematically, the Péclet number is defined as: where L is the characteristic length in the direction of flow. In the case of the finite element method, L is defined as the element length in the direction of fluid flow.

Modelling highly advective flows with Galerkin FE method
The Galerkin finite element method (GFEM) has been successfully adopted to simulate numerous problems in engineering. However, it is widely accepted that modelling of highly advective flow is often problematic, as GFEM produces numerical oscillations when solving the advection-diffusion equation in cases where the Péclet number is large. Most studies on this topic (e.g. [1][2][3]) involved a one-dimensional transient advective-diffusive flow along a bar with Dirichlet boundary conditions at both ends. It was shown that the steady state solution experiences oscillations, the magnitudes of which increase with increasing Péclet number.
Recently, Cui et al. [4] investigated the problem associated with the GFEM further by examining the effect of the finite element type, as well as the chosen boundary conditions. It was demonstrated that quadratic elements experience smaller oscillations than linear elements with the same Péclet number. Additionally, the behaviour of the problem with a fixed temperature boundary condition was shown to be different from that with the coupled thermo-hydraulic boundary condition. Indeed, in the case where a fixed temperature is applied at the position where the fluid leaves the mesh, the oscillations appear once the heat front reaches the extremity of the bar, with their amplitude increasing until steady state is achieved, at which point they remain constant. Conversely, when the coupled thermo-hydraulic boundary condition is prescribed, the oscillations also appear as the heat front reaches the boundary, however, their amplitude initially increases and finally reduces such that, at steady state, the temperature is uniform along the bar. Based on these findings, Cui et al. [4] obtained limits on the values of Péclet number which ensure non-oscillatory solutions for combinations of element types and boundary conditions, mitigating the impact of the limitations of the GFEM.
In order to model a problem with specific material properties (ρ f , C pf , k T ) and fluid velocity (v f ), the Péclet number (Equation (8)) can only be reduced by refining the finite element mesh (i.e. decreasing L). However, in some cases, such as modelling of open-loop GSES, the process of refining the mesh is often problematic as the high fluid velocity results in a large number of extremely small elements, therefore increasing the computational effort significantly. This paper explores a new implementation of the Petrov-Galerkin finite element method in order to remove this obstacle to computationally-efficient modelling of highly advective flows.

Petrov-Galerkin FE method for highly advective flows
To eliminate the spatial oscillation of temperature when modelling highly advective flows ( > Pe 1) in porous media, a Petrov-Galerkin finite element method (PGFEM), instead of the conventional GFEM, is adopted and investigated here. As the spatial oscillations only exist in the solution obtained for the temperature field, the PGFEM is applied to derive the FE formulation based on the law of energy conservation (Eq. (5)), while the FE formulation derived from the continuity condition (Eq. (3)) remains in the unchanged Galerkin form. Discretising Eq. (5) and applying the divergence theorem leads to: where W [ ] is the matrix of weighting functions, N [ ] is the matrix of temperature interpolation functions (or shape functions), the prime denotes the spatial derivatives, Ω represents the domain over which the integral is performed, and B represents the boundary of the solution domain. The conventional GFEM assumes while the PGFEM employs weighting functions which are different from the shape functions. To solve Eq. (9) the θ -method ( [28,30]) has been adopted as the time marching scheme in this paper. It should be noted that relevant modifications based on the PGFEM should also be made to the formulation of the residual (or out-of-balance) heat flux for the iterative solution process, where the weighting functions are no longer equal to the shape functions.

Linear elements
Huyakorn [10] proposed a PGFEM for solving 1D steady state advection-diffusion problems with the 2-noded linear isoparametric elements shown in Fig. 1. The weighting functions were defined as: where the subscripts denote the nodal number, N 1 and N 2 are the isoparametric nodal shape functions, s represents the natural ordinate of the 1D elements, α is the PG weighting factor and function f s ( ) was chosen as: The PG weighting factor, α, controls the level of upwinding. Its optimal value depends on the Péclet number and can be calculated as: Huyakorn and Nilkuha [14] attempted to apply the above PG weighting function to the solution of the 1D transient advection-diffusion equations. However, the original shape functions N [ ], instead of the new weighting functions W [ ], were used for the first term on the left-hand side of Eq. (9) which involves the derivatives of temperature with time. Therefore, the numerical examples they presented exhibited oscillations and over-diffused solutions.

Quadratic elements
The PG weighting functions for solving 1D steady state advectiondiffusion problems with the 3-noded quadratic elements shown in Fig. 2 were proposed by Heinrich and Zienkiewicz [9] as: where β 1 and β 2 are the PG weighting factors and the function g s ( ) was defined as: The optimal value of the PG weighting factors β 1 and β 2 can be obtained from: It should be noted that the above PGFEM was developed for steady state problems and its behaviour in solving 1D transient problems has not been shown elsewhere in the literature. In this paper, this PGFEM has been applied to Eq. (9), resulting in the PG FE formulation for transient advection-diffusion problems with quadratic line elements.

4-noded linear elements
The PG weighting functions for the 2D 4-noded linear element shown in Fig. 3 can be obtained by combining two 1D 2-noded linear element weighting functions ( [8,10]). However, the two 1D weighting functions now depend on the S and T natural coordinates, such that: It must be noted that the PG weighting functions shown above have only been used for steady state advection-diffusion problems in the literature (e.g. [8,10]). Following a similar approach to that described for the 1D element, a 2D PG FE formulation is derived by combining Eqs. (9) and (19)(20)(21)(22). This PG scheme requires the definition of the values of the Péclet number in the S and T directions, which will be discussed in later sections of this paper.

8-noded quadratic elements
Heinrich and Zienkiewicz [9] proposed the PG weighting functions for solving steady state advection-diffusion problems with the 8-noded quadratic elements shown in Fig. 4. However, when this PG scheme was   implemented and tested as part of the research presented in this paper, the obtained results were clearly unsatisfactory, exhibiting spatial oscillations even in the simplest case of 1D advection-diffusion with constant and uniform fluid velocity. Therefore, in this paper, new PG weighting functions are proposed for 8-noded quadratic elements based on the standard Garlerkin form of shape functions: It should be noted that with the new PG weighting functions shown above it is possible to simulate highly advective flow using a mesh which contains a mix of both 8-noded elements and 4-noded elements. The performance of the new PG weighting functions in both steady state and transient analyses will be shown in subsequent sections of this paper.

Determination of weighting factors for 2D finite elements
As explained above, the PG weighting functions for 2D elements depend on the optimal values of the weighting factors (α for linear elements and β 1 and β 2 for quadratic elements), which are formulated based on the Péclet numbers in the S and T local coordinate directions, Pe S and Pe T . For the proposed PG scheme, Pe S and Pe T were firstly evaluated at each Gauss point, which, according to Eq. (8), can be expressed as: = Pe where the subscripts S and T denote the directions of local coordinates with respect to global coordinates x and y, which can be defined by the vectors r { } S and r { } T : Furthermore, the unit vectors in the S and T directions (as shown in Based on the transformation from global to local coordinates, the velocities in the S and T directions, v fS and v fT , can be derived as: In this paper, the conductivity is assumed to be isotropic which implies that: , . In situations where the global conductivity is anisotropic, the transformation can be found in the Appendix A. It should be noted that, when the derived Pe S and Pe T at each Gauss point are directly used to evaluate the weighting factors, the PG scheme described above leads to accurate results when simulating 1D highly convective flows, with either constant or variable fluid velocities, using a 2D FE mesh with regular (square or rectangular) elements. However, invalid results were found when either multi-dimensional flows were simulated or distorted elements were used (see example 4 of 1D flows and the example of 2D flows shown later). Further research demonstrated that the performance of the proposed PG scheme when analysing such complex problems is drastically improved if elemental average Péclet numbers of Pe S and Pe T , as well as elemental average unit vectors of e S and e T , are used to determine the weighting factors. These are denoted Pe S , Pe T , e { } S and e { } T herein to distinguish these average quantities from those evaluated at each Gauss point. It should be noted that a minimum Péclet number, Pe min , is employed here (equal to 1.0 for linear elements and 2.0 for quadratic elements as proposed in Cui et al. [4]), below which Pe S and/or Pe T are set to zero, meaning that the weighting functions in the corresponding direction become identical to the element shape functions, as used in the GFEM. Additionally, for situations in which there is a variation in the directions of local axes S and T (e.g. when distorted elements are used) and/or there is a variation in the 2D advective flow directions, a correction factor, ω, was shown to be needed in order to avoid under-weighting the contribution of the less dominant component of the flow and obtain accurate results: Therefore, the modified elemental * Pe S and * Pe T , which are used to determine the weighting factors, are defined as: In the following sections, the performance of the PGFEM formulation as described above, is illustrated in a series of examples of increasing computational complexity, to which different components of the proposed method have varying degrees of importance. For example, in the simpler applications, the use of element-average Péclet numbers and correction factor ω, while active, does not affect the weighting functions. In effect, to demonstrate the accuracy resulting from the adoption of these two aspects of the formulation, specific examples featuring distorted elements (Example 5) and water flow with the spatiallyvarying direction (Example 6) were employed.

Numerical examples of 1D and 2D flows in a mesh with regular elements
The PGFEM described in Section 3 is implemented and applied to the modelling of 1D and 2D highly advective flows in a mesh with regular elements in this section, in order to validate the proposed approach for transient problems and explore its behaviour under various complex conditions (e.g. 1D flows with non-uniform velocities and uniform 2D flows). The finite element software ICFEP is employed in all analyses.

Example 1: 1D flow with uniform velocity
In order to validate the formulation for the proposed PGFEM for transient problems, several plane strain analyses of a problem involving 1D highly convective flow illustrated in Fig. 6 were performed with both 2D linear and quadratic elements. The advective-diffusive flow of heat from left to right in the x-direction, as indicated by the arrow, was simulated by applying pore fluid pressures, u, at both sides of the mesh (u 1 > u 2 ) to induce a fluid flow with a constant and uniform velocity, as well as a heat source in the form of a fixed temperature of 10°C at the left boundary. At the right boundary, either a fixed temperature of 0°C (i.e. no change from the initial temperature) or the coupled thermohydraulic BC were prescribed. Flows with Péclet numbers of = Pe 10 S , Fig. 6. Geometry of the analysed problem with initial and boundary conditions (Example 1).   Table 1. The θ-method with = θ 2/3 was adopted throughout. Figs. 7 and 8 show the nodal temperature distribution from left to right in the mesh with 4-noded linear elements and 8-noded quadratic elements, respectively, and the fixed temperature BC at the right boundary, whereas Figs. 9 and 10 show the results of the equivalent analyses with the coupled thermo-hydraulic BC at the right boundary. For comparison, the solution of the same analyses performed with the GFEM are also plotted on the figures. Cui et al. [4] observed that, in the case of the GFEM, the oscillations of nodal temperature occur only once the heat front reaches the right boundary of the mesh. Before that stage (i.e. the transient stage), the solution obtained using the GFEM is exact, and it can therefore be used to verify the accuracy of the PGFEM. The results of the transient stage (which was chosen at an arbitrary time and was the same in both the GFEM and PGFEM analyses) shown in Figs. 7 and 8 indicate that the PG formulation presented in this paper does not lead to any reduction in accuracy, independent of the magnitude of the Péclet number. Moreover, this accuracy of results highlights the importance of applying the proposed PG weighting functions to all terms of the transient advection-diffusion equation: (i) Huyakorn and Nilkuha [14], who applied the same linear weighting functions (i.e. Eqs. (1013)) but without modifying the time-dependent term, did not achieve accurate solutions when compared to those obtained by the GFEM in  Figs. 7 and 8; (ii) Equally, the PG formulation for 8-noded quadratic elements employed by Heinrich and Zienkiewicz [9] also leads to oscillatory solutions, which is not the case for the formulation proposed here, as shown in Fig. 8.
As shown in Cui et al. [4], the behaviour of the oscillations with the GFEM depends on the boundary condition. In the case of the fixed temperature, the oscillations increase once the heat front reaches the mesh boundary until a steady state is achieved with maximum amplitude of oscillations. Figs. 7 and 8 show that this problem is eliminated by employing the PG formulation. An exact steady state solution is obtained, which is a uniform temperature of 10°C at all nodes of the mesh with the exception of the right boundary, where a drop to 0°C satisfies the prescribed temperature BC. When the coupled boundary condition is applied, the GFEM produces oscillatory solution as the heat front reaches the right boundary, however, with time, the oscillations reduce and finally a steady state with a uniform temperature is obtained. Again, the results computed with the PGFEM are free of oscillations at all stages of the analysis, as shown in Figs. 9 and 10. These findings imply that the effectiveness of the newly proposed PG method is independent of the magnitude of the Péclet number, as well as the boundary conditions employed. It should also be noted that both reduced integration and full integration can be adopted when the proposed PGFEM is used to simulate 1D highly convective flows with uniform velocities.

Example 2: 1D flow with non-uniform velocity
This example is performed to investigate the performance of the proposed PG formulation in cases where the fluid velocity varies spatially. Here, a series of axisymmetric analyses was performed on the cylinder shown in Fig. 11. Similar boundary conditions as in Example 1 were prescribed to generate an advective-diffusive heat flux from right to left (the opposite direction to Example 1). Due to the axisymmetric nature of the problem, the fluid velocity increases along the radial direction with a minimum at the right boundary and a maximum at the left boundary. As all 400 elements in the mesh are of equal size, the Péclet number, which is directly proportional to the fluid velocity, increases in the same manner (from right to left) and is plotted in Fig. 12. Again, the two different boundary conditions were explored.
In this application, various combinations of linear and quadratic    Table 2. Although all combinations produced excellent results in Example 1, the problem involving non-uniform fluid velocity proved more challenging for cases 1, 2 and 3 which resulted in incorrect temperature distributions during both steady state and transient stages. Conversely, Case 4 produced the exact non-oscillatory solutions shown in Fig. 13 together with the corresponding GFEM results. The findings of this study suggest that the orders of the shape functions used to calculate the fluid velocity (i.e. the pore fluid pressure shape functions) and the PG weighting functions used to obtain the nodal temperature, must be the same to ensure accurate results.

Example 3: 2D flow with uniform velocity
The behaviour of the new PGFEM when simulating 2D uniform flow was studied by performing plane strain analyses with 2D linear elements. The mesh with dimensions of 10 m × 10 m and 20 elements in each direction, as well as the applied boundary conditions, are illustrated in Fig. 14     applied hydraulic and thermal boundary conditions, with both GFEM and PGFEM. Oscillations were observed during both transient (the time of which was chosen arbitrarily but was the same for all analyses) and steady state stages when GFEM is adopted, while they were eliminated by employing the PG formulation.

Comparison with existing methods
A set of benchmark tests involving highly-advective flows skew to the mesh, similar to Example 3 in Section 4.3, were adopted by Hughes et al. [22] and Tezduyar and Park [23], who claimed that these were very difficult problems which could be used to test the robustness of a method. Following the description presented in Hughes et al. [22] and Tezduyar and Park [23], the same tests were performed using the PGFEM proposed in this work. A mesh with dimensions of 10 m × 10 m and 20 square linear elements in each direction was employed. Linearly distributed pore fluid pressure boundary conditions were appropriately applied to induce flows which are 30°, 45°and 60°skew to the mesh respectively. Heat transfer is dominated by advection here in the sense that the element Péclet number is extremely large (both Pe S and Pe T are well in excess of 10000). For each case, as described in [22] and Tezduyar and Park [23], a fixed temperature of 10°C was prescribed at all the nodes on the bold line, while a fixed temperature of 0°C was prescribed at all the remaining boundary nodes (see Fig. 16). Fig. 17 shows the steady state temperature distributions in the mesh obtained using the proposed PGFEM. To better compare to the results with the existing methods shown in the literature (e.g. Hughes et al. [22]), 3D views, by plotting the nodal temperature solutions illustrated in Fig. 17 as a third dimension (vertical), are plotted in Fig. 18. This figure can be directly compared with the equivalent figures given in the literature [22], which are reproduced in the Appendix B. It can be seen that results obtained using the proposed PGFEM were significantly more accurate when compared to those by other existing methods (e.g. SU/PG and SC). It is noted that Tezduyar and Park [23] further modified the methods proposed by [22]. However, compared to the results shown in [22], limited improvement was observed in Tezduyar and Park [23].     An analysis of 1D convective flow, which is similar to the case with a Péclet number of 100 in Example 1, was carried out. However, a 2D mesh using elements with varying levels of distortion, as shown in Fig. 19, was used. Since distorted 2D elements are adopted in the mesh, local directions of S and T vary throughout the mesh, resulting in  significantly varying Péclet numbers Pe S and Pe T before the algorithm is applied. Based on the conclusion from Example 2, quadratic pore fluid pressure shape functions were used with linear temperature shape functions in this study. Fig. 20 shows the transient and steady state temperature distribution in the mesh for both the analysis using the proposed algorithm (i.e. * Pe S and * Pe T ), and one which only uses the average elemental Péclet numbers Pe S and Pe T without the correction factor of ω (Eq. (39)). It is noted that significant oscillations (over ± 50°C) were observed for both transient and steady state in the latter analysis. Fig. 21 compares the horizontal transient and steady state temperature distributions obtained in Example 1 and Example 4 (the analysis using * Pe S and * Pe T ). The very good agreement between the analysis using regular elements and distorted elements demonstrates the validity of the proposed algorithm for distorted meshes.
6.2. Example 5: 2D flow in a distorted mesh A more complex case with 2D convective flow and Péclet numbers of around 100 was studied. The 2D mesh using distorted elements shown in Fig. 22 was used. Constant pore pressure and temperature boundary conditions were prescribed at the top and bottom boundaries, while no water flow and no heat flux boundary conditions (i.e. boundaries were assumed to be both impermeable and insulated) were prescribed at the two lateral boundaries. This resulted in a 2D highly advective flow with varying direction and magnitude of fluid velocity as shown in Fig. 22. Similar to the previous examples, quadratic pore fluid pressure shape functions and linear temperature shape functions were adopted here. Fig. 23 compares the transient and steady state temperature distributions obtained in analyses with the proposed algorithm for 2D flow (i.e. using * Pe S and * Pe T ) and with the average elemental Péclet numbers without the correction factor of ω. Adequate heat transfer was observed for the first case, while large oscillations (over ± 100°C) were produced in the latter case.

Conclusions
This paper presents a new Petrov-Galerkin finite element formulation with continuous weighting functions which enables the accurate modelling of 2D transient and steady state highly advective flows in porous media. The scheme has been successfully implemented into a finite element software and applied in a series of numerical analyses to demonstrate the superiority of the proposed formulation. The ensuing conclusions are summarised as follows:  (1) The weighting functions of Heinrich et al. [8] and Huyakorn [10] for 2D 4-noded linear elements are found to be adequate, but need to be applied to all terms of the time-dependent advection-diffusion equation.
(2) The weighting functions for 2D 8-noded quadratic elements developed by Heinrich and Zienkiewicz [9] are found to produce oscillatory solutions when solving steady state advection-diffusion problems. The newly proposed weighting functions in this paper allow the simulation of highly advective flow using meshes where either only 8-noded elements, or a mixture of 4-noded and 8-noded elements are used. (3) A new algorithm is proposed to determine the weighting factors used for establishing weighting functions for 2D elements, accounting for both the variation in directions of local axes when distorted elements are used, as well as the variation in directions of 2D convective flows. (4) The results of numerical analyses using the proposed PG formulation demonstrate the ability of this method to eliminate the oscillations characteristic of the GFEM solutions and to accurately simulate both transient and steady state highly advective problems.
The proposed approach is also capable of dealing with 2D uniform, as well as varying fluid velocities in 2D meshes composed of both regular and distorted finite elements. This extends the applicability of the proposed methodology to the modelling of various problems encountered in geotechnical engineering, such as the open-loop GSES. (5) It is also highlighted that attention must be paid to the modelling approach in order to obtain accurate solutions. The authors recommend the use of higher order shape functions for the calculation of pore fluid pressures than for temperatures. (6) The robustness the proposed PGFEM has been demonstrated by comparing the obtained results against other exiting methods in the simulation of a widely used benchmark problem. Fig. A1. Comparison in temperature distribution of the benchmark test between different methods: (a) exact solutions from [22]; (b) SU/PG from [22]; (c) SC1 from [22]; (d) SC2 from [22]; (e) proposed PGFEM.