Proposing mathematical model for seawater intrusion phenomena in the coastal aquifer

Seawater intrusion is one of groundwater quality problem which in this problem, the mixing between freshwater and saltwater in the coastal aquifer occurs. Mathematical modelling can be formulated to describe the mechanism of this phenomena. The main objective of this research is to develop the mathematical model of groundwater flow and solute transport that applicable to seawater intrusion mechanism. This mechanism is arranged as a differential equation and distinguished into 3 equations. The first equation is groundwater flow equation in dependent-density. It means that the density of groundwater (ρρ) changes in spatial and temporal domain due freshwater and seawater are mixed in the coastal aquifer. The second equation is solute transport. Like as groundwater flow equation, in solute transport equation there is a change of solute concentration (CC) in the spatial and temporal domain. The last equation is the relationship between groundwater density (ρρ) and solute concentration (CC). Special case for the third equation, in which this equation is adopted from USGS Seawat model. The first equation and second equation are governed by Eulerian mass conservation law. The main theoretical consideration of governing groundwater flow equation is such as fluid and porous matrix compressibility theory, Darcy’s law for groundwater in motion theory and some properties of soil. In other hands, solute transport is involving advection transport and hydrodynamic dispersion transport. Hydrodynamic dispersion is arranged by diffusion Fick’s law and dispersion in porous media theory and it depends on transversal and longitudinal dispersivity. Using Jacob Bear’s theory which states that fluid density as temperature, concentration and pressure function, authors obtain three primary variables in this model. Those variables follow fluid density (ρρ), total head (h) and concentration (CC). In this model, isotropic and isobar condition is considered, hence fluid density (ρρ) is a function of concentration (CC) only. Finally, from this research, authors wish this mathematical model is applicable to modelling, describing and predicting seawater intrusion phenomena theoretically.


Introduction
As one of freshwater main supply, an aquifer has an important role as the main storage of groundwater, including coastal aquifer. In other hands, behind of that important role, the coastal aquifer has salinity problem, in which groundwater from coastal aquifer noticeable salty. That is caused by phenomena of the mixing between seawater and freshwater in the coastal aquifer. If groundwater balance is achieved in the coastal aquifer, freshwater will flow normally. It means freshwater will be disembogued into coastal or ocean. In other hands, the human activity which exploits groundwater in coastal aquifer using the pump as an example, will decrease freshwater head physically, therefore seawater which has the larger density than freshwater theoretically, will give pressure along aquifer. If this physical phenomenon occurs for a long time without re-balancing groundwater in the aquifer, seawater will replace freshwater in the aquifer. This condition is called seawater intrusion [1].
Seawater intrusion phenomena can be modelled as a mathematical model. These phenomena have two mechanisms follows groundwater flow and solute transport. This mechanism occurs in a single phase, therefore there is miscible between seawater and freshwater in the aquifer. The solute transport mechanism is transformed as groundwater density change trigger. It is because based on theory, solute concentration gives effect in groundwater density changes. If the model is considered in isotropic and isobar condition, groundwater density is the only function of concentration.
Mechanism of groundwater flow usually informing the changes of the total head (ℎ) of groundwater in steady condition (Example: Modflow Model). In other hands, the real condition is unsteady because there are the changes in groundwater density due to solute concentration. Therefore, there are three main variables in groundwater flow equation follows fluid density ( ), total head (ℎ) and concentration ( ). To obtain their values, three equations are necessary. The first equation is groundwater flow equations in dependent-density condition. Many researchers have developed many equations which related to this condition. Voss [2] has started developing Sutra (Saturated-Unsaturated Transport) model. His research obtained groundwater flow model in which water saturation (Sw), fluid pressure (p), fluid density ( ) and concentration (C) as primary variables. Another groundwater flow model was developed by Guo and Langevin [3] and the model was named as Seawat. Seawat model is one of the models which applicable to seawater intrusion case. Its primary variables follow fluid density ( ), total head (ℎ) and concentration ( ). Seawat was governed using fluid density as the function of pressure and concentration. The last example is Codesa model in which this model considers pressure head ( ), fluid density ( ) and concentration (C) as primary variables.
There are some examples of Sutra Model Application such as seawater intrusion modelling in Ernakulam Coast [4]. Hussain, Javadi and Sherif [5]  The objective of this paper is proposing another mathematical model especially groundwater flow model which coupled with solute transport and fluid densityconcentration relationship equation. Using fluid density as concentration function only or = ( ), this research will focus to obtain mechanism in alteration of groundwater head (ℎ), solute concentration ( ) and fluid density ( ). Authors wish can develop an alternative model dependent-density groundwater flow model which applicable to seawater intrusion problem.

Eulerian mass conservation law
The Eulerian mass conservation law is used as main consideration in governing equation of groundwater flow and solute transport. Mass conservation law can be written as From equation (1) we say that the changes of mass in time domain along control volume system is equal with net flux which across along control surface plane.

Fluid compressibility
Fluid compressibility ( ) is measure changes of volume when the substances are subjected changes in normal pressure [9]. In isothermal condition, equation (2) is a term of fluid compressibility ( ).
∀ is the volume of water and ∀ is the volume of pressure.

Porous matrix compressibility
In natural condition, porous media volume in the determined depth of aquifer has hydrostatic pressure (p) and external pressure ( ). Porous matrix compressibility ( ) is measure changes of porous media volume due to the changes in external pressure. Porous matrix compressibility ( ) written as follow [9].
∀ is the total volume of the porous matrix which arranged by volume of water ∀ and volume of void ∀ .

Groundwater motion
For condition groundwater in motion, Darcy's law is subjected to flux (q). Darcy's law is formulated as [10] = − ℎ (4) K is hydraulic conductivity, h is total head and l is the distance between two head points.

Advection transport
Advection transport is a constituent transport by mean velocity of the fluid [11]. Especially in porous medium fluid flow, porosity (n) is involved in advection transport equation. Bear and Cheng [12] proposed advection transport ( ) term as = (5)

Diffusion and hydrodynamic dispersion
Another mechanism of solute transport is molecular diffusion. Based on diffusion Fick's law, molecular diffusion can be defined as Dm is molecular diffusion coefficient and i is a spatial domain. Molecular diffusion in porous media is subject to attenuation due electricity effect and tortuosity ( ) [13]. Because of that condition, molecular diffusion of porous media can be defined as Rubin [13] proposed dispersion coefficient in porous media as the series follows Usually in porous media, two last terms of equation (9) are neglected due to large Pe value. Transforming a 1 as transversal dispersivity ( ) and a2 as longitudinal dispersivity ( ), dispersion in porous media can be defined follows in sequence Finally, hydrodynamic dispersion (Dh) is obtained by equation (12) ℎ = + D is the value of D L and D T depending on the vector of fluid velocity and spatial domain.
Equation (13) can be transformed using calculus rules and the result can be written as From equation (14), there are two equations distinguished as left-hand side and right-hand side equation. Using chain-rule theory, left-hand side equation can be rearranged The first term of right-hand equation (15) is developed using fluid compressibility and porous matrix compressibility theory. We assumed that the changes of water volume are affected by the changes of pressure volume and pore volume, hence The term of equation (16) is governed by fluid compressibility and porous matrix compressibility theory. Based on equation (2) and (3), we obtain ∀ is obtained from ∀ = ∀ + ∀ where ∀ = 0, because in the static system, the volume of the solid phase does not change. Consider = ℎ and assumed that and g are independent, finally, in equation (17) and (18) is transformed as is constant groundwater density, is the gravitational force, and ℎ is the term of the changes of total head. Considering, the volume of water ∀ is obtained from total volume ∀ which affected by porosity n factor, therefore ∀ = ∀ . Hence, we can transform equation (17) as Next, determining of in equation (18) is based on the illustration below If we assume that the external stress of porous media is neglected, then = − and equation (18) is transformed as Combining equation (20) and (21)  Groundwater motion is involved in right-hand equation (14). Using Darcy law based on equation (4) and assumed as isotropic condition, the right-hand equation (14) is transformed as

Solute transport equation
Solute transport equation is developed to obtaining the changes of solute concentration in the spatial and temporal domain. It starts from mass conservation law for solute transport [11].
Using calculus rules, equation (31) is transformed as In the left-hand equation, there is porosity parameter (n) because in porous media, the volume of water depends on porosity.
In equation (35) is sink and source term of solute concentration.

Equation of
The term of is adopted from Seawat model [3]. Based on their research, they state that

Hypothesis
Using the mathematical model that explained before, we wish the result of model approaching its theoretical behaviour. Theoretically, if groundwater is pumped, its piezometric head will decrease. In this case, seawater in which has larger density will give pressure into freshwater along aquifer. For a long time effect, this condition will give effect to groundwater quality due to miscible between seawater and freshwater. Because of that condition, solute transport and groundwater density is a function of the spatial and temporal domain, = ( , , ) and = ( , , ). If this condition occurred for a long time, seawater will replace freshwater in the aquifer. Theoretical behaviour of seawater intrusion can be described by figure (2)

Conclusion
From this paper, we obtain three main equation related to seawater intrusion. Those equations are equation (30), (37) and (38). Discretization of them can use finite difference method or finite element method. The main key in application to seawater intrusion problem is all of the equation have to couple in a simulation.
Sink and source term are related pumping rate (sink) and recharge rate (source) in equations. Using this term we should develop model not only for seawater intrusion problem but also solution for mitigation using the artificial recharge. Simulation this model to verify model capabilities is necessary. This paper is a research proposal to aim master engineering degree thesis in Civil Engineering Department, Universitas Indonesia. It's a continuation research from Ma'Ruffi Kurnia [14] in which last research was not simulated yet and in the present research, there is some revision especially in solute transport equation.
Further analysis is author will discretization this mathematical model using finite difference method. Alternating Direction Implicit (ADI) scheme will be used in that step. Simulation of the model will use a computer program with some benchmark problems and simple simulation.
Finally, special thanks to Herr Soeryantono, Ph.D and Dr. Dwinanti as author's supervisors who have given many idea and recommendation which is related to this research.