Preferential seawater intrusion in a heterogeneous alluvial aquifer that extends offshore: Insights from numerical modeling

The problem

Over one billion people living near coasts depend on groundwater (Ferguson and Gleeson, 2012; Zamrsky et al., 2024). When groundwater levels on land fall far enough, the hydraulic gradient at the coast slopes landward and seawater moves into the aquifer. This active seawater intrusion causes the most severe salinisation of onshore groundwater (Badaruddin et al., 2017; Werner, 2017).

The confined gravel aquifers of the Central Canterbury Plains supply drinking water to Christchurch. The city abstracts around 50 million cubic metres of groundwater a year (Stewart and van der Raaij, 2022), and confined-aquifer levels have been as low as 1.5 m below sea level (Scott and Wilson, 2012). However, there is no evidence of extensive seawater intrusion, most likely because fresh groundwater extends offshore (Knight et al., 2018).

The Canterbury gravels contain open framework gravel deposits, which make up 1% to 20% of the aquifer (Dann et al., 2008; Burbery et al., 2018; Sarris et al., 2018). Such spatial variation in hydraulic conductivity, called heterogeneity, has been found to increase intrusion rates under pumping in fluvio-deltaic (Yu and Michael, 2019a) and volcanic (Geng and Michael, 2020) aquifers. However, no earlier study had examined its effect on active seawater intrusion in alluvial aquifers. In this study, we quantify that effect.

What we did

We built 25 realisations of a confined coastal aquifer based on the Central Canterbury Plains. Each has a 20 m thick aquifer beneath a 20 m thick aquitard. We distributed three hydrofacies with a transition probability method based on drillers' logs (Sarris et al., 2023). Clean gravel, which represents clusters of open framework gravel, makes up about 15% of the aquifer and has an average hydraulic conductivity of 500 m/day. Sandy gravel makes up 30% at 60 m/day. Clay-bound gravel and sand make up the remaining 55% at 4 m/day.

We simulated variable-density groundwater flow, in which density increases with salt concentration, with MODFLOW6 (Langevin et al., 2020), and salt transport with the same code (Hughes et al., 2017). Each run first held the onshore head 0.75 m above the offshore head until steady state. We then lowered the onshore head to 1 m below the offshore head for 100 years. We defined the interface at 1% of the salinity of seawater, the concentration at which groundwater may become undrinkable (Ministry of Health, 2018).

We compared the results with two equivalent homogeneous models, each with a uniform conductivity equal to the effective value of the heterogeneous aquifer, around 50 m/day. The advective model uses analytical solutions, with a sharp interface and no dispersion. The advective-dispersive model is a two-dimensional numerical model.

What we found

At steady state, the interface lay around 1,500 m offshore on average. In one realisation it came within 100 m of the shoreline, and in another it lay more than 2,500 m offshore.

Over 100 years of intrusion, the interface moved landward in every realisation. It reached the shoreline in half of them, and the inland model boundary in a quarter. Intrusion varied along the coast and appeared greatest along connected pathways of high conductivity. We hypothesise that connected areas of clean gravel and sandy gravel form these pathways. However, the average rate of interface movement did not correlate strongly with a percolation metric of how well the high-conductivity facies connect (Renard and Allard, 2013).

On average, the interface in the heterogeneous realisations moved two to four times faster than in the advective model, and up to twice as fast as in the advective-dispersive model. At its fastest point along the coast, it moved more than five times faster than the advective-dispersive model in many realisations. Methods that assess vulnerability to intrusion from hydraulic conductivity, such as GALDIT (Ferreira and Chachadi, 2005), should therefore also consider preferential flow in alluvial aquifers.

At a hypothetical monitoring well on the shoreline, salinity passed 0.35 PSU (practical salinity units) after 20 years in one realisation and 40 years in another, and in no other realisation until after 90 years. Denser coastal monitoring networks may be required to detect such intrusion (Werner et al., 2013).

The conductivity contrast between our most and least conductive facies is about 100 times. Earlier models of fluvio-deltaic (Yu and Michael, 2019a) and volcanic (Geng and Michael, 2020) aquifers used contrasts of 10⁹ and 10⁵ times. Preferential intrusion still occurred at our lower contrast. This indicates that heterogeneity should be considered in alluvial aquifers, even where the interface lies offshore.

Limitations

Our modelled aquifer extends 4 km offshore, while the aquifers it represents extend more than 30 km. In a longer aquifer, the steady interface could lie further offshore. It would then have further to move before salinising onshore groundwater. If preferential flow paths are shorter than 4 km, intrusion rates over longer distances may be lower on average. Our geological data come from around 40 km inland, where aquifers are more permeable and the depositional environment may differ. We induced intrusion by lowering the onshore head. We would expect more severe preferential intrusion under more focused abstraction, such as pumping.

Larger coastal models typically use cells about 10 times larger than ours (Thomas et al., 2019; Micallef et al., 2020; Michael et al., 2016). To reproduce preferential intrusion, these models need an upscaling method, a way of representing small-scale variation in conductivity within large cells. Dual porosity models may provide one, because they can better reproduce transport through preferential flow paths (Baratelli et al., 2011, 2014; Bianchi and Zheng, 2016).

Glossary

An aquifer is a layer of rock or sediment that holds and yields usable groundwater.

A confined aquifer is an aquifer overlain by a layer of low permeability, called an aquitard, that restricts flow into and out of it.

Alluvial sediments are gravels, sands and silts laid down by rivers.

Fluvio-deltaic sediments are those laid down by rivers where they enter the sea.

Hydraulic conductivity is a measure of how easily water flows through a sediment, here given in metres per day (m/day).

Effective hydraulic conductivity is the single value that gives the same total flow through a heterogeneous aquifer.

The head is the level to which groundwater would rise in a well, and the hydraulic gradient is the slope of the head.

Seawater intrusion is the movement of seawater into a coastal aquifer.

The salt/fresh water interface is the boundary between fresh and salty groundwater, defined here as 1% of the salinity of seawater.

Open framework gravel deposits are gravel deposits up to 100 times more conductive than the gravel around them.

A hydrofacies is a category of sediment with similar hydraulic properties.

A transition probability method builds a geological model from the probability that one hydrofacies is followed by another at a given distance.

A realisation is one of many equally likely versions of a model, each with a different arrangement of hydrofacies.

An analytical solution is a calculation made from mathematical formulas, without numerical simulation.

MODFLOW6 is a computer code that simulates groundwater flow and salt transport.

A steady state is a condition in which flow and salinity no longer change with time.

An equivalent homogeneous model represents a heterogeneous aquifer with a single, uniform conductivity.

Advective transport is the movement of salt with the average flow of groundwater.

Dispersion is the spreading of salt by small-scale variations in flow velocity and by molecular diffusion.

A preferential flow path is a connected route of high conductivity along which groundwater moves faster than through the surrounding sediment.

PSU (practical salinity unit) is a unit of salinity.

GALDIT is a method for assessing how vulnerable an aquifer is to seawater intrusion.

A dual porosity model divides an aquifer into a mobile, fast-flowing part and a less mobile part that exchange salt.