water
Article
Imaging the Structure and the Saltwater Intrusion Extent of the
Luy River Coastal Aquifer (Binh Thuan, Vietnam) Using
Electrical Resistivity Tomography
Diep Cong-Thi 1,2, *, Linh Pham Dieu 1,2 , Robin Thibaut 1 , Marieke Paepen 1 , Huu Hieu Ho 2 , Frédéric Nguyen 3,4
and Thomas Hermans 1
1
2
3
4
*
Citation: Cong-Thi, D.; Dieu, L.P.;
Thibaut, R.; Paepen, M.; Ho, H.H.;
Nguyen, F.; Hermans, T. Imaging the
Structure and the Saltwater Intrusion
Extent of the Luy River Coastal
Aquifer (Binh Thuan, Vietnam) Using
Electrical Resistivity Tomography.
Water 2021, 13, 1743. https://
doi.org/10.3390/w13131743
Academic Editor: Eyal Shalev
Received: 19 May 2021
Accepted: 21 June 2021
Published: 24 June 2021
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims in
Laboratory of Applied Geology and Hydrogeology, Department of Geology, Ghent University,
Krijgslaan 281-S8, 9000 Ghent, Belgium; Linh.PhamDieu@Ugent.be (L.P.D.); Robin.thibaut@Ugent.be (R.T.);
Marieke.paepen@Ugent.be (M.P.); Thomas.Hermans@UGent.be (T.H.)
Vietnam Institute of Geosciences and Mineral Resources (VIGMR), Hanoi 100000, Vietnam;
hohuuhieu@yahoo.com
Urban and Environmental Engineering Research Unit, Liege University, Allée de la Découverte 9,
4000 Liège, Belgium; f.nguyen@uliege.be
Department of Earth and Environmental Sciences-Geology, KU Leuven, Celestijnenlaan 200e-bus 2411,
3001 Leuven, Belgium
Correspondence: Diep.CongThi@UGent.be
Abstract: With the growing population and the adverse effects of climate change, the pressure on
coastal aquifers is increasing, leading to a larger risk of saltwater intrusion (SI). SI is often complex
and difficult to characterize from well data only. In this context, electrical resistivity tomography
(ERT) can provide high-resolution qualitative information on the lateral and vertical distribution of
salinity. However, the quantitative interpretation of ERT remains difficult because of the uncertainty
of petrophysical relationships, the limitations of inversion, and the heterogeneity of aquifers. In
this contribution, we propose a methodology for the semiquantitative interpretation of ERT when
colocated well data are not available. We first use existing wells to identify freshwater zones and
characterize the resistivity response of clayey deposits. Then, we approximate the formation factor
from water samples collected in the vicinity of ERT data to derive a resistivity threshold to interpret
the saline boundary. We applied the methodology in the shallow aquifers of the Luy River in
the Binh Thuan province, Vietnam, where water resources are under pressure due to agricultural,
aquacultural, and industrial production. Twenty-one ERT profiles were collected and revealed a
much larger intrusion zone, compared to the previous study. Saltwater is present in lowland areas of
the left bank over almost the whole thickness of the aquifer, while the right bank is constituted of
sand dunes that are filled with freshwater. At a larger distance from the sea, a complex distribution
between fresh and saltwater is observed. Our methodology could be applied to other heterogeneous
aquifers in the absence of a dense monitoring network.
Keywords: saltwater intrusion; aquifer; resistivity; salinity; geophysics
published maps and institutional affiliations.
1. Introduction
Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
Saltwater intrusion (SI) occurs along shores and results from the interaction between
seawater and coastal aquifers through different processes, such as tidal mixing, water
density effects, surface–groundwater interaction, and variable recharge patterns [1]. The
complex combination of potential natural and anthropogenic impacts has placed coastal
aquifers under intense pressure by salinity contamination. On one hand, the natural state of
the system is often influenced by anthropogenic activities, such as groundwater extraction
or irrigation [2,3]. In many cases, several processes take place, and it might be very difficult
to discriminate the dominant processes [4].
4.0/).
Water 2021, 13, 1743. https://doi.org/10.3390/w13131743
https://www.mdpi.com/journal/water
Water 2021, 13, 1743
2 of 24
On the other hand, geological conditions, such as the presence of confining layers combined with paleohydrological conditions, can also yield complex saltwater distributions
related to ancient seawater trapped in sediments [5]. In addition, geological heterogeneity
also is identified as one of the most important factors affecting the spatiotemporal dynamics
of saltwater in coastal aquifers [1,6]. It not only affects the evolution of SI through creating
preferential flow paths but also makes the conceptualization of SI more difficult by preventing the generalization of the processes. One additional challenge for the management
of water resources is to collect data at a relevant scale. Often, monitoring well networks at
a regional or catchment scale is not sufficient to obtain a comprehensive view of SI, and the
relatively large distance between boreholes makes it difficult to characterize the underlying
spatially dependent processes.
Indeed, in geologically heterogeneous systems, data collected in wells might not
be sufficient to reveal the complexity of SI, both vertically and laterally. Geophysical
techniques offer a cost-effective alternative, by providing indirect information on the
aquifer properties at a high spatial resolution [7]. Geophysical data can be used as a
characterization tool but can also be more systematically integrated into the development
of groundwater conceptual models and the calibration of numerical models [6]. Electrical
and electromagnetic methods are particularly well suited for investigating SI because of
their sensitivity to the bulk electrical conductivity of the ground, which is directly related
to the salinity of pore water [8]. If airborne methods are becoming more popular to control
large areas [9,10], they generally require an investment of local authorities to cover the
large mobilization costs [11]. In addition, the airborne methods might not provide the
high-resolution necessary to capture shallow processes and are more challenging to apply
in a time-lapse manner for monitoring purposes.
An efficient alternative is to use electrical resistivity tomography (ERT) [12]. Numerous
applications have demonstrated the ability of ERT to characterize the salinity in coastal
aquifers at various spatial and temporal scales [9,13–16].
Nevertheless, the qualitative interpretation of ERT data requires the integration of
direct observation for validation, especially in heterogeneous systems. The presence of
conductive clay sediments can not only hide the presence of saltwater but can also span
the same range of magnitude as brackish water, making the unequivocal identification
of salinity more difficult [14,17,18]. Quantitative interpretation is even more challenging.
In a recent contribution, González-Quiros and Compte [6] studied the ability to delineate
salinity boundaries from inverted ERT data. They identified three main sources of errors
limiting the direct interpretation of ERT. First, the choice of the petrophysical relationship
is essential to link the bulk electrical resistivity with the salinity [19]. The petrophysical
relationship depends on many parameters such as the porosity, the cementation factor
(often grouped in the formation factor), and the clay content, which also influences surface
conductivity processes [20]. Despite recent efforts [19], there is no unified petrophysical
model that could be broadly applied without a necessary calibration step. This makes
it virtually impossible to quantify the salinity in absence of colocated data, i.e., ERT
data collected at the location of existing boreholes and wells to validate petrophysical
relationships. Second, the difficulty is further enhanced by the regularization used for the
ERT inversion, often leading to over smoothing [6] and the loss of resolution of surface
ERT at depth [14,15,21], which affects the petrophysical relationship [22,23]. It often
leads to discrepancies when comparing inverted results with borehole data [14,24,25].
The salinity is often overestimated in the low-salinity zone and underestimated in the
high-salinity zone [6]. This effect can be reduced but not eliminated by integrating crossborehole [25] or surface-to-hole measurements [14] or by choosing more appropriate
regularization [21]. Third, the heterogeneity of the deposits affects the ability to derive
accurate salinity estimation. As noted earlier, petrophysical relationships are affected by
the clay content, so that there is not a unique relationship between salinity and resistivity
in heterogeneous settings, introducing some errors if this is not accounted for during
interpretation [6]. However, accounting for this effect would require knowing the geological
Water 2021, 13, 1743
3 of 24
model in advance, which is impossible at the local scale and prone to uncertainty at a larger
scale [19]. This is not only important for salinity estimation, but heterogeneity also controls
the evolution of SI through its effect on hydraulic properties [26,27].
The above-described elements illustrate the difficulty to quantitatively interpret ERT
results in terms of salinity. An ideal study would integrate colocated data and laboratory
analysis to derive the petrophysical relationship, analyze the effect of spatially varying
resolution, and improve the inversion through incorporating prior information. In addition,
a geological model should be known to account for the heterogeneity. These elements
are rarely available in real-world field studies. There is, therefore, a strong need for
simple methodologies allowing a robust interpretation of ERT data. As an alternative, we
propose a two-step procedure for the qualitative characterization of salinity for aquifers
characterized by a high degree of heterogeneity in the presence of scarce colocated data
or in their absence. First, we identify freshwater zones from available water samples in
the aquifer to assess the resistivity response of clay lenses. Second, we use water samples
collected in the vicinity (but not colocated) of ERT profiles to approximate petrophysical
relationships and define a conservative, semiquantitative limit for the saline boundary.
We apply this methodology to derive the saline boundary in the coastal aquifer of the
Luy River in the Binh Thuan province (Vietnam). This aquifer is characterized by a large
heterogeneity, with the presence of numerous clay layers and lenses. Although shallow and
dug wells are available, agricultural constraints make it difficult to collect colocated data.
This region is thus particularly challenging for the interpretation of ERT data; however, our
approach can be broadly used in other study areas with high heterogeneity and limited
colocated measurements.
The paper is organized as follows. First, we present the geology and hydrogeology
of the study site. Then, the details of the methodology, including the ERT data analysis,
inversion, and depth of investigation estimation, are displayed. Finally, the results of our
two-step procedure are presented to assess the saline boundary and interpret the observed
salinity distribution.
2. Geology and Hydrogeology of the Study Site
Binh Thuan is one of the most arid provinces in Vietnam [28,29] with an average
monthly precipitation under 50 mm and a monthly evaporation of 140 mm [30] in the dry
season from November to April of the following year [28]. Therefore, this province is subject
to cyclic droughts putting the water resources under pressure. A further factor making
Binh Thuan vulnerable to SI is characterized by its long coastline of 192 km and seven river
systems flowing into the sea through estuaries affected by the tides with a magnitude in
the range from 0.8 to 2 m [16,31,32]. The Luy River catchment is also particularly impacted
by this phenomenon [26], as saline water was observed up to 10–15 km inland [33].
Understanding the dynamics of salt- and freshwater is thus essential but challenging
for the management of the coastal aquifer [8]. Previous investigations mapped the transition of saline and freshwater into two different areas based on shallow boreholes in the
downstream part of the Luy River. Two disconnected zones along the river were detected
(Figure 1) so that the river was identified as the main source for SI [34]. However, a recent
study estimated that seawater could penetrate about 7 km inland from the sea during the
highest tide [35]. This could not explain the presence of saltwater as far as 15 km from the
sea and the occurrence of several saline zones apparently disconnected from each other.
The distribution of salinity is, therefore, more complex than expected.
Water 2021, 13, 1743
4 of 24
Figure 1. Location of study area (red rectangle) and 21 collected ERT profiles (green lines; the yellow
points indicate the origin of the profile with the 1st electrode at 0m) and boreholes (red circles). The
red solid line indicates the limit between the upstream and downstream parts of our study area. The
two purple solid lines correspond to saline boundaries with a TDS limit = 1.5 g/L and 1.0 g/L from
well measurements obtained by NAWAPI (2015) [34]. Topographical isolines (in m above sea level)
are shown by brown curves.
The alluvial plain of the downstream part of the Luy River is a part of the Binh Thuan
coastal zone characterized by a U-shaped cross section and relatively homogeneous topography inclining gently toward the sea [31,36]. The study area has been governed by a part of
the Mesozoic continental active margin commonly called Da Lat Zone [36]. This is proven
by geological evidence consisting of the calc-alkaline magmatic arcs developed during the
late Jurassic–late Cretaceous, such as the Deo Ca complexes of I-type granitoids [36–38],
the Ankoet complex sourced from a mixture of I and S-type granitoid [39], and the Nha
Trang extrusive formation (K nt) [36–38,40] (Figure 2a). In the Mesozoic and Cenozoic, the
study site was affected by two perpendicular normal fault systems of NW-SE and NE-SW
strikes [16] and subparallel fault (Figure 2a). The NW-SE faults were considered to be
active in the late Mesozoic and reactivated in the Cenozoic, creating the Holocene and
Pleistocene coastal sediment plains and the Neogene–Quaternary sediment basins which
partly or totally overlay the older geological formations. The NE-SW faults are related to
Neogene–Quaternary basaltic eruptions that created the shorelines in the same direction as
the fault strike [37].
Sandy dunes are widely distributed on the right bank of the Luy river and along the
rest of the Binh Thuan coastal zone. They were formed by complex geodynamic systems
changing over time and depended on tide regime and sediment supplies. Generally, coastal
−
2a) from the
sand dunes and bars are part of the Phan Thiet formation (mb Q1 2−3 pt) (Figure
Pleistocene period originating from marine and marine-eolian depositional sources and the
weathering alteration of older bedrocks [38,41]. The sediment sequences are characterized
by ridges of sand composed mostly of quartz but bearing some heavy minerals such as
ilmenite, zircon, anatase, rutile, pyroxene, and monazite [36].
Water 2021, 13, 1743
5 of 24
(a)
(b)
Figure 2. (a) Geological map of the downstream part of the Luy River and its vicinity. Three fault systems (the dotted lines)
are covered by Quaternary sediments. (b) Two geological cross sections: AB and CD are approximatively parallel and
perpendicular to the shoreline, respectively.
On the left bank, the flat plains have a low topographical elevation, including tributaries of the Luy River [34,41]. The sediments originating from complex alluvial-marshmarine depositional cycles during the Holocene and Pleistocene are narrowly distributed
and adjacent to the white sandy dunes and low hills [42] (Figures 1 and 2). The sediment
sequence is characteristic of this depositional system with gritty sandy layers intercalated
with sand-rich clay or thin clay lenses that can locally acts as impermeable and confining
layers (Figure 2b). The mineralogical composition is made mainly of quartz, feldspar, mica,
and small rocky fragments from the bedrocks, such as granite, basalt, rhyolite, and dacite.
These Holocene and Pleistocene sediment sequences also correspond to the unconsolidated aquifers constituting the main groundwater resource in the area. The Holocene
aquifers (Q2 1−2 , Q2 2 , and Q2 2−3 , Figure 2b) have a thickness varying from 2 m to 20 m
on the left bank and increasing towards the sea [16,27,35]. However, they are considered
to be totally or partly absent on the right bank because of strong erosion and uplift processes [16,34,42,43]. Inversely, the Pleistocene aquifers (Q1 2 , Q1 2−3 , and Q1 3.2 , Figure 2b)
have a low thickness on the left bank, while they are relatively thick within the sand dunes
on the right bank. The recharge of the unconsolidated Holocene aquifer depends directly
on the resources of surface and rainwater. The storage capacity is limited and often insufficient for irrigation [16,34,42]. From previous studies, the Holocene aquifer is expected to be
mostly salinized by saltwater intrusion in the estuarial part and coastal low land. The soil
layers have also been degraded and salinized by agricultural activity [16,32,42]. In contrast,
the abundant storage capacity and productivity of the Pleistocene aquifer is proven by the
pumping test [30,32,34,42] and is considered to play an important role in irrigation and
industrial production activities. Depending on the presence or absence of clay layers and
lenses, the Holocene and Pleistocene aquifers can be separated by a confining layer or lying
on each other, in which case they cannot really be discriminated (see borehole description
in Section 5). Because of the tectonic activities, the bedrock formations are locally fractured
and can contain local aquifers; however, their capacity is limited.
Water 2021, 13, 1743
6 of 24
3. Methodology
3.1. Data Acquisition
Based on the existing saline boundary (Figure 1) identified through water samples
collected from shallow wells at a depth ranging from 3 m to 8 m [34], 21 ERT profiles
were collected using an ABEM® Terrameter LS. A multiple-gradient array [44] with a
separation factor s = 8 and an electrode spacing parameter a = 1, 2, 3, 4, 5, 6 was used for all
profiles. This multiple-gradient configuration was chosen to accommodate the high contact
resistance observed in most profiles due to dry soil conditions while maintaining a good
signal-to-noise ratio and a good compromise between vertical and lateral resolution [45].
The electrodes were buried and sprayed with saltwater to decrease the contact resistance.
The maximum injected current was limited to 200 mA to accommodate the high contact
resistance observed in most profiles due to these dry soil conditions.
Each profile was originally composed of 64 electrodes but was extended afterward
using the roll-along principle when possible [46]. An electrode spacing of 5 m was used to
image the whole thickness of the unconsolidated aquifers with sufficient resolution (see
Sections 3.4 and 5.3).
3.2. Data Quality
The recorded raw data were first analyzed for their error level. Data quality was
assessed using the contact resistance of the electrodes and the repeatability error. The latter
was obtained by repeating the same measurement sequence and observing the coefficient
of variation (expressed in %). All points with a coefficient of variation higher than 1% were
removed from the datasets prior to inversion. This threshold was chosen based on the
overall quality of the data.
The data quality was also assessed using reciprocal error [7,47,48]. Due to field
conditions, we focused on gathering one profile per day, and collecting reciprocal for long
profiles was too time demanding. Therefore, we have reciprocal for 8 out of the 21 profiles.
The reciprocal error is defined as
eN/R = RN − RR
(1)
where RN and RR are the normal and reciprocal resistance values, respectively. Slater et al. [47]
proposed an error model where this error increases linearly with the mean R of normal
and reciprocal resistances:
eN/R = a + bR
(2)
where a represents the minimum absolute error, and b defines the relative increase of the
error with the resistance. These parameters are determined by the envelope in an error
versus resistance plot that contains all the points after the removal of outliers.
Globally, the reciprocal error confirms high data quality. Except for a few outliers,
most data points have a relative reciprocal error Re lower than 5%, with the vast majority
even showing an error below 1% (Figure 3a).
Following the methodology of Slater et al. [47], we investigated if a relationship exists
between the error and the mean resistance (Figure 3b). For most of the profiles, no increase
of the reciprocal error with the resistance is observed. Globally, the reciprocal error ranges
between 0 and 0.001 Ohm, with similar values observed for the whole range of measured
resistance. Fitting the envelope error model of Slater et al. [47] therefore results in a very
low relative error (<1%) and an absolute error in the order of magnitude of 0.001 Ohm. A
constant error model would actually be suited in this case.
Water 2021, 13, 1743
7 of 24
ோ
(a)
(b)
Figure 3. (a) Histogram of the relative reciprocal error for representative profile 29. The large majority of points have a
reciprocal error lower than 1%. (b) A model of the reciprocal error for profile 29; an increase of the error with the resistance
is not observed. Similar distributions are observed for other profiles.
3.3. Inversion
To recover the electrical resistivity distribution m from the collected resistance data
d, an inverse problem must be solved. In ERT, m and d typically cover several orders of
magnitude and are often expressed in a log scale ensuring the positivity of the recovered
resistivity. From a distribution of electrical resistivity, the corresponding dataset can be
computed by solving the nonlinear forward problem noted:
d = f (m)
(3)
where f is the forward model operator.
The inverse problem is solved using a regularization approach, during which we
minimize an objective function of the form:
ψ(m) = ψ(d) + λψ(m)
(4)
where the first term on the right-hand side of the equation expresses the data misfit, the
second term represents the model functional with some assumed characteristic of the
model, and λ, the regularization parameter, balances these two terms. The second term
stabilizes the inversion process, in the sense that the sought model must verify some a
priori properties [49].
This objective function can also be expressed as
ψ(m) = kWd (d − f (m))k p + λkWm mkp
(5)
where Wd defines the data weighting matrix accounting for varying data quality; Wm
represents the model constraint matrix, often chosen as the roughness matrix, although
other options exist [21]; and the suffix p refers to the L p norm.
We use the software Res2Dinv [50] to retrieve the resistivity distribution of the subsurface explaining the collected data. The L1-norm is used on both the data to reduce
the potential effect of outliers and on the model to favor sharp resistivity contrasts. The
solution of the inverse problem is obtained by minimizing this objective function using an
iterative Gauss–Newton scheme, starting with an initial model m0 taken as the average
apparent resistivity in the data set [50].
The topography is included in the finite-element mesh to avoid any undesirable
effect of the surface on the recovered resistivity distribution. The topographical profile
Water 2021, 13, 1743
8 of 24
was extracted from the digital elevation model from the USGS (United States Geological
Survey) and corrected with GPS data and field observations.
We considered that the inversion process converged once the data misfit reached
a value below 3%, which was the case for all profiles. This is coherent with the error
estimated from reciprocal measurements.
3.4. Depth of Investigation Index
ERT data collected from the subsurface are more sensitive to shallow layers than to
deeper layers. The decrease of sensitivity and resolution with depth is well documented [51]
and must be accounted for to avoid the overinterpretation of resistivity structures. This
is even more crucial in saline grounds, as electrical flow lines tend to concentrate in high
conductivity zones, which can significantly reduce the depth of investigation [8].
We follow the methodology of Oldenburg and Li [52] to calculate the depth of investigation index (DOI). The DOI is computed based on a dataset inverted with two
different reference values (ρre f 1 and ρre f 2 ) yielding two different inversion models (ρ1 and
ρ2 ) (Equation (6)). This is done by adding a reference model in the objective function
(Equation (5)). ρre f 1 and ρre f 2 are generally chosen equal to 0.1 and 10 times the mean
measured apparent resistivity in the data set.
DOI =
log10 ρ1 − log10 ρ2
log10 ρre f 1 − log10 ρre f 2
(6)
A low value of the DOI approaching zero indicates that the discrepancy of the inverted
electrical resistivity between the two inversions is insignificant [52]. This is observed in wellresolved zones where the model is well constrained by the data. The inverted model is then
considered to be reliable. Conversely, high values of the DOI approaching unity indicate
that the inverted model is not well controlled by the data but rather by the regularization
procedure (regularization and reference model). Those zones should be interpreted with
caution. The DOI can be used to derive the maximum depth limit for which resistivity
values are driven by the data rather than through the constraint of the reference model [52].
A threshold value is often used to delimit the reliable part of the profile. DOI values
between 0.1 and 0.2 are often considered [52–54], but this threshold is rather subjective.
More advanced methodologies have been proposed to objectify the choice of the threshold [55]. In this paper, we follow the proposal of Caterina et al. [51] to use the gradient in
the DOI to define the threshold: when the DOI increases sharply with depth, it indicates
that the role of the reference model becomes predominant in the inversion [51].
4. Petrophysics and Salinity Threshold
In porous sediments, the bulk electrical resistivity mostly depends on the resistivity
of the pore fluid, the degree of saturation, and the presence of clay. The petrophysical
relationship is classically expressed using Archie’s law [56,57]:
σf
1
= σb =
ρb
F
(7)
where ρb is the bulk electrical resistivity; σb and σ f represent the bulk and fluid electrical
conductivity, respectively; and F refers to the formation factor, depending mostly on the
porosity φ and two empirical parameters: m, the cementation exponent, varying between 1
and 4, and a, an empirical constant close to 1:
F=
a
φm
(8)
If clayey sediments are present, surface conduction related to the electrical double
layer takes place and an additional term, the surface conductivity σs , which depends on
the surface mobility of the ions in the electrical double layer and the excess surface charges
Water 2021, 13, 1743
9 of 24
taking part in surface conduction per pore volume unit, is needed in the petrophysical
relationship [10]:
σf
σb =
+ σs
(9)
F
In the study area, the Holocene and Pleistocene aquifers are composed of a mixture of
clayey and sandy sediments and are affected by the presence of saltwater. Both clay and
salinity reduce the bulk electrical resistivity, and unequivocally interpreting low resistivity
values as saltwater intrusions might be misleading.
Nevertheless, it should be noted that according to Equation (9), the effect of clay
dominates the response when the fluid conductivity is low, i.e., in freshwater conditions. At
high salinity, clay sediments only relatively contribute to the conductivity response, so that
the lowest resistivity values can generally be linked to the presence of saltwater. Therefore,
we use the following considerations in the interpretation of the resistivity profiles.
1.
2.
3.
4.
Several profiles—for instance, 20 and 29 (Figure 4)—were collected in both freshwater
conditions and clay-rich sediments as deduced from nearby water samples (Figure 6)
and lithological descriptions of contiguous boreholes (e.g., Figure 4). On those profiles,
the lowest inverted resistivity is in the range from 10 Ohm.m to 15 Ohm.m (Figure 4),
which will be considered as the threshold for clayey sediments in freshwater conditions. Resistivity lower than the threshold will be considered as indicating the
presence of salt or brackish water (Figure 6).
The formation factor can be estimated using either colocated or laboratory measurements from Equation (7). Given the heterogeneity of the aquifer, the latter options
would require collecting many samples and carrying out extensive laboratory experiments and bear the risk of the scale representability [14]. As an alternative, we identify
the wells for which water samples were collected in close vicinity (but not colocated)
of the ERT profiles, assuming that the observed salinity remains representative of
the area. The well location was perpendicularly projected on the profile, and the
bulk inverted resistivity at the corresponding depth was extracted. To account for
the uncertainty of the projection, we extended the analysis to the neighborhood of
the projected location in a range corresponding to the distance of the borehole to
the profile.
A fractional correction was also used for the temperature compensation, assuming
that the increase of electrical conductivity is 0.02 ◦ C−1 , in order to compare the
conductivity at the reference temperature of 25 ◦ C [58–60]. The temperature of the
groundwater at the time of the survey was 30 ◦ C.
A linear regression (Figure 5) shows that an average formation factor from 2.2
(Equation (7)) to 2.8 (Equation (9)) explains the observed trend, consistent with values
typically observed in sandy sediments.
Water 2021, 13, 1743
10 of 24
Figure 4. Relationship between resistivity and the lithological composition of boreholes. The low resistivity features are
characterized by higher clay content in both saturated and fractured zones. Note that Qh, Qp are Holocene and Pleistocene
aquifers, respectively; K refers to Deo Ca Cretaceous granite complex. The same notation is used in the following figures.
Figure 5. Formation factor of unconsolidated sediment in saturated conditions. Blue points represent
extracted values from ERT measurement corresponding to the individual borehole, and red points
are mean values. Dotted lines are average linear trends adjusted to Archie’s law (Equation (7), blue)
and Waxman and Smits (Equation (9), red).
Water 2021, 13, 1743
11 of 24
Figure 6. Distribution of the water samples in the shallow boreholes collected in the study area and their TDS content.
Note that the samples were collected after the salinity map [34] was produced. The distribution of salinity in the upstream
part is very heterogeneous (green lines correspond to ERT profiles, black lines correspond to ERT profiles that failed, and
the two purple solid lines indicate the previously saline boundaries corresponding to TDS limit = 1.5 g/L and 1.0 g/L,
respectively) [34]. Topographical isolines (in m above sea level) are shown by blackish brown curves.
Based on those considerations, it is now possible to estimate a resistivity threshold
for the saline boundary. Based on the fluid electrical conductivity measured on the water
samples (Figure 6), we approximate the relationship between fluid electrical conductivity
(in S/m) and total dissolved solid content (TDS in mg/L) as linear, with σ f ≈ TDS
α [61].
The conversion factor α is commonly estimated between 5000 and 6500 [61]. The electrical
resistivity of water at total dissolved solids (TDS) of 3 g/L, interpreted as the saltwater transition, would correspond approximately to 1.6 Ohm.m in the former case and 2.2 Ohm.m
in the latter. Considering a formation factor between 2.2 and 2.8, the corresponding bulk
electrical resistivity should be comprised of between 3.5 Ohm.m and 6 Ohm.m. Note that
we used Archie’s law for the conversion, as the estimation of the surface conductivity from
a limited number of points is inaccurate and its role for high salinity is limited.
According to these considerations, we can qualitatively interpret a resistivity value
below 12 Ohm.m as characteristic of brackish (more than 1.5 g/L) to saline water. Similarly,
resistivity values over 15 to 20 Ohm.m are indicative of freshwater (see color scale of
Figure 4). To account for possible error in the estimation of the conversion factor and in the
formation factor, we set our limit for the saline boundary (3 g/L) at a conservative limit of
7 Ohm.m.
Note that the actual threshold does not significantly affect the interpretation. The
uncertainty related to our petrophysical analysis (in the range of 2 Ohm.m) is low compared
to the uncertainty resulting from the inversion. Inverted resistivity values commonly
overestimate the actual resistivity in saline conditions [6,18] but might enlarge the footprint
of the low resistivity zone because of smoothing occurring during the regularized inversion
process. Using conservative thresholds with a grey zone corresponding to resistivity in the
range of 7 to 20 Ohm.m allows us to define saline and freshwater zones in a robust way.
Water 2021, 13, 1743
12 of 24
In the next section, we use a unique color scale to display all the ERT results across
the study area (Figure 4). According to the above considerations and related uncertainty,
we chose a conservative threshold resistivity value of 7 Ohm.m to characterize the saltwater transition, with saline water occurring in blue areas with a resistivity lower than
7 Ohm.m. Freshwater is present in light yellow-brown areas with a resistivity higher than
approximately 20 Ohm.m. Brackish water is found in the medium color areas ranging
between bright yellow and bright light blue. Note that this range of value can also be
characteristic of the presence of clayey sediments in fresh conditions and is thus the most
difficult to interpret.
5. Results and Discussion
Considering the different topographical and sedimentary features, we divide our
study site into two regions, namely the upstream and downstream parts. The limit between
the two zones is chosen at a distance from the sea of about 7.5 to 8 km (Figure 1). This
distance corresponds to the upstream saltwater intrusion limit in the river estuary during
high tide [35]. The vertical axis on the figure always refers to the elevation expressed in m
above sea level. Note that according to our thresholds, the term saline water will be used for
resistivity below 7 Ohm.m. Brackish and freshwater correspond to higher resistivity values.
Except if specified otherwise, these terms do not correspond to the salinity measured in
water samples; they are simply used to ease the interpretation.
5.1. Upstream Zone
The upstream part is composed of eight profiles placed on cultivated fields with an
alternation of dragon fruit and rice plantations (Figure 1), lying on the complex mixture of
alluvial-marsh-marine alluvial deposits from the Holocene.
Figure 7 shows the inverted results of four profiles in the upstream part. The distribution of saltwater (ρ < 7 Ohm.m) is discontinuous and seems to be limited to localized
areas, as can be seen in profiles 19 and 22. The variation of resistivity values ranging
from 1 Ohm.m to 7 Ohm.m at shallow depths indicates that saline water dominates the
major part of profile 21. Brackish to fresher water lenses are commonly observed near the
surface, characterized by higher resistive features above 10 Ohm.m. In this area, saline and
brackish water seem to occupy mostly areas with low elevation along the Luy River. In
the areas with slightly higher elevation, the lateral variation of resistivity features ranges
from 10 Ohm.m to 20 Ohm.m, likely indicating that freshwater is dominantly present
as observed in profile 20 (Figure 4), where clay lenses are probably responsible for the
intermediate resistivity values. Profile 25 is adjacent to the river with the end of the profile
corresponding to the river bank. Here, the intermediate range of resistivity might not only
relate to clayey deposits but also point towards brackish water resulting from the direct
interaction between the river and the aquifer. However, according to recent investigations,
the part of the river directly affected by seawater intrusions during high tides would be
limited to the downstream part [28]. This observation would rule out recent infiltration
from the river to explain a higher salinity. Nevertheless, shallow water samples have
revealed strong salinity variations at a small scale as well as the presence of saltwater in
the vicinity of profiles 21 and 25 and further upstream (Figure 6). This is in accordance
with our ERT profiles that indicate the presence of saline-water lenses within an almost
fresh aquifer.
Water 2021, 13, 1743
13 of 24
Figure 7. Resistivity distribution of representative profiles and correlation with lithological description in the upstream part.
The variation in color indicates resistive characteristics. Blue corresponds to saline water. Freshwater is present in light
yellowish-brown areas and brackish water in intermediate colors. Qh, Qp: Holocene and Pleistocene aquifers; K: Deo Ca
Cretaceous granite complex.
A reasonable assumption for this distribution could be the remnants of a past transgressive systems tract (TST) having taken place during the Pleistocene and Holocene [42,62,63].
Although these saline-water lenses would be expected to freshen with time through natural
recharge, those processes could be slowed down in clay-rich layers in limited recharge
conditions as observed in this arid region. The role of heterogeneity could explain the
apparently discontinuous distribution of saline water in the aquifer of the upstream part.
However, this should be confirmed by further investigations.
The sharp vertical variation of resistivity at depth indicates the transition to the underlying rock basement. This transition occurs at different elevations varying significantly
from −10 to −30 m. As will be demonstrated in Section 5.3, this part of the profile is still
located in the resolved part of the tomogram. Inhomogeneous structures are observed
in the granite bedrock with resistivity values reaching from 500 Ohm.m to 1000 Ohm.m
(Figure 7c,d). The transition between the unconsolidated layer and the bedrock is smooth,
but this could be an effect of inversion, even if a blocky model constraint (L1-norm) was
Water 2021, 13, 1743
14 of 24
used to generate sharp contrasts during inversion. This could also correspond to a weathered zone in the upper part of the fresh bedrock. By comparison with the limited boreholes
reaching the bedrock in the area, it was estimated that the depth of the bedrock was best
corresponding with the 35-Ohm.m resistivity isoline.
As can be seen, the bedrock limit is not horizontal and displays lateral variations in the
bulk resistivity. By combining it with borehole logs in the area [34], several specific features
in the inversion results located in the basement are considered to correspond to small
blocks of younger basement related to past magmatic activity (Figure 7b,c) [36,39]. We
hypothesize that basaltic fluids erupted and penetrated through the older granitic bedrock
that could be deformed and cracked. These assumptions are coherent with the presence of
basaltic fragments observed in the nearby borehole LK15-BT and a basaltic bedrock layer
overlying the granite in borehole LK9S (Figures 1 and 7), located at approximately 2 km
and 700 m of the profiles, respectively. Lower resistivity values ranging from 30 Ohm.m to
100 Ohm.m in the bedrock might be related to the increase of secondary porosity through
fractures and the presence of secondary clay minerals originating from the weathering of
the granite bedrock. Both features were also observed in boreholes.
5.2. Downstream Zone
In the downstream part, the flat plains at low elevations on the left bank of the Luy
River are inclined toward the sea, with a slope ranging from 0 to 5◦ [31,43]. Seawater can
intrude along the river in high tide conditions which can directly impact the aquifers. Nine
ERT profiles were collected on the left bank, with five of them being in the direct vicinity
(<700 m) of the Luy River or its tributaries. Figure 8 shows extremely low resistivity in
all profiles, fluctuating between 0.5 and 5 Ohm.m, which points towards the presence
of high salinity values. Higher resistivity values varying from 7 Ohm.m to 35 Ohm.m
are locally found near the surface. Saline water thus mostly occupies both the Holocene
and Pleistocene unconsolidated aquifers. Freshwater is only present in lenses (Profile
27) or a thin layer (Profile 17) near the surface and is probably directly controlled by the
infiltration of rainwater. In comparison with Figure 1, where the previously estimated
saltwater boundary is depicted, the extent of saltwater migration into the aquifers seems to
expand to the northeast and the north on the left bank of the river. Moreover, extremely
high salinity is locally observed and indicated by a bulk resistivity lower than 1 Ohm.m
in profile 17 (Figure 8b), which is close to a tributary of the Luy River (Cau Nam River,
Figure 1). According to the petrophysical relationship, such a low value indicates a salinity
close to that of seawater.
Closer to the Luy River, the succession of three profiles (namely 11, 23, and 12; Figure 9)
also shows a complex distribution of salinity similar to that observed in the upstream part.
Profiles 11 and 12 bear saline to brackish water with resistivity ranging from 1 Ohm.m
to 10 Ohm.m. In contrast, profile 23 seems to have brackish to fresh water conditions as
demonstrated by higher resistivity values fluctuating between 10 Ohm.m to 30 Ohm.m.
This further illustrates that the fresh–saline conditions can vary rapidly within the aquifer.
The salinity measured in shallow boreholes also visibly depicts this complexity in the fresh
and saline water distribution (Figure 6).
Water 2021, 13, 1743
15 of 24
−
−
Figure 8. Resistivity distribution of representative profiles on the left bank of the upstream part. Qh, Qp: Holocene and
Pleistocene aquifers; K: Deo Ca Cretaceous granite complex.
Figure 9. Resistivity distribution of the three consecutive profiles on the left bank of the upstream part. The distance
between LK03-BT and the profiles is small enough (approximately 500 m) to compare the borehole log with ERT results. Qh,
Qp: Holocene and Pleistocene aquifers; K: Deo Ca Cretaceous granite complex.
Water 2021, 13, 1743
16 of 24
As for the upstream zone, the transition to the underlying granite basement is well imaged by an increase in resistivity at an elevation of around −20 to −30 m. Lateral variations
of electrical resistivity illustrate the heterogeneous structures of this granite bedrock, which
is likely related to fractures and alterations. In Figure 8c, the resistivity in the basement
reaches 500 Ohm.m, possibly indicating unaltered granite [64–66]. Inversely, quite low
resistivity ranging between 10 Ohm.m and 30 Ohm.m is observed in the upper part of
the bedrock. We assume that the electrical resistivity decrease of the bedrock indicates
the presence of secondary clay formed during weathering processes. The unexpected low
resistivity locally observed in the bedrock could point towards the presence of saline water
that would have infiltrated from the unconsolidated aquifers through the hydraulically
conductive fractures of the basement due to the density effect. However, these hypotheses
could only be confirmed with direct observations.
On the right bank of the river, the sandy dunes correspond to a high elevation zone
(50 m to 170 m above sea level, Figure 6). In the upper part of these deposits, high electrical resistivity values up to 5000 Ohm.m are imaged (Figure 10), indicating unsaturated
conditions. At a larger depth, lower resistivity features ranging between 50 Ohm.m to
350 Ohm.m show that the remnants of the Holocene aquifer and the upper part of the
Pleistocene aquifer are filled with freshwater. Deeper, the underlying Pleistocene aquifer
contains cleaner and finer sands, probably in a more compacted state [67]. It is characterized
by a sharp variation in resistivity ranging from 280 Ohm.m to 800 Ohm.m (Figure 10a). At
an elevation of around −10 m, the electrical resistivity decreases significantly to a minimum
value of 10 Ohm.m (Figure 10a,c). This elevation corresponds to the depth of the bedrock.
We hypothesize that the top of the basement is either composed of clay-rich minerals in a
strongly weathered condition or filled with saline water. However, the resolution in this
zone is not sufficient, as the bedrock is located deeper below the surface compared to other
investigated zones (see Section 5.3).
Figure 10. Resistivity distribution of representative profiles on the right bank of the upstream part. The consecutive drilled
boreholes (LK 19/26, LK20/31, and LK19/27) are close to profile 4, and LK16-7/10 coincides with profile 10. The variation
in resistivity of lithological layers is similar to the discrepancy in the state and structure of the sedimentary sequence and
bedrock. Qh, Qp: Holocene and Pleistocene aquifers; K: Cretaceous intrusive and extrusive rocks; N2 : Late Neogene Lien
Huong formation.
The dune area constitutes an important freshwater lens larger than 50-m thick, as
observed in profiles 4, 9, and 10 (Figure 10). As it constitutes a preferential recharge zone,
Water 2021, 13, 1743
17 of 24
a freshwater lens has developed in the dunes and creates a natural hydraulic gradient
towards the river and the sea. This possibly explains why less occurrence of saline water is
observed on the right bank. It also seems to prevent the direct intrusion of saltwater from
the sea, as illustrated in profile 4. The boundary of the saline–fresh interface is thus close to
the sea on the right bank.
5.3. Depth of Investigation
Figure 11 shows the DOI for one representative profile in each of the identified zones,
as other profiles within the same zone show similar patterns. The DOI maps show the
expected patterns, with the DOI approaching zero close to the surface down to about 35
to 50 m below the surface. At this depth, a zone with intermediate DOI values (0.2 to 0.4)
visibly indicates that the two inversion models have slight discrepancies. For profiles in the
upstream zone and on the left bank of the downstream zone (profiles 11 and 22), these zones
correspond to the more resistive conditions of the bedrock. Typically, higher resistivity
values are less well resolved, as current lines flow preferentially in conductive media.
Therefore, the bedrock is well constrained as a resistive body, but its actual resistivity is
uncertain. This is confirmed by the fact that the DOI increases sharply only at a larger depth
between 75 and 85 m, which we consider as our threshold for determining the maximum
depth of investigation. The same is observed in the right bank (profile 10), but the higher
DOI values correspond to the depth of the transition to the bedrock. Below this limit, the
results are not well constrained by the data. Therefore, both the unconsolidated aquifer
and the upper part of the bedrock are located in well-resolved areas and can be interpreted,
except on the right bank where the bedrock depth corresponds to a higher DOI, indicating
that interpretation of the bedrock should be made with caution, although the decrease in
resistivity itself seems consistent in all representative profiles.
5.4. Advantages and Limitations of the Methodology
In this paper, we have tackled the challenging task of interpreting ERT data in terms
of salinity with sparse colocated data. As with any attempts to quantitatively interpret ERT
results, our methodology suffers from limitations that are discussed below.
As can be seen from Figure 6, the estimation of the formation factor is limited by the
amount of available data. Because of the absence of colocated data, we had to consider
the uncertainty related to the distance between the ERT and wells, so that only five wells
were considered sufficiently close, and a large spread of values is observed. Some points
also deviate from the trend, likely indicating that a single petrophysical relationship might
not be sufficient to characterize the study area, as expected from the heterogeneity of the
sediments. On the one hand, it is clear that the methodology could benefit from a more
robust estimation of the formation factor. On the other hand, the few available points are
not only enough to validate the expected trend but also derive a formation factor coherent
with values found in the literature.
It is well known that the bulk resistivity values obtained after inversion can deviate
strongly from the real bulk resistivity values. If advanced methodologies incorporating
prior information can help to reduce the discrepancy, these often require additional data.
It is clear that low resistivity values are likely overestimated, high resistivity values are
underestimated, and the transition zone might be an effect of smoothing rather than the
truth. This can obviously alter the interpretation of ERT images.
Water 2021, 13, 1743
18 of 24
Figure 11. Three representative inverted images and corresponding DOI values. Threshold values to define the maximum
penetration depth are chosen from 0.5 to 0.6. In each profile, the bottom image is a DOI map, while the top image is the
ERT inversions filtered using the DOI. The well-resolved regions lie within low DOI index (dark blue); inversely, high DOI
values indicate zone poorly constrained by the data (orange, red, and pink).
Nevertheless, our methodology allows addressing both limitations. Using a conservative limit for the saline and fresh boundaries, we limit ourselves to a semiquantitative
interpretation (saline, brackish, fresh), but we avoid introducing errors in the interpretation.
Indeed, the choice of the thresholds was conservative to avoid the influence of errors made
in the petrophysical relationship and left a “grey zone” between saline and fresh classes,
accounting for the likely presence of clay lenses in the study area, which is difficult to
discriminate from brackish water. The interpretation in terms of boundaries is still relevant,
as the difference between 5 and 10 g/L makes little difference for farmers: in any case, the
water is above the useful limit. It also allows integrating geophysical data in conceptual
hydrogeological models or even more broadly in the calibration of numerical models.
In addition, the systematic use of the DOI map allows limiting the interpretation
to reliable zones within the inverted images. In our case, it allows results down to a
depth reaching 80 m to be interpreted, which is sufficient for the characterization of
unconsolidated sediments.
The complete ERT inversion models produce a quite clear overview of the link between
the bulk resistivity values and both the stratigraphic units, especially of the bedrock, and
the salinity.
6. The New Fresh–Saline Groundwater Map
The ERT results are used to update the previous saline boundary map from NAWAPI [34].
In Figure 12, we summarized every ERT profile as a two-layer log displaying the average
Water 2021, 13, 1743
19 of 24
resistivity value of the aquifer and the bedrock as well as the depth of the bedrock. The blue
scale indicates the salinity in the unconsolidated sediments, while the red corresponds to
freshwater conditions. The boundary (yellow solid line) corresponds to our interpretation
based on ERT, corresponding to an estimated TDS value of 3 g/L in the unconsolidated
aquifers. Dotted segments indicate a lack of data to conclude with certainty. Local variations within the zones (presence of salt/fresh water lenses) are possible as described in
Section 5. Note that the original map used a threshold of 1.5 g/L; the extension of 3 g/L
based on shallow water samples would have been even more restricted.
Figure 12. Three-dimensional view of the survey area with a color scale cropped to the saline (7 Ohm.m) and fresh (20
Ohm.m) thresholds. Note that the elevation axis appears vertically on the figure, and the scale is exaggerated. The updated
saline boundary in the Luy River catchment is shown in yellow, while the old boundaries are shown in purple.
Generally, the boundary of saline water expands toward the north and northeast up
to at least 1.5 km compared to the initially mapped saltwater extension. This observation
is mostly based on the profiles collected on the downstream part of the left bank of the
river. The extension could even be larger as neither water samples nor ERT profiles were
collected further in the north. On the right bank, our ERT profiles show that the dune area
and its freshwater lens constitute a protection against SI, as the hydraulic head is larger and
groundwater flow is directed towards the river and the sea. Therefore, the zone of saltwater
intrusion is narrower in the southeast than previously estimated. If this zone could be seen
as a freshwater reserve, SI could eventually occur if this lens was uncontrollably exploited.
Water 2021, 13, 1743
20 of 24
The discrepancy with the previous map is probably related to the fact that ERT
investigates the whole aquifer thickness, while the old map was based on shallow boreholes.
Our ERT investigations show that a freshwater lens can exist on top of the saline water. This
observation is also corroborated by the experience of farmers who exploit shallow wells
and by our latest investigation during the period from 2020 to the first half of 2021, during
which an increase in salinity was observed in the dry season because of overexploitation.
An alternative solution would be to accumulate freshwater during the rainy season for
consumption in the dry season.
In the upstream zone and part of the downstream zone, the distribution of saline
water is extremely heterogeneous and complex as shown by both the high variability of
salinity from water samples and ERT profiles at short distances. The presence of SI cannot
be explained only by interaction with the river, as previously thought, but might also
result from the interplay between past transgressions with connate saltwater filling the
sediment pores. The presence of clay layers and lenses acting as flow barriers and trapping
connate water, complicated recharge patterns influenced by the presence of clay near the
surface, and agricultural practices including water extraction and irrigation could prevent
freshening processes. This is further confirmed by the presence of saline water towards
the western part of the river catchment (Figure 7), but no ERT profiles were collected in
that zone.
The origin of saltwater is thus probably a complex combination of different processes.
Recent intrusions related to the flow regime of the Luy River likely play a notable role
in the downstream part, including along tributaries, where high tides in areas with low
elevations can result in direct infiltration from the river into the aquifer. The uncontrolled
exploitation of the aquifer for irrigation and industrial production can further accelerate
the migration of saline water.
However, ERT data shed new light on the distribution of saltwater in the area. The
presence of freshwater lenses of a limited extent on top of the more saline and deeper
aquifer could indicate that the development of a freshwater aquifer, as expected from
natural recharge, is slow and likely perturbed by human activities, namely water extraction
and irrigation practices, aquaculture, and salt production leading to a depletion of these
limited freshwater resources. Moreover, the complex distribution of saltwater, including
in the upstream part of the catchment, far from the sea and the river, seems to indicate
the presence of connate saltwater buried in the sediments that has not been freshened yet.
These hypotheses cannot be resolved from ERT data only but require additional analysis,
such as hydrochemistry and flow and transport modeling [68].
Moreover, saltwater intrusions could also reach deeper zones of the bedrock under the
sea level. If the presence of saltwater in the basement in the downstream part of the river on
both banks, including under the sand dunes, was confirmed, this would likely indicate that
the origin of almost all the saltwater is related to past transgressions. Indeed, if freshening
processes are expected, for example, under the dune areas, these can be impeded by low
hydraulic conductivity zones slowing down groundwater flows, as could occur in fresh
bedrock, weathered zones containing a significant content of secondary clay, or in clay
sedimentary deposits. But this should be confirmed by direct borehole information, as the
resolution of ERT at this depth is not sufficient to make an unequivocal conclusion.
7. Conclusions
In this paper, we developed a new methodology to interpret ERT data in terms of salinity in a heterogeneous aquifer in the absence of colocated data. We defined a conservative
threshold of resistivity to interpret semiquantitatively saline and fresh boundaries.
We applied the methodology to the Luy River water catchment, Binh Thuan province,
Vietnam. The estuaries of the rivers in Binh Thuan are assumed to be the key driver for
saltwater intrusion along the coastal areas during high tides [28]. However, our investigation revealed that the zone affected by SI is much larger than previously expected
from shallow water samples. In particular, the bottom part of the unconsolidated aquifer
Water 2021, 13, 1743
21 of 24
systematically bears saltwater on the left bank of the downstream part of the Luy River
system. The high elevation dune complex constitutes a freshwater reserve that partly
protects the right bank of the river from SI. ERT also reveals a complex distribution of SI
in the upstream part of the river, falsifying the hypothesis that the river is the only driver
for SI in the area. Other processes, including paleohydrological conditions, geological
heterogeneity, and agricultural practices, likely play a preponderant role. We plan to build
on this preliminary study to develop a numerical model of the aquifer, test the different
hypotheses, and propose a sustainable management of the water resources.
As demonstrated in our study, ERT is a suitable tool to collect semiquantitative data
on saline water intrusion, even in the absence of extensive colocated data for validation.
It is cost effective and easy to apply both in rural and semiurban areas and should be
generalized for the characterization of seawater intrusion in complex environments.
Author Contributions: D.C.-T. conceptualized the survey plan, processed the ERT data, applied the
methodology, and was responsible for writing the paper. L.P.D. and R.T. contributed to the fieldwork
and data processing. M.P. contributed to data processing. H.H.H. participated in the fieldwork and
the conceptualization of the survey plan. F.N. and T.H. conceptualized the survey plan and the
methodology and supervised the study. All authors have read and agreed to the published version
of the manuscript.
Funding: This research is funded by VLIR-UOS and the Belgian development cooperation through
the grant VN2019TEA494A103, supporting the fieldwork and the PhD scholarships of the first two
authors.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: The data presented in this study are available on request from the corresponding author. The data are not publicly available yet due to the terms of the research agreement.
Acknowledgments: We deeply thank the VIGMR staff for their support in the field and NAWAPI for
data sharing. We would like to thank David Caterina from Liege University, who assisted us in the
analysis of reciprocal errors.
Conflicts of Interest: The authors declare no conflict of interest.
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
Werner, A.D.; Bakker, M.; Post, V.E.A.; Vandenbohede, A.; Lu, C.; Ataie-Ashtiani, B. Seawater intrusion processes, investigation
and management: Recent advances and future challenges. Adv. Water Resour. 2013, 51, 3–26. [CrossRef]
Saeed, M.M.; Bruen, M.; Asghar, M.N. A review of modeling approaches to simulate saline-upconing under skimming wells.
Nord. Hydrol. 2002, 33, 165–188. [CrossRef]
Stigter, T.Y.; Van Ooijen, S.P.J.; Post, V.E.A.; Appelo, C.A.J.; Carvalho Dill, A.M.M. A hydrogeological and hydrochemical
explanation of the groundwater composition under irrigated land in a Mediterranean environment, Algarve, Portugal. J. Hydrol.
1998, 208, 262–279. [CrossRef]
Turner, I.L.; Acworth, R.I. Field measurements of beach face salinity structure using cross-borehole resistivity imaging. J. Coast.
Res. 2004, 20, 753–760. [CrossRef]
Gossel, W.; Sefelnasr, A.; Wycisk, P. Modelling of paleo-saltwater intrusion in the northern part of the Nubian Aquifer System,
Northeast Africa. Hydrogeol. J. 2010, 18, 1447–1463. [CrossRef]
González-Quirós, A.; Comte, J.C. Relative importance of conceptual and computational errors when delineating saltwater
intrusion from resistivity inverse models in heterogeneous coastal aquifers. Adv. Water Resour. 2020, 144, 103695. [CrossRef]
Binley, A.; Hubbard, S.S.; Huisman, J.A.; Revil, A.; Robinson, D.A.; Singha, K.; Slater, L.D. The emergence of hydrogeophysics for
improved understanding of subsurface processes over multiple scales. Water Resour. Res. 2015, 51, 3837–3866. [CrossRef]
Paepen, M.; Hanssens, D.; De Smedt, P.; Walraevens, K.; Hermans, T. Combining resistivity and frequency domain electromagnetic
methods to investigate submarine groundwater discharge (SGD) in the littoral zone. Hydrol. Earth Syst. Sci. 2020. [CrossRef]
Goebel, M.; Knight, R.; Halkjær, M. Mapping saltwater intrusion with an airborne electromagnetic method in the offshore coastal
environment, Monterey Bay, California. J. Hydrol. Reg. Stud. 2019, 23, 100602. [CrossRef]
Viezzoli, A.; Munday, T.; Cooper, Y.L. Airborne electromagnetics for groundwater salinity mapping: Case studies of coastal and
inland salinization from around the world. B Geofis. Teor. Appl. 2012, 53, 581–600. [CrossRef]
Water 2021, 13, 1743
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
22 of 24
Delsman, J.; Van Baaren, E.; Vermaas, T.; Karaoulis, M.; Bootsma, H.; De Louw, P.; Pauw, P.; Oude Essink, G.; Dabekaussen, W.;
Van Camp, M.; et al. TOPSOIL Airborne EM Kartering Van Zoet en Zout Grondwater in Vlaanderen (FRESHEM Vlaanderen).
Vlaamse Milieumaatschappij & Deltares. 2019. Available online: https://www.vmm.be/publicaties/airborne-em-kartering-vanzoet-en-zout-groundwater-in-vlaanderen (accessed on 1 October 2019).
Loke, M.H.; Chambers, J.E.; Rucker, D.F.; Kuras, O.; Wilkinson, P.B. Recent developments in the direct-current geoelectrical
imaging method. J. Appl. Geophys. 2013, 95, 135–156. [CrossRef]
De Franco, E.; Biella, G.; Tosi, L.; Teatini, P.; Lozej, A.; Chiozzotto, B.; Giada, M.; Rizzetto, F.; Claude, C.; Mayer, A.; et al.
Monitoring the saltwater intrusion by time lapse electrical resistivity tomography: The Chioggia test site (Venice Lagoon, Italy). J.
Appl. Geophys. 2009, 69, 117–130. [CrossRef]
Nguyen, F.; Kemna, A.; Antonsson, A.; Engesgaard, P.; Kuras, O.; Ogilvy, R.; Gisbert, J.; Jorreto, S.; Pulido-Bosch, A. Characterization of seawater intrusion using 2D electrical imaging. Near Surf. Geophys. 2009, 7, 377–390. [CrossRef]
Kuranchie, F.A.; Shukla, S.; Habibi, D.; Zhao, X.; Kazi, M. Studies on Electrical Resistivity of Perth Sand. Int. J. Geotech. Eng. 2014,
8, 449–457. [CrossRef]
National Center for Water Resources Planning and Investigation (NAWAPI). Project: Mapping of Groundwater Resources at Scale
1:200.000 for the Whole Provinces in Vietnam. 2018. Available online: http://nawapi.gov.vn/index.php?option=com_content&
view=article&id=4712%3A-an-lp-bn-a-cht-thy-vn-t-l-150000-cac-tnh-ninh-thun-va-binh-thun-nhng-kt-qu-t-c&catid=62%3
Anhim-v-chuyen-mon-hoan-thanh&Itemid=134&lang=vi (accessed on 15 April 2018). (In Vietnamese)
Attwa, M.; Günther, T.; Grinat, M.; Binot, F. Evaluation of DC, FDEM, and IP resistivity methods for imaging perched saltwater
and a shallow channel within coastal tidal flat sediments. J. Appl. Geophys. 2011, 75, 656–670. [CrossRef]
Ussher, G. Understanding the resistivities observed in geothermal systems. In Proceedings of the World Geothermal Congress,
Kyushu-Tohoku, Japan, 28 May–10 June 2000.
Delsman, J.; Van Baaren, E.S.; Siemon, B.; Dabekaussen, W.; Karaoulis, M.C.; Pauw, P.; Vermaas, T.; Bootsma, H.; De Louw, P.G.B.;
Gunnink, J.L.; et al. Large-scale, probabilistic salinity mapping using airborne electromagnetics for groundwater management in
Zeeland, the Netherlands. Environ. Res. Lett. 2018, 13, 084011. [CrossRef]
Revil, A.; Coperey, A.; Shao, Z.; Florsch, N.; Fabricius, I.L.; Deng, Y.; Delsman, J.R.; Pauw, P.S.; Karaoulis, M.C.; De Louw, P.G.B.;
et al. Complex conductivity of soils. Water Resour. Res. 2017, 53, 7121–7147. [CrossRef]
Hermans, T.; Vandenbohede, A.; Lebbe, L.; Martin, R.; Kemna, A.; Beaujean, J.; Nguyen, F. Imaging artificial salt water infiltration
using electrical resistivity tomography constrained by geostatistical data. J. Hydrol. 2012, 438–439, 168–180. [CrossRef]
Hermans, T.; Irving, J. Facies discrimination with electrical resistivity tomography using a probabilistic methodology: Effect of
sensitivity and regularisation. Near Surf. Geophys. 2017, 5, 13–25. [CrossRef]
Day-Lewis, F.D.; Singha, K.; Binley, A. Applying petrophysical models to radar travel time and electrical resistivity tomograms:
Resolution-dependent limitations. J. Geophys. Res. 2005, 110, B08206. [CrossRef]
Goebel, M.; Pidlisecky, A.; Knight, R. Resistivity imaging reveals complex pattern of saltwater intrusion along Monterey coast. J.
Hydrol. 2017, 551, 746–755. [CrossRef]
Palacios, A.; Ledo, J.J.; Linde, N.; Luquot, L.; Bellmunt, F.; Folch, A.; Marcuello, A.; Pilar Queralt, P.; Pezard, P.A.; Martínez,
L.; et al. Time- lapse cross-hole electrical resistivity tomography (CHERT) for monitoring seawater intrusion dynamics in a
Mediterranean aquifer. Hydrol. Earth Syst. Sci. 2020, 24, 2121–2139. [CrossRef]
Kerrou, J.; Renard, P. A numerical analysis of dimensionality and heterogeneity effects on advective dispersive seawater intrusion
processes. Hydrogeol. J. 2010, 18, 55–72. [CrossRef]
Yu, X.; Michael, H.A. Mechanisms, configuration typology, and vulnerability of pumping-induced seawater intrusion in
heterogeneous aquifers. Adv. Water Resour. 2019, 128, 117–128. [CrossRef]
Hoa, T.T.; Campbell, J.B.; Wynne, R.H.; Yang, S.; Son, P.V. Drought and Human Impacts on Land Use and Land Cover Change in
a Vietnamese Coastal Area. Remote Sens. 2019, 11, 333. [CrossRef]
Hien, L.T.T.; Gobin, A.; Huong, P.T.T. Spatial indicators for desertification in south-east Vietnam. Nat. Hazards Earth Syst. Sci.
2019, 19, 2325–2337. [CrossRef]
Department of Natural Resources and Environment of Binh Thuan (DONRE BT). Planning for Allocation and Protection of
Underground Water Resources in Coastal Sandy Areas of Binh Thuan Province, DONRE, Binh Thuan. 2018. Available online: https://thuvienphapluat.vn/van-ban/tai-nguyen-moi-truong/nghi-quyet-47-nq-hdnd-2018-phan-bo-bao-ve-tai-nguyennuoc-duoi-dat-cat-ven-bien-binh-thuan-380156.aspx?v=d (accessed on 1 December 2018).
Vietnam Institute of Meteorology Hydrology and Climate change (IMHEN). Project: Determining Adaptation and Prevention
under Climate Change for Binh Thuan Province. Binh Thuan 2012. (In Vietnamese)
Department of Natural Resources and Environment of Binh Thuan (DONRE BT). Mission of River Basin Planning in Binh Thuan
Province, DONRE. Binh Thuan 2013. (In Vietnamese)
Cao, D.G.; Thang, T.T.; Cuong, T.H.; Hiep, N.Q.; Khanh, N.Q.; Hung, D.M.; Hong, N.B.; Cuong, B.M.; Duong, P.T. Assessment
of the Current Situation of Salinity and Contamination in the Groundwater and the Capacity of Domestic Water Supply in the
Central Coastal Areas from Binh Dinh to Ba Ria-Vung Tau.VIGMR, Hanoi 2017, 128p. Available online: http://www.vigmr.vn/
2017/10/10/de-tai-khoa-hoc (accessed on 31 December 2017). (In Vietnamese)
Water 2021, 13, 1743
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.
60.
61.
23 of 24
National Center for Water Resources Planning and Investigation (NAWAPI). Project: Hydrogeological Mapping at Scale 1:50000
in Ninh Thuan and Binh Thuan Provinces, NAWAPI, Hanoi 2015. Available online: http://nawapi.gov.vn/index.php?option=
com_content&view=article&id=4712%3A-an-lp-bn-a-cht-thy-vn-t-l-150000-cac-tnh-ninh-thun-va-binh-thun-nhng-kt-qu-tc&catid=62%3Anhim-v-chuyen-mon-hoan-thanh&Itemid=134&lang=vi (accessed on 1 October 2015). (In Vietnamese)
Quang, D.N.; Ngan, V.H.; Tri, M.C.; Cong, M.V. On the control of saltwater intrusion: A case study for Binh Thuan province,
Vietnam. In Proceedings of the 10th International Conference on Asian and Pacific Coasts, Hanoi, Vietnam, 25–28 September 2019.
Tran, V.T.; Vu, K. Geology and Earth Resources of Vietnam; General Department of Geology and Minerals of Vietnam, Hanoi
Publishing House for Science and Technology: Hanoi, Vietnam, 2011.
Nguyen, D.T. Geological Maps and Explanatory Notes of B’Lao, Da Lat—Cam Ranh, Gia Ray—Ba Ria and Phan Thiet, Scale 1:
200,000. Cent. Inf. Arch. J. Geol. 1999. C-48-VI, C-48-XII, XVIII, C-49-I, II and C-49-VII. Available online: http://www.idm.gov.vn/
1P3NPIT/vi-VN/Ban-Do-Dia-Chat.aspx (accessed on 2 January 2000). (In Vietnamese)
Nguyen-Thuy, D.; Ta, P.H.; Nguyen-Van, H.; Van Dinh, H.; Van Dang, B.; Dang, N.H.; Do, H.T.T.; Nguyen, A.T.K.; Tran, T.D.; Van
Bui, V.; et al. Evaluation of Geological Heritage of Geosites for a Potential Geopark in Binh Thuan–Ninh Thuan Coastal Zone,
Vietnam. Geoheritage 2018, 11, 689–702. [CrossRef]
Tran, D.L.; Nguyen, X.B. Geology of Vietnam, the second chapter: Magmatic part. Cent. Inf. Arch. J. Geol. 1998, 251–258.
Available online: http://idm.gov.vn/nguon_luc/Xuat_ban/Anpham/magma2/HTM/t251.htm (accessed on 31 December 1998).
(In Vietnamese)
Tran, D.L.; Nguyen, X.B. Geology of Vietnam, the second chapter: Magmatic part. Cent. Inf. Arch. J. Geol. 1998, 241–250.
Available online: http://idm.gov.vn/nguon_luc/Xuat_ban/Anpham/magma2/HTM/t241a.htm (accessed on 31 December
1998). (In Vietnamese)
Tran, N.; Nguyen, D.D.; Dinh, V.T.; Vu, V.V.; Ma, K.K.; Trinh, N.T. Environment and forming mechanism of Phan Thiet Red sands.
Vietnam J. Geol. 1998, 245, 10–20. (In Vietnamese)
Vietnam-Netherlands Center for Water and Environment (VINWATER) and the Sub-institute of Hydrometeorology and Climate
Change (SIHYMECC). Project: Integrated Water Resources Management and Urban Development Related to Climate Change in
Binh Thuan Province, DONRE, Binh Thuan, 2017.
Hoang, P. Geological and Mineral Map of Phan Thiet, scale 1:50.000. Cent. Inf. Arch. J. Geol. 1997, C-49-37-A. Available online:
http://idm.gov.vn/1P1NPIT/vi-VN/Ban-Do-Dia-Chat.aspx (accessed on 1 December 1998). (In Vietnamese).
Dahlin, T.; Zhou, B. Multiple-gradient array measurements for multichannel 2D resistivity imaging. Near Surf. Geophys. 2006, 4,
113–123. [CrossRef]
Dahlin, T.; Zhou, B. A numerical comparison of 2D resistivity imaging with 10 electrode arrays. Geophys. Prospect. 2004, 52,
379–398. [CrossRef]
Overmeeren, R.A.V.; Ritsema, I.L. Continuous verticalelectrical sounding. First Break 1988, 6, 313–324. [CrossRef]
Slater, L.; Binley, A.M.; Daily, W.; Johnson, R. Cross-hole electrical imaging of a controlled saline tracer injection. J. Appl. Geophys.
2000, 44, 85–102. [CrossRef]
LaBrecque, D.J.; Miletto, M.; Daily, W.; Ramirez, A.; Owen, E. The effects of noise on Occam’s inversion of resistivity tomography
data. Geophysics 1996, 61, 538–548. [CrossRef]
Tikhonov, A.N.; Arsenin, V.A. Solution of Ill-Posed Problems; Winston & Sons: Washington, DC, USA, 1977.
Loke, M.H. Tutorial: 2-D and 3-D Electrical Imaging Surveys. 2011. Available online: http://geoelectrical.com (accessed on 1
September 2011).
Caterina, D.; Beaujean, J.; Robert, T.; Nguyen, F. A comparison study of different image appraisal tools for electrical resistivity
tomography. Near Surf. Geophys. 2013, 11, 639–657. [CrossRef]
Oldenburg, D.W.; Li, Y. Estimating depth of investigation in DC resistivity and IP surveys. Geophysics 1999, 64, 403–416. [CrossRef]
Parsekian, A.D.; Claes, N.; Singha, K.; Minsley, B.J.; Carr, B.; Voytek, E.; Flinchum, B. Comparing measurement response and
inverted results of electrical resistivity tomography instruments. J. Environ. Eng. Geophys. 2017, 22, 249–266. [CrossRef]
Marescot, L.; Loke, M.H.; Chapellier, D.; Delaloye, R.; Lambiel, C.; Reynard, E. Assessing reliability of 2D resistivity imaging in
mountain permafrost studies using the depth of investigation index method. Near Surf. Geophys. 2003, 1, 57–67. [CrossRef]
Deceuster, J.; Tanguy, R.; Nguyen, F.; Kaufmann, O. A modified DOI-based method to statistically estimate the depth of
investigation of dc resistivity surveys. J. Appl. Geophys. 2014, 172, 172–185. [CrossRef]
Archie, G.E. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 1942, 146, 54–62.
[CrossRef]
Glover, P.J.W. A modify Achie’Law for n phases. Geophysics 2010, 75, 247–265. [CrossRef]
Hermans, T.; Nguyen, F.; Robert, T.; Revil, A. Geophysical methods for monitoring temperature changes in shallow low enthalpy
geothermal systems. Energies 2014, 7, 5083–5118. [CrossRef]
Hermans, T.; Wildemeersch, S.; Jamin, P.; Orban, P.; Brouyère, S.; Dassargues, A.; Nguyen, F. Quantitative temperature monitoring
of a heat tracing experiment using cross-borehole ERT. Geothermics 2015, 53, 14–26. [CrossRef]
Hayashi, M. Temperature-electrical conductivity relation of water for environmental monitoring and geophysical data inversion.
Environ. Monit. Assess. 2004, 96, 119–128. [CrossRef]
Keller, G.V.; Frischknecht, F.C. Electrical Methods in Geophysical Prospecting; Pergamon Press: Oxford, UK, 1966.
Water 2021, 13, 1743
62.
63.
64.
65.
66.
67.
68.
24 of 24
Hoang, N.T.D.; Do, Q.T. Sea level changes in quaternary period revised from effectiveness of active tectonic subsidence in Quang
Nam coastal plain. Hue Univ. J. Sci. Earth Sci. Environ. 2018, 127, 29–36. [CrossRef]
Nguyen, V.T. Characteristics of Sedimentary Facies and Sequence Stratigraphy of the Quaternary Sand Formations in the Coastal
Areas of Ninh Thuan-Binh Thuan Provinces. Ph.D. Thesis, VNU University of Sciences, Vietnam National University, Hanoi,
Vietnam, 2019.
Ebong, E.D.; Abong, A.A.; Ulem, E.B.; Ebong, L.A. Geoelectrical Resistivity and Geological Characterization of Hydrostructures
for Groundwater Resource Appraisal in the Obudu Plateau, Southeastern Nigeria. Nat. Resour. Res. 2021, 30, 2103–2117.
[CrossRef]
Kumari, A.; Kumar, D.; Warwade, P. Application of Multi-Criteria Decision Making (MCDM) and Electrical Resistivity Tomography (ERT) Techniques for Identification of Groundwater Recharge zone(s) in Granitic Hard Rock Aquifer. J. Earth Syst. Sci. 2021,
130. [CrossRef]
Sreeparvathy, V.; Kambhammettu, B.V.N.P.; Peddinti, S.R.; Sarada, P.S.L. Application of ERT, Saline Tracer and Numerical Studies
to Delineate Preferential Paths in Fractured Granites. Groundwater 2019, 57, 126–139. [CrossRef] [PubMed]
Thao, T.V.; Cuong, P.H.; Chi, P.D.; Duong, P.S.; Huy, D.Q.; Hanh, T.V.; Lien, D.H.; Manh, N.T.; Nhan, T.N.; Phuong, T.Q.; et al.
Project: Investigation and Assessment of Titanium-Zircon Mineral Potential in the Reddish Sandy Layer in Provinces: Ninh
Thuan, Binh Thuan, and Northern Ba Ria—Vung Tau. Cent. Inf. Arch. J. Geol. 2008. Available online: http://idm.gov.vn/6P1
NPIT/vi-VN/Bao-Cao-Dia-Chat.aspx (accessed on 31 December 2008).
Sarker, G.; Berrens, R.; Von Arx, J.; Pelczar, P.; Reik, W.; Wolfrum, C.; Peleg-Raibstein, D. Transgenerational transmission of
hedonic behaviors and metabolic phenotypes induced by maternal overnutrition. Transl. Psychiat. 2018, 8, 195–211. [CrossRef]
[PubMed]