Discontinuous Finite Element Methods For Seawater Intrusion Problem

Seawater intrusion occurs widely in many coastal countries and regions, it is the complicated movement of saline water into a freshwater aquifer, which is generally caused by changes of natural water environment and the social and economical development. It is well known that the density of the freshwater is lower than that of the seawater, so the fresh water column would flow over the column of the seawater in the aquifer, and they will arrive dynamic equilibrium of pressures on the interface. In order to satisfy the demand of the water for people's lives and economical uses, the coastal counties and regions often mine the groundwater excessively, which will induce the problem of seawater intrusion. When fresh water is withdrawn at a faster rate than it can be replenished, the water table is drawn down as a result. The decline of the water table reduces the pressure on the interface, which leads to the destruction of balance. To reach the new balance of pressures, seawater is pulled into the freshwater aquifer. When the pumping well is contaminated owing to the seawater invasion, the water can not be potable and used for production, and this brings great damages to people's lives and social and economical development of the coastal countries and regions. So far seawater intrusion problem has been found in hundreds of areas in scores of countries and regions, thus it is widely concerned by international community, and is under active research and government.In recent decades, there is some literature on the mathematical model and its numerical approximation of seawater intrusion problem. Based on the as-sumption that the interface, between the two fluid is a transition zone, Henry [51] obtained the analytic solution for the steady flow in the case of homogeneous medium and simplified boundary condition. In 1970, Pinder [67] first developed the transition zone model of seawater intrusion, and they also proposed a finite element methods for Henry's model. Segol et al. [77,78] developed a cut-plane two dimensional finite-element model with the water table and concentration cou-pling with each other. Huyakorn et al. [54] proposed the groundwater equation which depend on concentration and the solute transport equation, and consid-ered finite-element model of three-dimensional case. However, the research results above which are obtained under some precise assumptions, cannot really reflect the seawater intrusion process. Early in 1990s, Xue et al. [85,86] designed the first three-dimensional finite-element seawater intrusion model. They considered the case that the fluid velocity varies with the concentration, but they did not give the numerical prediction for the protection projects avoiding seawater intru-sion; for detail see [47] and citations therein. Since the end of the last century, Yuan ct al. deeply studied the mathematical model of the seawater intrusion, which has the characteristic of convection-dominated diffusion and nonlinearity, and proposed characteristics finite element methods [88] and modified method of upwind with finite difference fractional steps procedure [89.90,91,92] for seawater intrusion to simulate the transport of the mixing fluid after the protec-tion projects is carried out. Similar to the methods developed for oil reservior problems in [89,90,91,92], these methods avoid numerical oscillation which is the disadvantage of finite difference or finite element methods. However, they also have drawbacks:excessive diffusion introduced by upwind methods affects the accuracy; the characteristic methods are inconvenient for nonlinear problems [45]. Furthermore, they cannot do well with the problems with discontinuous solutions. In our paper, after observing the characteristic of the coupled equa-tions, we developed schemes combining discontinuous Galerkin with mixed finite clement methods for seawater intrusion to overcome the drawbacks above.The first discontinuous Galerkin method was introduced in 1973 by Reed and Hill [69] for the neutron transport equation which is a time-independent linear hyperbolic equation, and since then discontinuous Galerkin methods for hyperbolic equation have been widely concerned and actively developed. Many different numerical methods were put forward by incorporating discontinuous Galerkin methods with other methods. For example, Johnson and Pitkaranta [55] developed methods for hyperbolic equation, which is obtained by combining discontinuous Galerkin methods with upwind method; Cockburn el al. incorpo-rated so called 'high resolution'techniques (such as upwind, numerical flux and slope limiters; see [57] to know the ideas for detail.) for finite difference or fi-nite volume methods into the the framework of discontinuous Galerkin methods, and proposed Runge-Kutta discontinuous Galerkin (RKDG) methods in their serial papers [24,23,22,20,25], which are the major development of discon-tinuous Galerkin methods for hyperbolic equations [21]. To handle elliptic and convection-diffusion equations, based on the ideas of Runge-Kutta discontinuous Galerkin methods, Cockburn and Shu [26] and Castillo and Cockburn [15] pre-sented local discontinuous Galerkin (LDG) methods. Also in 1970s, but indepen-dently, discontinuous Galerkin finite element methods for elliptic and parabolic equations (also called interior penalty (IP) methods) were proposed and studied [37,7,84,2]. In recent years, many new types of interior penalty methods are de-veloped, including symmetric interior penalty Galerkin (SIPG) methods [2], the method of Baumann and Oden [9,64], nonsymmetric interior penalty Galerkin (NIPG) methods [74,72,75], incomplete interior penalty Galerkin (IIPG) meth-ods [34], and so on; See summary literature [3,14,4,70] for more information.Because discontinuous Galerkin methods can conserve mass locally, be suit-able for hp adaptivity and approximate solution with discontinuities better com-paring with continuous Galerkin methods, they are often employed to deal with practical problems. Computational fluid dynamic is the most widely used field, for examples, incompressible [80,81] and compressible [16,30,31] miscible dis-placement problems, gas dynamic problems [17], semiconductor device simulation [17,60], Navier-Stokes problems [8,48,49,50,65], two-phase flow [43,41,42], inviscid compressible fluid [35], fluid in porous media [73,82] and poroelasticity media [66], and so on. Additionally, they are also applied to other physical prob-lems, for example elasticity [61,71]. Observing that the models mentioned above are mostly of convection-diffusion type, so the numerical schemes given in these literature are the combination of the two classes of discontinuous Galerkin meth-ods, which handle the diffusion terms with the interior penalty methods and the convection terms with discontinuous Galerkin methods for hyperbolic equations. Moreover, to approximate better, techniques for finite difference and finite vol-ume methods, such as upwind [55,79,52,81,82,31,70], slope limiters [70] and numerical flux [35,36], are brought into the framework of discontinuous Galerkin methods to deal with convection-diffusion, especially convection-dominated equa-tions.Mixed finite element (MFE) methods, which have experienced comprehen-sive and deep research since they were presented [11,68,62,63], are often applied to solve partial differential equations of elliptic and parabolic types [56,13,12]. To carry out mixed finite element methods, we first change the original equation into a mixed formulation, by introducing a new variable (certain linear combi-nation of some derivatives of the unknown function), then apply finite element methods to the mixed variational formulation. In this way the original and new variables can be obtained at the same time, and they are approximated both at the same optimal order. Thanks to this advantage, since 1908s mixed finite element methods are combined with other methods to solve coupled systems of partial differential equations, for instance the model for miscible displacement in porous media [38,44]. In recent years, along with the development of discon-tinuous Galerkin methods, this new methods and mixed finite element methods are connected together to handle complicated coupled systems, for instance, in-compressible [80] and compressible [16,30,31] miscible displacement problems, fluid in linear poroelasticity media [66]. Noting that the mathematical model of seawater intrusion is also a coupled system consisting of a parabolic equation and a convection-diffusion equation, we handle it with similar methods, that is, mixed finite element methods for water head equation and discontinuous Galerkin methods for concentration equation.Under the aborative guidance of Professor Hongxing Rui, the author has finished this dissertation. In this dissertation, the author deeply study the math-ematical model of seawater intrusion problem considered by Yuan et al. in the literature [88,89,90,91,92], and develop some schemes to solve this model. Mixed finite element method and discontinuous Galerkin finite element methods are used, and iterative method is also used in this paper since the mathematical model is a coupled system of partial differential equations. Using these numerical methods, the author propose continuous-in-time scheme and implicit discrete-in-time scheme. For these schemes, we theoretically analyze the approximation property of the schemes, and obtain a priori error estimate. In order to improve computation rate and efficiency, the author develop a semi-implicit schemes by modifying the fully implicit discrete-in-time scheme, and finally we give a do-main decomposition procedure to solve the modified schemes in parallel. The dissertation is divided into four chapters.In chapter 1, we do some preliminary work. Firstly, we give the mathe-matical model which is put forward and studied in [88,89,90,91,92]. The model is a coupled system consisting of a parabolic equation and a convection-diffusion equation. Two different boundary value conditions are imposed on the system, and corresponding schemes for them are provided in the following chap-ters. Secondly, some definitions and notation are given, including Sobolev space and broken Sobolev space and their norms, partition and its properties, jump and average of scalar and vector. Finally, some theorems and inequalities required in error estimates are given in the form of lemmas. These lemmas were mentioned and proved in many books and papers.In chapter 2, we develop a continuous-in-time scheme to solve the mathemat-ical model equipped with the first boundary value condition. Mixed finite element method are used in order to solve the parabolic equation, while discontinuous Galerkin finite element methods (more precisely, SIPG) are employed to handle the convection-diffusion equation, where the convection term is treated only by inner product. After giving the variation formulations of the two equations, we propose the iterative continuous-in-time scheme. Optimal error estimates are ob-tain both in L2-norm and H1-norm. Results in this chapter are mainly from the author's paper [58], which had been accepted by Journal of Systems Science and Complexity (SCI).For the continuous-in-time scheme, we investigate some methods to discretize time, including explicit method, implicit method and high-order methods. Their stability and convergence are studied briefly, and their merits and drawbacks are also shown. In chapter 3, we mainly consider the model equipped with the second bound-ary value condition. The parabolic equation is also treated by mixed finite element method. For the convection-diffusion equation, upwind discontinuous Galerkin finite element methods are used to approximate the convection term better. Ad-ditionally, by introducing parameters A andθ, we consider two different interior penalty methods:SIPG and NIPG and two time-discretization methods:back-ward Euler's method and Crank-Nicolson in an identical expression. Then we obtain the fully implicit discrete-time scheme. Both hp a priori error estimate and results of the designed numerical experiments show that the scheme is ef-fective and reliable. Results in this chapter are mainly from the author's paper [59].It is well known that fully implicit scheme for nonlinear problems is of low efficiency, while explicit scheme is conditional stable and the time increment is limited. To find the way out of the dilemma, we develop a semi-implicit scheme, where linear diffusion terms are treated implicitly, and nonlinear convection terms are treated explicitly. Theoretical analysis and numerical examples both indicate that the accuracy of semi-implicit scheme and fully implicit scheme are almost the same.In chapter 4, we investigate a domain decomposition algorithm based on discontinuous Galerkin finite element methods for parabolic problems. As the semi-implicit scheme given in the previous chapter can be regard as the fully implicit scheme for parabolic equation, the parallel methods in this chapter can effectively shorten its computational time. Taking a standard parabolic equation as an example, we develop the domain decomposition algorithm, and discuss its stability and convergence. Finally numerical examples are done to verify our conclusion. Results in this chapter come from the paper [76].