J. Nig. Soc. Phys. Sci. 5 (2023) 1512

Journal of the
Nigerian Society

of Physical
Sciences

Simulation of the Movement of Groundwater in an Aquifer

Suha Ibrahim Salih AL-Ali, Nihad Jalal Kadhem AL-Awsi

Computer Science and Mathematics Department, College of Computer Science and Mathematics, Tikrit University, Tikrit, Iraq

Abstract

This study investigates the impact of extracting fresh water from areas where salt water and fresh water meet in tropical regions. Traditionally,
fresh water is expected to be found above salt water in the ocean or underground. To carry out the investigation, Green’s Function method is used,
and a numerical chart is presented that includes an equation derived from Green’s II matching. The study computes the shape of the interface
during water withdrawal and flows through the basins and sources of the line. In addition, this study obtains an analytical solution to the linear
problem for the non-withdrawal scenario. Finally, the study identifies the maximum rate of water withdrawal before the initial breakthrough of
salt water for different density ratios.

DOI:10.46481/jnsps.2023.1512

Keywords: Free Surface, Tropical Island, Porous Medium, Freshwater Lens, Island Tub, Green’s Function.

Article History :
Received: 21 April 2023
Received in revised form: 19 May 2023
Accepted for publication: 22 May 2023
Published: 22 June 2023

c© 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license

(https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

Communicated by: B. J. Falaye

1. Introduction

The provision of drinking water in tropical islands often re-
lies on a freshwater lens trapped in the soil beneath the surface.
However, this supply is threatened by various factors, such as
the interface between freshwater and saltwater layers, surface
channels, and pumping practices. To ensure continued access
to freshwater, it is crucial to conserve the freshwater lens and
use efficient extraction practices to prevent saltwater intrusion
and depletion of the freshwater supply [1]. Figure (1) provides
a visual representation of this process. According to the United
Nations, around 2.2 billion people lack access to safe drinking
water, and this number is expected to increase due to population
growth and climate change [2]. Small island developing states

Email address: suhaibrahim3@tu.edu.iq (Suha Ibrahim Salih AL-Ali
)

(SIDS) are particularly vulnerable to water scarcity because of
their limited resources, geographic isolation, and susceptibility
to natural disasters [3]. Climate change exacerbates the water
scarcity crisis in SIDS by increasing the frequency and inten-
sity of extreme weather events, causing saltwater intrusion, and
reducing precipitation [4]. For example, the Pacific Island na-
tion of Tuvalu faces severe water shortages due to droughts and
over-extraction from the freshwater lens, which has led to soil
and groundwater salinization [5]. Some island communities
are responding to these difficulties by implementing creative
solutions such as rainwater collection, desalination, and water
conservation techniques. These solutions, however, necessitate
money and technical skill, both of which may be limited in
SIDS [6, 7]. Furthermore, policies that promote water secu-
rity and sustainability in these vulnerable areas are required [8,
9]. The interface between low-density top liquids and higher-
density liquids in a two-layer water object changes shape when

1



AL-Ali & AL-Awsi / J. Nig. Soc. Phys. Sci. 5 (2023) 1512 2

Figure 1. A diagram indicate the focus of research-pulling freshwater
over the underground seawater on an island

liquid is preferentially removed from one layer, according to
research [10]. One of the primary goals of this research is to
validate the deformed interface shape.

The goal of this study is to look into the variables that could
cause the interface to advance towards the extraction port, re-
sulting in liquid extraction from both levels and mixing. The re-
search focuses on a single island with defined top left and right
bounds that contains a freshwater lens atop a saltwater layer
[11]. To address this issue, academics employ the Darcy Act
as a foundation, which contains boundary conditions on both
sides and the island’s bottom. This study is in the field of hy-
drology, and it focuses on saltwater intrusion in coastal aquifers,
which can cause a variety of environmental and economic prob-
lems [12]. The authors assume the presence of a freshwater lens
above a saltwater layer, which is a common condition in real-
world scenarios like Aruba. The findings of the study can be
used to develop policies and management methods to reduce
saltwater intrusion in coastal aquifers [13].

Section 3 of the study employs an analytical technique [14]
to address the issue raised in Section 2. The Fourier series and
its properties are used as the principal tool for problem analysis
[15, 16]. It was explained in detail how the series’ resulting
factors are implemented in MATLAB to generate a solution.
The Fourier series is a mathematical approach that uses sine and
cosine functions to describe periodic functions as a sum of the
sin and cos functions. It is frequently employed in a variety of
domains, including signal processing, communication systems,
and control theory [17]. The Fourier series’ handling qualities
are critical to solving the problem stated in Section 2, as they
allow the authors to analyse the system’s periodic nature.

2. Formulation of free surface case problem

On the island, there is a porous medium with dimensions of
-L < x < L , which contains fluid of varying densities, specifi-
cally two layers of freshwater and saltwater. The properties of

these fluids can be described using the following equations.

Ωi = pi + ρigy, y = ξ(x), (1)

where i = 1, 2 and y = ξ (x). Equation (1) represents the sum
of total pressures, which remains constant even when gravity
changes. This is called the compression of heads in the two
layers. The variables, ρi andpi represent the density and pres-
sure of each layer of fluid, respectively. The interface between
the two fluids is denoted by y = ξ (x), and the surface of the
less dense fluid is assumed to be in contact with air.

Since the air is assumed to be in a stable state, the pres-
sure at the horizontal surface of the less dense fluid must be
equal to the atmospheric pressure, which is denoted by p1 = 0.
Additionally, we assume that there is no flow through the fluid
surface interface and no flow in the lower layer, i.e.,

Ωi = ρ1gy. (2)

When dealing with a non-persistent pressure situation, it is
physically impossible to obtain a fixed solution. Thus, the pres-
sure along the interface of two fluid layers must match, specifi-
cally at the point where y = ξ (x). This means that the pressures
of both fluid layers, denoted by p1 and p2, must be equal i.e.,

p1 = p2, (3)

at this interface point. This condition ensures that the fluids are
in hydrostatic equilibrium and that there is no net force acting
on the interface. Assuming no fluid flow across the surface in-
terface and in the lower layer (i.e., Ω2 = 0) and with y = ξ (x),
we can obtain:

p2 = −ρ2gy. (4)

By matching the pressures at the fluid interface, we can con-
clude the following:

Ω1 = (ρ1 −ρ2) gy, (5)

and by dividing (5) by ρ1g to be Ω1 = ρ1gzω1 we will have the
following equation:

ω1 = (1 −ψ). (6)

In order to simulate the process of extracting liquids from a
freshwater lens on a tropical island, a tub was positioned at co-
ordinates (xs, ys). To achieve this, we applied the following
formula to Ω1:

Ω1 →
m
4π

ln[(x−xs)
2

+ (y + ys)
2]. (7)

We then applied the equation Ω1 = p1gzω1, taking into account
that µ = zp1 g as (x, y) approaches (xs, ys).

Ω1 →
µ

4π
ln[(x−xs)

2
+ (y + ys)

2]. (8)

Assuming that the tub pressure force and density ratio between
the two layers are represented by the balanced ratios of L and
Ψ, and the island’s aspect ratio is represented by the balanced
ratio of µ, we can rewrite the given equations as follows: Using
the expression ∇Ω1 =

ρ1 gz
z ∇ω1, where q = −κ∇ω1 represents

2



AL-Ali & AL-Awsi / J. Nig. Soc. Phys. Sci. 5 (2023) 1512 3

the speed, we can use Darcy’s law to find the relation between
q, z, and ∇ω1:

q
z
T

= −κρ1g∇ω1. (9)

Since T = z
ρ1 g

, we can simplify equation (9) as follows:

q = −∇ω1. (10)

Equation (10) shows that q is directly proportional to −∇ω1.
The parameters L, ψ, and µ have balanced ratios that represent
the tub pressure force, the density ratio between the two layers,
and the island’s aspect ratio, respectively.

3. Freshwater lens on an island with No-withdrawal

To solve the linear problem, Laplace’s equation is used with
the conditions specified in equation

∇
2ω1 =

∂2ω1

∂x2
+
∂2ω1

∂y2
= 0 −L < x < L, ζ (x) < y < 1, (11)

which is subject to a range of conditions:ω1 = 1, y = 1 −L < x < Lω1 = y, x = L 0 < y < 1.
We assume ω1 = ω and use the formula given in equation

ω (x, y) = y+
N∑

k=0

S k sinh
[(

k
L

)
π (y − 1)

]
sin

[(
k
L

)
πk

]
.(12)

and equation ω = (1 −ψ) y = 1 y = ζ (x)ωy − ζ′ (x) ωx = 0, (13)
defines additional conditions that must be met, including ω =
(1 −ψ)y = 1 and ωy = 0 at y = 0, where we take a small value
of ζ assume a large value of ψ.

ωy = 1+
N∑

k=0

S k

(
k
L

)
π cosh

[(
k
L

)
π (y − 1)

]
sin

[(
k
L

)
πx

]
.(14)

At y = 0, equation (14) becomes:

−1 =
N∑

k=0

S k

(
k
L

)
π cosh

(
−

k
L

)
π sin

[(
k
L

)
πx

]
. (15)

To apply Orthogonal, we multiply both sides of (15) by
sin

(
j
L

)
πx, integrate over the range−L to L, and get:

∫ L
−L

sin
( j

L

)
πx d x =

∫ L
−L

N∑
k=0

−S k

(
k
L

)
π cosh

(
−

k
L

)
π sin2

( j
L

)
πx d x, (16)

assume that Ak = S k
(

k
L

)
π cosh

(
−

k
L

)
π, i.e.

S k =
2[(−1)k − 1]

kπ
.

From which we can conclude that:

Ak =
2L[(−1)k − 1]
(kπ)2 cosh( kL )π

. (17)

Therefore, the Fourier series approximation is:

Ω1 (x, y) = y +
N∑

k=0

S k

[(
k
L

)
π (y − 1)

]
sin

[(
k
L

)
πk

]
. (18)

Using (13) and evaluating it at y = 0, we get an approximation
of ζ = ω1−ψ :

ζ =
ω

1 −ψ

 N∑
k=0

S k shin
(

k
L

)
π sin

[(
k
L

)
πx

] . (19)
The linear island problem involves the interface between fresh-
water and saltwater, and the depth of this interface is determined
by the density rate. The Ghyben-Herzberg relation, which is
based on hydrological principles, provides a formula for deter-
mining the depth of the interface in a system that is in a constant
state of balance. The depth of the interface, according to this
relationship, is proportional to the ratio of the freshwater head
to the total head, where the total head is the sum of the fresh-
water head and the depth of the interface below sea level [18].
Therefore, if the freshwater head increases, the depth of the in-
terface will also increase, and if the freshwater head decreases,
the depth of the interface will decrease as well. The resulting
approximation of ζ is given in Figure (2) shows the resulting ap-
proximation of ζ for an island with width L=100 and a density
ratio of 1.01. The measurements were conducted with varying
intensity ratios, and it was discovered that the best measure-
ment was achieved at a density value of ψ = 1.01 with n = 300
points. To analyse the interface between freshwater and saltwa-

Figure 2. A Diagram determine w-plane and z-plane for free surface
condition

ter in the linear island problem, researchers use the Fourier se-
ries approximation, a mathematical technique that breaks down

3



AL-Ali & AL-Awsi / J. Nig. Soc. Phys. Sci. 5 (2023) 1512 4

complex functions into simpler functions that can be more eas-
ily studied. According to the Fourier series approximation, the
depth of the interface is inversely proportional to the density
rate, which means that as the density rate increases, the inter-
face depth decreases, and vice versa. The Ghyben-Herzberg re-
lation, which states that ζ = ω1−ψ , also supports this relationship
between the interface depth and the density rate. In a system
in constant balance, when freshwater is present above the in-
terface, the depth of the interface above sea level is equal to
the depth of the interface below sea level. Both the Fourier se-
ries approximation and the Ghyben-Herzberg relation are cru-
cial tools for understanding the interface between freshwater
and saltwater in the linear island problem [19, 30]. These prin-
ciples have practical applications in fields such as hydrology
and groundwater management, where researchers can use them
to better understand groundwater systems and make informed
decisions about resource management and conservation [20].

3.1. Integration of the free surface problem equation

To solve a problem numerically, we need to find an appro-
priate solution that satisfies the governing equation and bound-
ary conditions. In the case of Laplace’s equation, we use
Green’s function F to find a single solution that meets the nec-
essary boundary requirements for the area of concern [21].
Green’s function F can be derived using the free surface con-
dition, which is a solution to the equation:

∇
2 F

(
x, y, x0, y0

)
= δ (x − x0, y − y0) , (20)

where the function is subject to the boundary conditions:

F (±L, y, x0, y0) = 0, 0 < y < 1. (21)

F (x, 1, x0, y0) = 0, −L < y < L. (22)

To solve this equation using integral equations or finite element
methods, we only need to find a solution for on the boundary.
This is achieved by using Green’s second identity and applying
the boundary clauses of equations (21) and (22), which ensures
that the integration line of remains on the boundary [22].

The appropriate Green’s function F for this problem needs
to satisfy the conditions of equations (20), (21), and (22) [23].
One way to find this function is by using the techniques of port-
folio conversions of angles in composite variables and applying
the solution on the w-plane to the physical (real) plane. We
can consider logarithms to be individual at w = w0, where
w0 = u0 − iv0 and w̄0 = u0 − iv0. We need a function that
evaluates to zero on the real axis, and so we test:

F = Re [ln (w − w0) − ln(w − w̄0)] . (23)

There are different methods to find the appropriate Green’s
function F for a particular problem, depending on the bound-
ary conditions and the governing equation. Some examples of
these methods include the method of images, separation of vari-
ables, and the method of integral equations.

f (z) = A
∫ z

0
(τ− x1)

2L1 (τ− x2)
2L2 . . . (τ− xn)

2Ln dτ+C.(24)

To convert from the w-plane to the z-plane for each Li as shown
in Figure (3), the conversion should be applied to the inner cor-
ners of the polygon. Functions of form (24) were used to trans-
form the real axis into a polygonal path. Beginning with the
assumption that w = 1 corresponds to z = i + L, the appropriate
values were offset in (23) to obtain the desired outcome.

Figure 3. Diagram determine w-plane and z-plane for free surface con-
dition

z = A
∫ w

1
(w − 1)

−1
2 (w + 1)

−1
2 dw + i + L

= A
∫ w

1
(w − 1)

−1
2 dw + i + L. (25)

And by using the definitions cosh2θ− sinh2θ = 1 we obtain

z = A
∫ w

1
(sinh2θ)

−1
2 sinh θ dθ + i + L. (26)

To clarify the given expression, let us rewrite it as follows: Let
w = cosh

[
z−(i+L)

A

]
. Since cosh(D) cannot be zero unless D

is complex, we can use cosh y = cos iy and set cos (−iLA ) = 0,
which yields A = −iL2π .
Using this, we can rewrite w as:

w = cosh
[
z − (i + L)
−i2Lπ

]
= cosh

[
π (z − i)

2L
−
π

2

]
= cos i

[
π (z − i)

2L
−
π

2

]
= cos i

(
z −

π

2

)
= sin z, (27)

and when z = x + iy we will obtain

w = sin
πx
2L

+ cosh
π(y − 1)

2L

+i sinh
π(y − 1)

2L
cos

πx
2L
. (28)

Using the Schwarz-Christoffel transformation, we can obtain
the following function: [missing function].

F = ln
∣∣∣∣∣ w − w0(w − w̄0

∣∣∣∣∣ , (29)
and now let us say that

f = sin
πx
2L

cosh
π(y − 1)

2L
, (30)

g = sinh
π(y − 1)

2L
cos

πx
2L
, (31)

4



AL-Ali & AL-Awsi / J. Nig. Soc. Phys. Sci. 5 (2023) 1512 5

We will obtain a

w − w0
w − w̄0

=
( f − f0) + i(g − g0)
( f − f0) + i(g + g0)

, (32)

and when you take the real part of F, we find that (25) simplified
to the follows:

F =
1
2

ln
( f − f0)

2
+ (g − g)2

( f − f0)
2

+ (g + g)2
. (33)

3.2. Derivation of the free surface problem equation
We will now derive the complementary equation for the un-

defined interface, which needs to satisfy the following equation,
as mentioned earlier:

∇
2ω1 =

∂2ω1

∂x2
+
∂2ω1

∂y2
= 0 −L < x < L, ζ (x) < y < 1.(34)

These equations are subject to the following boundary condi-
tions:

f (x)

ω1 = 1, y = 1,−L < x < Lω1 = y, x = ±L, 0 < y < 1. (35)
We also impose the condition that the pressure must match
along the interface between saltwater and freshwater where
y = ζ (x). Furthermore, we assume that there is no flow through
the interface, leading to the following conditions:

∂ω1
∂n ,−L < x < L
ω1 = y, x = ±L, 0 < y < 1

(36)

When considering fluid dynamics, it is important to keep in
mind the effect of negative pressure which results from the
pulling of liquids. To address this issue, Green’s second iden-
tity is commonly used as a tool in solving related equations [24,
225].

By utilizing the condition that the Laplacian of a function
∇2ω1 is equal to zero (∇2ω1 = 0), we can simplify the equation
and arrive at the expression below:

ω1 →
µ

2π
ln

√
(x − xs)

2 + (y − ys)
2. (37)

Here, F is a function that satisfies Laplace’s equation, i.e.,
∇2 F = 0, except at a specific point (x0, y0). Along the free
surface, ∂ω/∂n = 0 and similarly, along the other three bound-
aries (top, right side, and left side), F = 0. As a result, the
second term in Equation (1) falls away, leading to:∫

Γ

ω
∂F
∂n
−ω

∂ω

∂n
= 0. (38)

However, we must also consider what happens along the bound-
aries s1, sE1 sE2 , wheres1 is the limit along the boundary and
the top, sE1 is the ring around the source (i.e., where the liquid
passes through), and sE2 is the ring around a single point on the
surface. Therefore, we obtain:∫

S 1
ω
∂F
∂n

d s =
∫

sE1

F
∂ω

∂n
d s =

∫
sE2

ω
∂F
∂n

d s = 0. (39)

Assuming that a single point is enclosed within the ring sE1 ,
which causes all the liquid to pass through it, we can consider
its specific behaviour.

ω →
1

4π
ln

[
( fs − f0)

2
+ i(gs − g0)

2

( fs − f0)
2

+ i(gs + g0)
2

]
. (40)

We start by considering a function of a point (x0, y0). Now,
let’s consider ∂ω

∂n as an integral of sE1 . This integral represents
the speed of the liquid out of the tub. We can express this as
follows:∫

S 1
−F

∂ω

∂n
d s = −F

∫
sE1

∂ω

∂n
d s

= F
∫

sE2

µ

2π
.

1
r

d s =
Fµ
2π

∫ 0
2π

r.
1
r

dθ = Fµ, (41)

where r = [(x−xs)2 + (y−ys)2]
1
2 is the distance from the source

point (xs, ys) to the point (x, y). In other words, the flow from
the tub has force µ, so the contribution from the integration of
sE1 is µF(xs, ys). Therefore, we can write:

µ

4π
ln

[
( f0 − fs)

2
+ (g0 − gs)

2

( f0 − fs)
2

+ (g0 + gs)
2

]
+

∫
s1
ω
∂F
∂n

d s =
∫

sE2

F
∂ω

∂n
d s = 0, (42)

where ( fs, gs) is the location of the source point, and ( f0, g0) is
the location of the point at which we are evaluating the function.
We can think of the integral along the sE1 line as the integration
of a function with a constantly changing value along this line
where (x, y) → (x0, y0).

∫
∂F
∂n d s is ”the flow” of ”the source” due

to the effect of F. Since the loop around the tub point (x0, y0)
is a semicircle, the flow volume due to the source the tub F is
Q = π. Thus, the second integration in (42) could be ωπ(x0, y0).
Therefore, we can write:

µ

4π
ln

[
( f0 − fs)

2
+ (g0 − gs)

2

( f0 − fs)
2

+ (g0 + gs)
2

]
+

∫
s1
ω
∂F
∂n

d s −ωπ(x0, y0) = 0. (43)

Note that there is a single point of integration on s1 such as
(x, y) → (x0, y0), which is the surface point. To deal with this
single point, we extract from the method of collecting and sub-
tracting ω0 under integration, giving us the following formula:

µ

4π
ln

[
( f0 − fs)

2
+ (g0 − gs)

2

( f0 − fs)
2

+ (g0 + gs)
2

]
−ωπ(x0, y0)

+

∫
s1

(ω−ω0)
∂F
∂n

d s + ω0

∫
sE1

∂F
∂n

d s = 0. (44)

4. Outcomes And discussion

The study examined the impact of density ratios on the for-
mation of the interface and depth of the freshwater lens in free

5



AL-Ali & AL-Awsi / J. Nig. Soc. Phys. Sci. 5 (2023) 1512 6

surface models. It used both analytical and numerical methods
to solve the problem of individual free surfaces. The Fourier
series and its orthogonal properties were employed to calculate
the interface without the need for pumping or pulling on the
island. The study also explored the impact of the presence or
absence of a source/tub on the island on the maximum pulling
rates at different density ratios. With no statistically significant
difference between the analytical and numerical solutions, the
results demonstrated that the lowest intensity situation resulted
in the greatest pulling rate. The findings were consistent with
previous research [26, 27]. Finally, the study assessed the ef-
ficiency of the numerical chart and found that the execution of
the program with a low n-value produced results that were just
as accurate as those with a high n-value. The research provides
valuable insights into the impact of density ratios on free sur-
face models and offers practical solutions for calculating inter-
faces and freshwater lens depth. The study’s methodology and
results can be used in various applications, including environ-
mental management and coastal engineering.

5. Conclusion

The study aimed to examine the impact of withdrawing wa-
ter from the freshwater layer of an island that has fixed bound-
aries at its bottom, left, and right. The unknown interface be-
tween two layers of fluids of different densities on an island
with consistent boundaries was calculated by integrating rele-
vant parameters such as the island’s pulling rate and density ra-
tio. The study used an analytical approach based on the Fourier
series to calculate a solution to the linear problem, and it was
found to provide a good solution to the non-linear problem for-
mulated. The height of the calculated interface through the an-
alytical approach was consistent with the height predicted by
the Ghyben-Herzberg relation [29, 30, 31], which describes
the equilibrium interface between freshwater and saltwater in
a coastal aquifer system. Interestingly, the study also found
that there is a maximum pulling rate for different density ra-
tios after which fixed solutions cease to exist. These findings
are consistent with the results of previous studies conducted on
similar systems [32]. The study highlights the importance of
understanding the complex hydrological processes that occur
in coastal aquifer systems and the potential impact of human
activities on these systems.

References

[1] W. Li, X. Chen, L. Xie, Z. Liu & X. Xiong, “ Reply to Comments on: Li
et al.(2019)” Bioelectrochemical Systems for Groundwater Remediation:
“The Development Trend and Research Front Revealed by Bibliometric
Analysis”, Water 11 (2020) 1532. https://doi.org/10.3390/w12061603.

[2] World Health Organization (WHO), United Nations Department of Eco-
nomic and Social Affairs, Progress on sanitation and drinking water: 2015
update and MDG assessment. World Health Organization, 2015.

[3] A. E. Urai & C. Kelly, “Rethinking academia in a time of climate crisis”,
ELife 12 (2023) 84991.

[4] G. Worster, K. Moffatt & G. Batchelor (Eds.) Perspectives in Fluid Dy-
namics: A Collective Introduction to Current Research, Cambridge, UK:
Cambridge University Press, 2000.

[5] World Health Organization, Health and climate change: country
profile 2019: small island developing states initiative: Tuvalu No.
WHO/CED/PHE/EPE/19.3. 3. World Health Organization, 2019.

[6] L. Chambers, S. Lui, R. Plotz, D. Hiriasia, P. Malsale, R. Pulehetoa-
Mitiepo, M. Natapei, N. Sanau, M. Waiwai, L. Tahani & A. Willy, “Tra-
ditional or contemporary weather and climate forecasts: reaching Pacific
communities”, Regional Environmental Change 19 (2019) 1521.

[7] M. A. Mycoo, Achieving SDG 6: “water resources sustainability in
Caribbean Small Island Developing States through improved water gov-
ernance, In Natural Resources Forum (Vol. 42, No. 1, pp. 54-68). Oxford,
UK: Blackwell Publishing Ltd.

[8] Pacific Community, Water Security Toolkit for Pacific Island Countries,
2019.

[9] H. Canton, Canton, Helen, ed. “The Europa Directory of International
Organizations: 2021”. Routledge, 2021.

[10] S. Al-Ali, G. C. Hocking, D. E. Farrow & H. Zhang, “A spectral mod-
elling approach for fluid flow into a line sink in a confined aquifer”, Eu-
ropean Journal of Applied Mathematics 33 (2022) 960.

[11] A. D. Werner, & M. Bakker, “Saltwater intrusion in coastal aquifers”,
Wiley Interdisciplinary Reviews: Water 6 (2019) 1388.

[12] E. Custodio, “Aquifer overexploitation: what does it mean?”. Hydrogeol-
ogy Journal (2002)10:254–277.

[13] L. Martinez-Perez, L. L. Carrera, J. Marazuela, M. A. Goyetche, T. M.
Pool & A. Folch, “ A multidisciplinary approach to characterizing coastal
alluvial aquifers to improve understanding of seawater intrusion and sub-
marine groundwater discharg”, Journal of Hydrology, 607,(2022)127510.

[14] M. O. Ogunniran “A Class of Block Multi-derivative Numerical Tech-
niques for Singular Advection Equations”, Journal of the Nigerian Soci-
ety of Physical Sciences 1 (2019) 62

[15] S. Al-Ali, G. C. Hocking, & D. E. Farrow “Critical surface coning due to
a line sink in a vertical drain containing a porous medium”, The ANZIAM
Journal 61 (2019) 249.

[16] O. O. Olaiya, R. A. Azeez, & M. I. Modebei, “Efficient Hybrid Block
Method For The Numerical Solution Of Second-order Partial Differential
Problems via the Method of Lines”. Journal of the Nigerian Society of
Physical Sciences, (2021) 26.

[17] R. N. Bracewell, & R. N. Bracewell, The Fourier transform and its appli-
cations New York: McGraw-Hill Vol. 31999, 1986.

[18] C. K. Wentworth, “Storage consequences of the Ghyben-Herzberg theor”,
Eos, Transactions American Geophysical Union 23 (1942) 683.

[19] W. Badon-Ghyben, Notes on the probable results of well drilling near
Amsterdam, Tijdschrift van het Koninklijk Inst. van Ing. Den Haag, 1888.

[20] D. K. Todd, & L. W. Mays, Groundwater Hydrology, Third Edition. John
Wiley and Sons, Inc. (2005).

[21] P. K. Kundu, I. M. Cohen & D. R. Dowling, Fluid mechanics, Academic
press, 2015.

[22] L. N. Trefethen & D. Bau,“Numerical linear algebr”, Siam 181 (2022).
[23] H. T. Davis, “Introduction to nonlinear differential and integral equa-

tions”, US Atomic Energy Commission, 1960.
[24] K. Aziz, & A. Settari, Petroleum reservoir simulation, 1979.
[25] C. A. Fletcher, Computational techniques for fluid dynamics 2: Specific

techniques for different flow categories Springer Science & Business Me-
dia, 2012.

[26] Z. Ling, L. Shu, Y. Sun, R. Wang & Y. Li, “Impact of island urbanization
on freshwater lenses: A case study on a small coral island”, Water 13
(2021) 3272.

[27] A. D. Werner & C. T. Simmons, “Impact of sea-level rise on sea water
intrusion in coastal aquifers”, Groundwater 47 (2009) 197.

[28] A. D. Werner, M. Bakker, V. E. Post, A. Vandenbohede, C. Lu, B. Ataie-
Ashtiani & D. A. Barry, “Seawater intrusion processes, investigation and
management: recent advances and future challenges”, Advances in water
resources 51 (2013) 3.

[29] H. I. Essaid, “A comparison of the coupled fresh water-salt water flow and
the Ghyben-Herzberg sharp interface approaches to modeling of transient
behavior in coastal aquifer systems”, Journal of Hydrology 86 (1986) 169.

[30] H. Bouwer & H. Bouwer, Groundwater hydrology, New York: McGraw-
Hill, 1978.

[31] A. J. Raudkivi & R. A. Callander, Analysis of groundwater flow, No.
Monograph, 1976.

[32] W. B. Ghyben, Nota in verband met de voorgenomen putboring nabij
Amsterdam, Tijdschrift van Let Koninklijk Inst. Van Ing., 1888.

6