COMSOL 6.4 - Seawater Intrusion in a Coastal Aquifer

Seawater intrusion is a significant concern in coastal regions, where the balance between freshwater sources and saltwater infiltration is constantly under threat. This problem is complex because it not only threatens the supply of fresh water but also the quality of groundwater.

This example illustrates how to set up a model for saltwater intrusion in a coastal aquifer when a pumping well is positioned at some distance from the shoreline. The inspiration for this approach arose from the research paper titled “Preferential Flow Enhances Pumping-Induced Saltwater Intrusion in Volcanic Aquifers.” In this research, a range of methods were utilized to replicate the process of saltwater intrusion in a volcanic aquifer distinguished by the presence of highly conductive “lava tubes”, contrasted with other geological formations displaying comparatively lower conductivity levels.

Given the distinctive characteristics of this scenario, where a low-conductivity aquifer intersects with highly conductive tubes, our model incorporates both the homogenized porosity approach as well as the dual porosity approach to capture these conditions.

Model Definition

Figure 1 illustrates the model’s geometry and scenario, with the sea level represented by a hydraulic head condition. The model assumes a consistent rate of freshwater recharge, and the land surface is subjected to precipitation.

Figure 1: Model geometry and conditions.

Freshwater has a density of ρf = 1000 kg/m3, while saltwater has a higher density of ρs = 1025kg/m3. These density differences are essential and are incorporated into the model by considering the transport of salt within the aquifer and utilizing a linearized density relationship:

(1)

Saltwater typically has a salinity of approximately 35% or 35 g/l. In this model, the specific value of cs (maximum salt concentration) is not critical, as all effects are contingent upon the relative salt concentration c. Nonetheless, realistic salt concentration values are utilized by assuming a molar mass of 58.44 g/mol, resulting in cs = 598.9mol/m3.

The transport through the porous aquifer is governed by the conservative formulation of the diffusion–convection equation which also includes dispersion.

(2)

here, Dd and De are the dispersion and diffusion coefficients., εp is the porosity , and on the right hand side Qs denotes the source term. The convective velocity u is calculated using Darcy’s law:

(3)

together with the continuity equation

(4)

Here, K denotes the hydraulic conductivity, g and g are the gravity constant and vector, respectively, while Qm represents a source term. Note that the density ρ is given by Equation 1.

In Ref. 1 the researchers conducted a comparison of different approaches, including a homogenized approach and various heterogeneous approaches to represent diverse facies with distinct conductivities. They specifically set the conductivity for lava tubes to be two orders of magnitude larger than the next highest facies, sparking the idea of investigating a dual porosity approach.

In this example, a homogenized approach is initially considered assuming an anisotropic hydraulic conductivity where horizontal conductivity exceeds vertical conductivity. Then the solution is compared to that of a dual porosity approach.

Dual porosity implies that the flow occurs within the highly conductive area only, namely the lava tubes, which constitute the macroporous part of the system. In contrast, flow within all other facies is stagnant due to their significantly lower conductivity, forming the microporous system. However, these facies can still exchange fluids with the tubes, acting as sources or sinks. Consequently, Equation 3 and Equation 4 are employed for the volume fraction of the macropores (θM). The source term Qm represents the interporosity flow and depends on the pressure difference between the macro- and micropores

(5)

The fluid transfer function αw (s/m2)is influenced by various factors such as the facies structure and characteristics, and it is typically not known with accuracy. A more detailed discussion can be found in Ref. 2. Smaller values of αw indicate longer fluid exchange times. Therefore, when observation periods are sufficiently long and a steady state regime is achieved, the precise value becomes less significant. For the volume fraction of the micropores only an ordinary differential equation (ODE) needs to be solved:

Above equations are incorporated using the feature within the Darcy’s Law interface. To account for the dual porosity characteristics in salt transport, this is manually integrated by introducing an additional ODE for the micropores:

in conjunction with a mass source term within the transport equation for the macropores (Equation 2) with the mass source being:

Like for the fluid transfer function αw, the mass transfer function αs (1/s) is not known and discussed in Ref. 2.

Results and Discussion

The model is calculated over a duration of one year. Notably, the concentration distribution remains relatively stable after approximately 100 days. The entire year is simulated to account for varying precipitation rates (Figure 2), although it demonstrates that their impact is relatively minor.

Figure 2: Varying precipitation rate over one year.

Figure 3 shows the pressure.

Figure 3: Pressure distribution after one year.

The concentration plot shows the salt concentration for the dual porosity approach. A contour line marks the location of the interface between saltwater and freshwater, where the concentration is 17.5 g/mol, representing half the concentration found in saltwater in the homogenized approach.

Figure 4: Concentration distribution and total concentration flux (arrows) after one year for the dual porosity approach. The black contour line indicates the position of the freshwater–saltwater interface for the homogenized approach.

The outcome of the dual-porosity simulation is that even though the volume fraction of the highly conductive lava tubes is small compared to the aquifer’s low-conductive constituents, they can have a large impact on the saltwater intrusion characteristics. The results show that the dual-porosity approach predict a larger saltwater intrusion as compared to the homogenized approach.

Note that the parameters in this example are chosen arbitrarily. To predict saltwater intrusion in real aquifers it is essential to have good experimental data of the aquifer composition.

References

1. X. Geng and H.A. Michael, “Preferential flow enhances pumping-induced saltwater intrusion in volcanic aquifers,” Water Resour. Res., vol. 56, 2020; doi.org/10.1016/j.jhydrol.2022.127835.

2. T. Vogel, H.H. Gerke, R. Zhang, and M.Th. Van Genuchten, “Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties,” J. Hydrol., vol. 238, nos. 1–2, pp. 78–89, 2000. doi.org/10.1016/S0022-1694(00)00327-9.

Subsurface_Flow_Module/Solute_Transport/seawater_intrusion

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Start by loading parameters and the geometry sequence into the model.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file seawater_intrusion_parameters.txt.

Geometry 1

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file seawater_intrusion_geom_sequence.mph.

In the toolbar, click

Definitions

Next, define a variable for the density to account for density changes due to salt concentration.

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
rho	rhof+(rhos-rhof)*c/cs	kg/m³	Seawater density

Darcy’s Law (dl)

In the window, under click .

In the window for , locate the section.

Select the checkbox.

Gravity 1

The sea level is at y = 0m. To accurately account for the influence of gravity effects, the reference position is modified to align with the surface elevation for locations above sea level (y > 0) and is set to 0 for locations below sea level (y < 0).

In the window, under > click .

In the window for , locate the section.

Select the checkbox.

Specify the rref vector as