# Classification of the Group Invariant Solutions for Contaminant Transport in Saturated Soils under Radial Uniform Water Flows

Authors: M. M. Potsane; R. J. Moitsheki
Journal: Advances in Mathematical Physics
 (2014)

---

## Abstract

The transport of chemicals through soils to the groundwater or precipitation at the soils surfaces leads to degradation of these resources. Serious consequences may be suffered in the long run. In this paper, we consider macroscopic deterministic models describing contaminant transport in saturated soils under uniform radial water flow backgrounds. The arising convection-dispersion equation given in terms of the stream functions is analyzed using classical Lie point symmetries. A number of exotic Lie point symmetries are admitted. Group invariant solutions are classified according to the elements of the one-dimensional optimal systems. We analyzed the group invariant solutions which satisfy the physical boundary conditions.


---

## Body

## 1. Introduction

Prediction of water and contaminant movements in soils is important to the theory of soil salinity and underground water pollution. Among others, the sources of contamination are movement of agricultural and industrial contaminants. Some exact solutions are constructed for constraint assumptions; as such, investigation and study of these problems are extremely difficult and challenging [1] even when the problem is given in terms of linear partial differential equations [2].Broadbridge et al. [2–4] applied the classical Lie point symmetries to analyze the convection-dispersion equations arising in solute transport theory. A compendium of exact (group invariant) solutions is given in [5]. Other researchers used different techniques to obtain exact solutions (see, e.g., [6–8]). Chen et al. [9] constructed analytical solutions for two-dimensional advection-dispersion equation with transverse dispersivities depending linearly on the spatial variable. Yadav et al. [10] constructed analytical solutions for solute transport in semiinfinite porous domain. Symmetry reductions and group invariant solutions were constructed for models describing motion of a polytropic gas [11].Exact analytical solutions for contaminant transport in porous media have been constructed, for example, in [12–15]. Chrysikopoulos and Sim [12] developed the one-dimensional stochastic model describing virus transport in homogeneous, saturated semifinite porous media. The ideas in [12] are extended to three-dimensional problems in [13, 14]. Experimental investigations of acoustically enhanced colloid transport in water-saturated packed columns are carried out in [15]. In this paper we focus on theoretical macroscopic deterministic models given in terms of partial differential equations.Numerical models are able to simulate complex reactive transport phenomena but can be time consuming to construct and subject to numerical discretization errors [7]. Also, available packages have significant disagreement in their prediction of solute transport [16]. Therefore, exact solutions are very important because they are needed both as validation tests for numerical schemes and also provide insight into the water and solute transport processes.In this paper, we construct group invariant solutions and classify them according to the one-dimensional optimal systems. This is a significant advancement of the work by [4, 6]. We consider contaminant transport under radial uniform water flows. The resulting convection-dispersion equation given in terms of stream function is analyzed when dispersion coefficient is proportional to the water pore velocity. It turns out that the Lie point symmetries are admitted when the dispersion coefficient is a constant in one case and depending on water pore velocity in the other. In Section 2, we provide the mathematical models. We briefly discuss the symmetry methods in Section 3. The problem at hand is analyzed in Section 4 and the group invariant solutions are classified according to the one-dimensional optimal systems in Section 5 and some physically realistic solutions are discussed. Lastly, some concluding remarks are given in Section 6.


## 2. Mathematical Model

The theoretical background in this section is obtained from Hillel [17]. The macroscopic deterministic model, which is based on local conservation laws, is given by
(1)
∂
(
θ
c
)
∂
t
=
∇
·
(
θ
D
(
v
)
∇
c
)
-
∇
(
c
V
)
,
where D
(
v
) denotes the coefficient of hydrodynamic dispersion and ∇ is given in Cartesian coordinates. Experimental and theoretical observations show that the dispersion coefficient can take the power law form D
(
v
)
=
λ
v
p with λ being a proportionality constant and 1
≤
p
≤
2 (see, e.g., [18]). Here, v is the pore water velocity. Given uniform water flow in saturated soils, then the continuity equation is given by ∇
·
V
=
0. Uniform saturated water flow implies that the soil water content θ
=
θ
s, where θ
s is the water content at saturation. By Darcy’s law, which states that water flux is proportional to the gradient of the hydraulic pressure head, that is, V
=
-
k
s
∇
Φ, where Φ is the total hydraulic pressure head and k
s is the hydraulic conductivity at saturation (see, e.g., [17]), as such we obtain Laplace’s equation ∇
2
Φ
=
0. Under these assumptions, (1) then becomes
(2)
∂
c
∂
t
=
∇
·
(
D
(
v
)
∇
c
)
+
k
∇
Φ
·
∇
c
,
where k
=
k
s
/
θ
s and v
=
|
k
∇
Φ
|. Note that this problem becomes extremely difficult to solve exactly when v must be the modulus of the potential flow velocity field for an incompressible fluid (see, e.g., [2]). However, one may transform (1) from Cartesian coordinates (
x
,
y
) to the streamline coordinates (
ϕ
,
ψ
) using Laplace preserving or conformal transformations [19]; see also [2]. The resulting solute transport equation is given by
(3)
∂
c
∂
t
=
v
2
∇
¯
·
[
D
(
v
)
∇
¯
c
]
+
v
2
∂
c
∂
ϕ
,
where ∇
¯
=
(
∂
/
∂
ϕ
,
∂
/
∂
ψ
). The velocity potential is ϕ and ψ is a conjugate harmonic stream function such that ϕ
x
=
ψ
y and ϕ
y
=
-
ψ
x. Also, the functions ϕ and ψ satisfy the Laplace equation. In radial water flows, the velocity potential is given by ϕ
=
-
(
Q
/
θ
s
)
log
r and the stream function ψ
=
-
(
Q
/
θ
s
)
arctan
(
y
/
x
), where Q is the source strength (pumping rate). Introducing the normalized concentration, velocity potential, and time given by C
=
c
/
c
s
,
ϕ
=
-
log
R, and τ
=
t
/
t
s, respectively, with c
s
,
t
s, and R being the concentration, time at soil saturation, and the distance or radius from the point source, we may write (3) as
(4)
∂
C
∂
τ
=
v
2
∇
¯
·
[
D
(
v
)
∇
¯
C
]
+
v
2
∂
C
∂
ϕ
.
We refer to (4) as the governing equation. Equation (4) is susceptible to symmetry analysis (see, e.g., [3, 4]).


## 3. Algebraic Techniques for Symmetry Reduction

In brief, a symmetry of a differential equation is an invertible transformation of the dependent and independent variables that does not change the original differential equation. Symmetries depend continuously on a parameter and form a group: theone-parameter group of transformations. This group can be determined algorithmically by hand or by computer software programs such as YaLie [20], Reduce [21], and Dimsym [22]. The theory of application of Lie groups to differential equations may be found in texts such as those of [23–25]. Given a second order partial differential equation such as (4) describing contaminant transport under radial water flows, we seek transformations of the form
(5)
τ
¯
=
τ
+
ϵ
ξ
1
(
τ
,
ϕ
,
C
)
+
𝒪
(
ϵ
2
)
,
ϕ
¯
=
ϕ
+
ϵ
ξ
2
(
τ
,
ϕ
,
C
)
+
𝒪
(
ϵ
2
)
,
C
¯
=
C
+
ϵ
η
(
τ
,
ϕ
,
C
)
+
𝒪
(
ϵ
2
)
,
generated by the vector field
(6)
X
=
ξ
1
(
τ
,
ϕ
,
C
)
∂
∂
τ
+
ξ
2
(
τ
,
ϕ
,
C
)
∂
∂
ϕ
+
η
(
τ
,
ϕ
,
C
)
∂
∂
C
.
Note that the transformations in (5) are equivalent to the one-parameter Lie group of transformations that leaves the equation in question invariant. Since (4) is second order, then one may prolong the symmetry generator (6) accordingly. The invariance criterion is then given by
(7)
X
[
2
]
(
Equation
(
4
)
)
|
Equation
(
4
)
=
0
,
where X
[
2
] is the second prolongation given by
(8)
X
2
=
ξ
1
∂
∂
τ
+
ξ
2
∂
∂
ϕ
+
η
∂
∂
C
+
ζ
τ
∂
∂
C
τ
+
ζ
ϕ
∂
∂
C
ϕ
+
ζ
ϕ
ϕ
∂
∂
C
ϕ
ϕ
.Hereζ
τ
,
ζ
ϕ, and ζ
ϕ
ϕ are the extended or prolonged infinitesimals (see, e.g., [24]). The invariance criterion results in the overdetermined system of linear homogeneous partial differential equations known as the determining equations, which may be solved even by interactive programs such as YaLie [20] and Reduce [21].


## 4. Lie Point Symmetry Analysis of (4)

We consider contaminant transport under steady radial water flow in saturated soils. In this case, the relevant normalized point source, the water velocity, and the dispersion coefficient are given byϕ
=
-
log
R
,
v
=
e
ϕ, and D
(
v
)
=
λ
e
p
ϕ, respectively (see, e.g. [3]). One may simply omit the dependence of contaminant concentration on the clockwise polar angle coordinate. In the initial symmetry analysis, (4) with arbitrary p and λ admits the time translation and the scaling of C, that is, X
1
=
∂
/
∂
τ and X
2
=
C
(
∂
/
∂
C
), and the infinite symmetry generator
(9)
X
ϑ
=
ϑ
(
τ
,
ϕ
)
∂
∂
C
,
where ϑ is any solution of (4). Extra symmetries are admitted only when p
=
0 and p
=
2.


### 4.1. Constant Dispersion Coefficient

Givenp
=
0, then the dispersion coefficient becomes a proportionality constant, λ. Given dispersion coefficient as an arbitrary constant λ
≠
0, (4) admits finite two extra symmetries
(10)
X
3
=
2
τ
∂
∂
τ
-
∂
∂
ϕ
,
X
4
=
τ
2
∂
∂
τ
-
τ
∂
∂
ϕ
-
1
λ
(
λ
τ
-
τ
2
+
e
-
2
ϕ
4
)
C
∂
∂
C
.
Further, symmetries may be admitted when λ
=
1 and λ
=
-
1. Note that λ
=
-
1 implies negative dispersion coefficient which is not physically realistic, and as such we focus on the case λ
=
1. Given p
=
0 and λ
=
1, then (4) admits finite four extra symmetries given by
(11)
X
3
=
e
ϕ
∂
∂
ϕ
,
X
4
=
2
τ
∂
∂
τ
-
∂
∂
ϕ
,
X
5
=
τ
∂
∂
ϕ
+
(
e
-
2
ϕ
2
)
C
∂
∂
C
,
X
6
=
τ
2
∂
∂
τ
-
τ
∂
∂
ϕ
-
(
τ
2
+
e
-
2
ϕ
4
)
C
∂
∂
C
.


### 4.2. Velocity-Dependent Dispersion Coefficient

The casep
=
2 is in agreement with solute transport theory. This implies that the dispersion coefficient is now given in terms of the water pore velocity. In this case, (4) with arbitrary proportionality constant λ admits extra four finite Lie point symmetries given by
(12)
X
3
=
e
2
ϕ
∂
∂
ϕ
,
X
4
=
1
λ
(
e
-
4
ϕ
4
-
τ
e
-
2
ϕ
2
)
C
∂
∂
C
+
τ
∂
∂
ϕ
,
X
5
=
τ
∂
∂
τ
-
1
4
∂
∂
ϕ
+
1
λ
(
e
-
2
ϕ
8
-
τ
4
)
C
∂
∂
C
,
X
6
=
τ
2
∂
∂
τ
-
τ
2
∂
∂
ϕ
+
1
λ
(
τ
e
-
2
ϕ
4
-
e
-
4
ϕ
16
-
λ
τ
2
-
τ
2
4
)
C
∂
∂
C
.
Note that the admitted symmetry structure and number are not affected by the constant λ; that is, we obtain the same symmetries up to the specified λ value.


## 4.1. Constant Dispersion Coefficient

Givenp
=
0, then the dispersion coefficient becomes a proportionality constant, λ. Given dispersion coefficient as an arbitrary constant λ
≠
0, (4) admits finite two extra symmetries
(10)
X
3
=
2
τ
∂
∂
τ
-
∂
∂
ϕ
,
X
4
=
τ
2
∂
∂
τ
-
τ
∂
∂
ϕ
-
1
λ
(
λ
τ
-
τ
2
+
e
-
2
ϕ
4
)
C
∂
∂
C
.
Further, symmetries may be admitted when λ
=
1 and λ
=
-
1. Note that λ
=
-
1 implies negative dispersion coefficient which is not physically realistic, and as such we focus on the case λ
=
1. Given p
=
0 and λ
=
1, then (4) admits finite four extra symmetries given by
(11)
X
3
=
e
ϕ
∂
∂
ϕ
,
X
4
=
2
τ
∂
∂
τ
-
∂
∂
ϕ
,
X
5
=
τ
∂
∂
ϕ
+
(
e
-
2
ϕ
2
)
C
∂
∂
C
,
X
6
=
τ
2
∂
∂
τ
-
τ
∂
∂
ϕ
-
(
τ
2
+
e
-
2
ϕ
4
)
C
∂
∂
C
.


## 4.2. Velocity-Dependent Dispersion Coefficient

The casep
=
2 is in agreement with solute transport theory. This implies that the dispersion coefficient is now given in terms of the water pore velocity. In this case, (4) with arbitrary proportionality constant λ admits extra four finite Lie point symmetries given by
(12)
X
3
=
e
2
ϕ
∂
∂
ϕ
,
X
4
=
1
λ
(
e
-
4
ϕ
4
-
τ
e
-
2
ϕ
2
)
C
∂
∂
C
+
τ
∂
∂
ϕ
,
X
5
=
τ
∂
∂
τ
-
1
4
∂
∂
ϕ
+
1
λ
(
e
-
2
ϕ
8
-
τ
4
)
C
∂
∂
C
,
X
6
=
τ
2
∂
∂
τ
-
τ
2
∂
∂
ϕ
+
1
λ
(
τ
e
-
2
ϕ
4
-
e
-
4
ϕ
16
-
λ
τ
2
-
τ
2
4
)
C
∂
∂
C
.
Note that the admitted symmetry structure and number are not affected by the constant λ; that is, we obtain the same symmetries up to the specified λ value.


## 5. Classification of Group Invariant Solutions

In general it is possible to reduce the number of independent variables by one using any linear combination of the admitted base vectors such as those in (11) and (12). In other words, for each l-parameter subgroup (or equivalently l-dimensional subalgebra) of the full symmetry group (or symmetry algebra) there is a corresponding family of group-invariant solutions, which may be infinite [23]. Thus one needs a systematic means to classify these solutions such that none can be derived from the other. In order to classify these solutions one needs to construct a set of optimal systems.Definition 1 (see [23]).
An optimal system ofl-parameter group-invariant solutions to a differential equation (or system of differential equations) is a collection of solutions with the following properties.  (i)
Each solution in the list is invariant under somel-parameter symmetry group of the differential equation (or system of differential equations).
(ii)
If there exists another solution which is invariant underl-parameter symmetry group, then there is a further symmetry generator admitted by the equation (or system) which maps this old solution to the new one.Clearly, an optimal system is a set of elements which lead to symmetry reductions that are not equivalent by any transformation. Since we aim to classify the group invariant solutions of (4), we construct the optimal system for the Lie point symmetries in (11) and (12).


### 5.1. Construction of the Optimal System of Subalgebras

In this section we adopt the method in [23] to construct the one-dimensional system of subalgebras of the algebra spanned by the base vectors in array of (11). To construct the optimal system we first need to determine the commutators of the admitted symmetries.Definition 2 (see [26]).
Given the generators(13)
X
1
=
ξ
1
i
(
x
,
u
)
∂
∂
x
i
+
η
1
(
x
,
u
)
∂
∂
u
,
X
2
=
ξ
2
i
(
x
,
u
)
∂
∂
x
i
+
η
2
(
x
,
u
)
∂
∂
u
,
admitted by a kth-order partial differential equation
(14)
u
(
k
)
=
F
(
x
,
u
,
u
(
1
)
,
…
,
u
(
k
-
1
)
)
,
then the commutator of X
1 and X
2 is defined by
(15)
[
X
1
,
X
2
]
=
(
X
1
(
ξ
2
i
)
-
X
2
(
ξ
1
i
)
)
∂
∂
x
i
+
(
X
1
(
η
2
)
-
X
2
(
η
1
)
)
∂
∂
u
.
One may list the commutators of symmetries in (11) together with the time translation and concentration scaling in Table 1.Table 1
Commutators of the admitted symmetries.
[X
i
,
X
j]
X
1
X
2
X
3
X
4
X
5
X
6
X
1
0
0
0
2
X
1
2
X
3
4
X
4
-
2
X
2
X
2
0
0
0
0
0
0
X
3
0
0
0
X
3
-
X
2
2
X
5
X
4
-
2
X
1
0
-
X
3
0
X
5
2
X
6
X
5
-
2
X
3
0
X
2
-
X
5
0
0
X
6
2
X
2
-
4
X
4
0
-
2
X
5
-
2
X
6
0
0Furthermore we construct a set of one-dimensional subalgebras which are equivalent to a unique element of the set under some element of the adjoint representation given by(16)
Ad
(
e
ϵ
X
i
)
X
j
=
∑
n
=
0
∞
ϵ
n
n
!
(
Ad
X
i
)
n
X
j
=
X
j
-
ϵ
[
X
i
,
X
j
]
+
ϵ
2
2
!
[
X
i
,
[
X
i
,
X
j
]
]
-
⋯
,
where the commutator of X
i and X
j is defined above. The adjoint representation of base vectors in (11) is given in Table 2. Let
(17)
X
=
a
1
X
1
+
a
2
X
2
+
a
3
X
3
+
a
4
X
4
+
a
5
X
5
+
a
6
X
6
.
It remains to use the adjoint table to simplify as much as possible the constants in (17). The key here is the recognition of the function η
(
X
)
=
a
4
2
-
4
a
1
a
6, which is invariant under the full adjoint action [23]. To begin the simplification we concentrate on the constants a
1
,
a
4
,
a
6. If X is given as in (17), then
(18)
X
~
=
∑
i
=
1
6
a
i
~
X
i
=
Ad
(
exp
(
α
X
1
)
)
∘
Ad
(
exp
(
β
X
6
)
)
X
has the coefficients
(19)
a
1
~
=
a
1
-
2
α
a
4
+
4
α
2
a
6
,
a
4
~
=
a
4
-
4
α
a
6
+
4
β
(
a
1
-
2
α
a
4
+
4
α
2
a
6
)
,
a
6
~
=
a
6
+
2
β
(
a
4
-
4
α
a
6
)
+
4
β
2
(
a
1
-
2
α
a
4
+
4
α
2
a
6
)
.
Three cases arise.Table 2
Adjoint representation for the base vectors.
Ad
X
1
X
2
X
3
X
4
X
5
X
6
X
1
X
1
X
2
X
3
-
2
ϵ
X
1
+
X
4
-
2
ϵ
X
3
+
X
5
4
ϵ
2
X
1
+
2
ϵ
X
2
-
4
ϵ
X
4
+
X
6
X
2
X
1
X
2
X
3
X
4
X
5
X
6
X
3
X
1
X
2
X
3
-
ϵ
X
3
+
X
4
ϵ
X
2
+
X
5
-
ϵ
2
X
2
+
X
6
-
2
ϵ
X
5
X
4
e
2
ϵ
X
1
X
2
e
ϵ
X
3
X
4
e
-
ϵ
X
5
e
-
2
ϵ
X
6
X
5
X
1
+
2
ϵ
X
3
-
ϵ
2
X
2
X
2
-
ϵ
X
2
+
X
3
X
4
+
ϵ
X
5
X
5
X
6
X
6
X
1
-
2
ϵ
X
2
+
4
ϵ
X
4
+
4
ϵ
2
X
6
X
2
X
3
+
2
ϵ
X
5
X
4
+
2
ϵ
X
6
X
5
X
6Case 1.
Ifη
(
X
)
>
0, then we choose α to be a real root of a
1
-
2
α
a
4
+
4
α
2
a
6
=
0 and β
=
a
6
/
(
8
α
a
6
-
2
a
4
). This implies a
1
~
=
a
6
~
=
0 and a
4
~
=
η
(
X
)
≠
0. Thus, X is equivalent to the vector X
~
=
X
4
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
5 and X
3, namely, Ad
(
exp
(
-
a
5
~
X
5
)
) and Ad
(
exp
(
a
3
~
X
3
)
), as such X
5 and X
3 in X
~ vanish. No further simplifications are possible; therefore X is equivalent to a multiple of X
4
+
a
X
2 for some a
∈
ℝ provided η
(
X
)
>
0.Case 2.
Ifη
(
X
)
<
0, we set α
=
0 and β
=
-
a
4
/
4
a
1 so that a
4
~
=
0. One may then assume both the coefficients of X
1 and X
6 to be unity. Thus, X is equivalent to the vector X
~
=
X
1
+
X
6
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
5 and X
3, namely, Ad
(
exp
(
(
-
a
3
~
/
2
)
X
5
)
) and Ad
(
exp
(
(
a
5
~
/
2
)
X
3
)
), as such X
5 and X
3 in X
~ vanish. No further simplifications are possible; as such X is equivalent to X
1
+
X
6
+
a
X
2
,
a
∈
ℝ given η
(
X
)
<
0.Case 3.
Considerη
(
X
)
=
0. Two subcases arise. If not all the coefficients a
1, a
4, a
6 are zero, then we are free to choose α and β such that a
1
~
≠
0 and a
4
~
=
a
6
~
=
0. In this case X is equivalent to a multiple of X
~
=
X
1
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
1 and X
3, namely, Ad
(
exp
(
(
a
3
~
/
2
a
5
~
)
X
5
)
) and Ad
(
exp
(
(
-
a
2
~
/
a
5
~
)
X
3
)
)
,
a
5
~
≠
0, as such X
2 and X
3 in X
~ vanish. Therefore X is equivalent to a multiple of X
1
±
X
5. If a
5
~
=
0, then acting on X
~ by Ad
(
exp
(
-
(
a
3
~
/
2
)
X
5
)
), so X is equivalent to a multiple of X
1
+
a
X
2
,
a
∈
ℝ. If all the coefficients a
1
,
a
4
,
a
6 are zero, assuming a
3
≠
0 (say a
3
=
1), then acting on X by Ad
(
exp
(
a
2
~
X
5
)
) and Ad
(
exp
(
(
-
a
5
~
/
2
)
X
5
)
) yields X
3, a multiple of X. If a
3
=
0 but a
5
≠
0, acting by any group generated by X
1, namely, Ad
(
exp
(
ϵ
X
1
)
), gives a nonzero coefficient in front of X
3 implying that this case is similar to the case when a
3
≠
0. Thus, the only remaining vectors are multiples of X
2. The set of one-dimensional optimal system is given by
(20)
{
X
4
+
a
X
2
,
X
1
+
X
6
+
a
X
2
,
X
1
±
X
5
,
X
1
+
a
X
2
,
X
2
,
X
3
}
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
a
∈
ℝ
.Repeating these calculations, we construct the one-dimensional optimal system for the Lie point symmetries in (12) together with the time translation and the scaling of C and obtain
(21)
{
X
5
+
a
X
2
,
X
1
+
X
6
+
a
X
2
,
X
1
±
X
4
,
X
1
+
a
X
2
,
X
2
,
X
3
}
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
a
∈
ℝ
.


### 5.2.   Construction of Group-Invariant Solutions

The group-invariant solutions are constructed in this section. We further classify the group-invariant solutions according to the elements of the one-dimensional optimal systems. The reductions and group invariant solutions are listed in Tables3 and 4.Table 3
Group invariant solutions for (4) with constant D.
Symmetry
Reductions
a
∈
ℝ
Group invariant solutions
X
4
+
a
X
2
C
=
τ
a
/
2
G
(
γ
), where γ
=
τ
e
2
ϕ and G satisfies γ
3
G
′
′
+
1
4
γ
(
6
γ
-
1
)
G
′
′
-
1
8
a
G
=
0.
a
<
0
C
=
τ
a
/
2
[
1
2
1
τ
e
ϕ
k
2
1
F
1
(
a
+
1
2
,
3
2
;
-
1
4
τ
e
2
ϕ
)
+
k
1
1
F
1
(
a
2
,
1
2
;
-
1
4
τ
e
2
ϕ
)
].
a
=
0
C
=
[
k
2
-
2
π
k
1
erf
(
1
2
τ
e
2
ϕ
)
].
a
>
0
C
=
τ
a
/
2
[
1
2
1
τ
e
ϕ
k
2
1
F
1
(
1
-
a
2
,
3
2
;
-
1
4
τ
e
2
ϕ
)
+
k
1
1
F
1
(
-
a
2
,
1
2
;
-
1
4
τ
e
2
ϕ
)
].
X
1
+
a
X
2
C
=
e
a
τ
G
(
γ
), where γ
=
ϕ and G satisfies
G
′
′
+
G
′
-
a
e
-
2
ϕ
G
=
0.
a
<
0
C
=
e
a
τ
[
k
1
cos
a
e
-
ϕ
-
k
2
sin
a
e
-
ϕ
].
a
=
0
C
=
-
e
-
ϕ
k
1
+
k
2.
a
>
0
C
(
τ
,
ϕ
)
=
e
a
τ
[
k
1
cosh
(
a
e
-
ϕ
)
-
i
k
2
sinh
(
a
e
-
ϕ
)
].
X
1
+
X
5
C
=
e
(
τ
e
-
ϕ
+
(
2
/
3
)
τ
3
)
G
(
γ
), where γ
=
τ
2
+
e
-
ϕ and G satisfies
G
′
′
-
γ
G
=
0.
—
C
=
exp
(
τ
e
-
ϕ
+
2
3
τ
3
)
(
Ai
(
τ
2
+
e
-
ϕ
)
k
1
+
Bi
(
τ
2
+
e
-
ϕ
)
k
2
).
X
1
+
X
6
+
a
X
2
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
tan
-
1
(
2
τ
)
]
G
(
γ
),  where γ
=
e
ϕ
1
+
4
τ
2 and G satisfies γ
6
G
′
′
+
2
γ
5
G
′
+
(
1
-
a
γ
2
)
G
=
0.
a
<
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
tan
-
1
(
2
τ
)
]
w
w
w
w
w
w
w
w
w
w
×
[
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
2
1
F
1
(
3
+
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
w
w
w
w
w
w
w
w
w
w
w
w
w
+
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
1
U
(
3
+
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
]
.
a
=
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
]
w
w
w
w
w
w
w
w
=
×
[
1
2
e
ϕ
1
+
4
τ
2
J
-
1
/
4
(
1
2
e
2
ϕ
(
1
+
4
τ
2
)
)
k
1
Γ
(
3
4
)
w
w
w
w
w
w
w
w
e
e
w
-
+
1
2
e
ϕ
1
+
τ
2
J
1
/
4
(
1
2
e
2
ϕ
(
1
+
4
τ
2
)
)
k
2
Γ
(
5
4
)
]
.
a
>
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
arctan
(
2
τ
)
]
w
w
w
w
w
w
w
w
w
×
[
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
2
1
F
1
(
3
-
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
w
w
w
w
w
w
w
w
w
e
w
=
+
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
1
U
(
3
-
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
]
.Table 4
Group invariant solutions for (4) given D
(
v
)
=
λ
e
2
ϕ.
Symmetry
Reductions
a
∈
ℝ
Group invariant solutions
X
5
+
a
X
2
C
=
exp
[
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
8
λ
ln
(
τ
)
]
G
(
γ
), where γ
=
e
ϕ
τ
4 and G satisfies
λ
2
γ
6
G
′
′
-
1
4
γ
λ
(
1
-
12
γ
λ
4
)
G
′
-
a
8
G
=
0.
a
<
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
ln
(
τ
)
8
λ
)
×
[
k
2
e
-
2
ϕ
4
λ
τ
1
F
1
(
4
+
a
8
λ
,
3
2
;
-
e
-
4
ϕ
16
λ
τ
)
+
k
1
1
F
1
(
a
8
λ
,
1
2
;
-
e
-
4
ϕ
16
λ
τ
)
]
.
a
=
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
)
[
k
2
-
λ
π
k
1
erf
(
e
-
2
ϕ
4
λ
τ
)
]
.
a
>
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
8
λ
ln
(
τ
)
)
W
W
W
W
W
W
W
W
W
W
W
W
W
W
W
W
×
[
k
2
e
-
2
ϕ
4
λ
τ
1
F
1
(
4
-
a
8
λ
,
3
2
;
-
e
-
4
ϕ
16
λ
τ
)
+
k
1
1
F
1
(
-
a
8
λ
,
1
2
;
-
e
-
4
ϕ
16
λ
τ
)
]
.
X
1
+
X
4
C
=
exp
[
8
λ
τ
3
3
+
τ
e
-
2
ϕ
-
τ
2
]
G
(
γ
), where γ
=
4
τ
2
+
e
-
2
ϕ and G satisfies
G
′
′
-
1
2
λ
G
′
-
γ
4
λ
G
=
0.
—
C
=
exp
(
8
λ
τ
3
3
+
τ
e
-
2
ϕ
-
τ
2
)
×
{
e
4
τ
2
+
e
-
2
ϕ
4
λ
[
k
1
Ai
(
1
+
4
(
4
τ
2
+
e
-
2
ϕ
)
(
4
λ
)
4
/
3
)
W
W
W
W
W
W
W
+
k
2
Bi
(
1
+
4
(
4
τ
2
+
e
-
2
ϕ
)
(
4
λ
)
4
/
3
)
]
}
.
X
1
+
a
X
2
C
=
e
a
τ
G
(
γ
), where γ
=
ϕ and G satisfies
G
′
′
+
(
2
+
e
-
2
ϕ
λ
)
G
′
+
a
λ
e
-
4
ϕ
G
=
0.
a
<
0
C
(
τ
,
ϕ
)
=
e
a
τ
[
exp
(
-
1
4
λ
(
-
1
-
1
-
4
a
λ
)
e
-
2
ϕ
)
k
1
+
exp
(
-
1
4
λ
(
-
1
+
1
-
4
a
λ
)
e
-
2
ϕ
)
k
2
].
a
=
0
C
=
[
k
2
+
exp
(
e
-
2
ϕ
2
)
k
1
].
a
>
0
C
=
e
a
τ
exp
[
(
1
4
λ
(
1
-
1
+
4
a
λ
)
e
-
2
ϕ
)
k
1
+
exp
(
1
4
λ
(
1
+
1
+
4
a
λ
)
e
-
2
ϕ
)
k
2
].
X
1
+
X
6
+
a
X
2
Leading to reductions which are unsolvableDefinition 3 (see [24]).
u
=
H
(
x
) is a group-invariant solution of a kth-order PDE corresponding to the appropriate admitted symmetry generator if and only if u
=
H
(
x
) satisfies
(22)
X
(
u
-
H
(
x
)
)
=
0
,
when
u
=
H
(
x
)
.
That is,
(23)
ξ
i
(
x
,
H
(
x
)
)
∂
H
∂
x
i
=
η
(
x
,
H
(
x
)
)
.
The solution of (23) is obtained from the characteristics equation
(24)
d
ξ
1
d
x
1
=
d
ξ
2
d
x
2
=
⋯
=
d
η
u
.
Classifications of the group-invariant solutions by the elements of the optimal system in (20) and (21) are listed in Tables 3 and 4, respectively. Wherever they appear, k
1 and k
2 are arbitrary constants, 1
F
1
(
b
,
c
;
z
) and U
(
b
,
c
;
z
) are the confluent hypergeometric functions, Ai
(
z
) and Bi
(
z
) are the Airy functions, while J
n
(
z
) and Γ
(
z
) represent the Bessel function of  the first kind and Euler gamma function, respectively, and erf
(
z
) represents an error function. A well-documented review of such functions is presented by Abramowitz and Stegun [27].


### 5.3.  Some Physical Examples

#### 5.3.1. Given Constant Dispersion Coefficient

Example  a. Suppose a concentration C
0 of a solute is supplied to a single point in an instant of time. We require to determine the subsequent concentration of the pollution at various distances from where it was released (see, e.g., [28]). We would expect concentration to vanish at large distance; that is,
(25)
C
(
ϕ
,
τ
)
⟶
0
,
as
ϕ
⟶
-
∞
.The generatorX
4 in (11) (also, an element of the one-dimensional optimal systems given a
=
0) leads to an invariant solution in functional form given by
(26)
C
=
G
(
γ
)
,
where γ
=
τ
e
ϕ and G satisfies the equation
(27)
2
γ
3
G
′′
+
(
4
γ
2
-
1
)
G
′
=
0
,
and hence
(28)
G
=
c
1
+
c
2
erf
(
1
2
γ
)
.
Here erf
(
z
)
=
(
2
/
π
)
∫
0
z
‍
e
-
ν
2
d
ν is the error function [27]. In terms of the original variable and subject to the boundary conditions, we obtain
(29)
C
=
C
0
erfc
(
e
-
ϕ
2
τ
)
.
Here, the erfc(·) is the complement error function defined by (
1
-
erf
(
·
)
). Solution (29) is depicted in Figure 1.Figure 1
Contaminant concentration profile along the radius. HereC
0
=
5.Total flux acrossR
=
R
a is given by
(30)
C
+
1
R
∂
C
∂
R
|
R
=
R
a
=
C
0
erfc
(
R
a
2
τ
)
-
C
0
e
-
R
a
2
/
4
τ
π
τ
R
a
.
The contaminant flux (30) is depicted in Figure 2. Total flux across R
a
=
4 increases and flattens at large time.Figure 2
Solute flux across a fixed radius. HereR
a
=
4 and C
0
=
10.Example  b. The X
6-invariant solution is given in functional form as
(31)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
G
(
γ
)
,
where
(32)
γ
=
τ
e
ϕ
,
G
satisfies
the
ODE
γ
G
′′
+
G
′
=
0
.
We impose the boundary conditions
(33)
C
⟶
0
,
ϕ
⟶
-
∞
,
C
=
ω
(
τ
)
,
ϕ
⟶
∞
.
Infinite concentration at the origin implies that there is a high supply of contaminants at this point. Furthermore, contaminant concentration vanishes when time evolves. In terms of the original variables we obtain the exact (group-invariant) solution given by
(34)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
.Solution (34) is depicted in Figures 3 and 4. In Figure 3, a sharp peak of concentration is observed shortly after τ
=
0 and decreases at later stage. This may be interpreted as an injection of contaminants at a single point; that is, the concentration at a single point increases but due to diffusion at larger time it smoothes out. Note that here we have restricted our analysis using symmetry generator X
6. This symmetry generator leads to simpler and realistic exact solution. In Figure 4, we observe that concentration at the origin decreases with time. Furthermore, this concentration vanishes at large distances.Figure 3
Contaminant concentration profile as time evolves.Figure 4
Contaminant concentration profile at various fixed times.


#### 5.3.2. Given Velocity-Dependent Dispersion Coefficient

Example
a (steady-state solution). The time translation leads to the analysis of the steady-state contaminant transport. Steady-state solutions may be constructed and subjected to the following imposed boundary conditions:
(35)
C
=
C
0
,
R
=
0
,
(36)
C
+
1
R
d
C
d
R
=
0
,
R
=
R
a
.
The boundary condition (35) implies that pollutants are supplied at the origin, and boundary condition (36) correspond to the assumption that pollutants are not carried though at some distance R
a, rather it accumulates here. We obtain the exact solution
(37)
C
=
C
0
{
(
1
-
1
Δ
)
+
1
Δ
exp
(
R
p
λ
p
)
}
,
where Δ is given by
(38)
Δ
=
1
-
(
1
+
R
a
p
-
2
λ
)
exp
(
R
a
p
λ
p
)
.
The solution (37) is depicted in Figures 5, 6, and 7. We observe in Figure 5 that concentration starts decreasing and converges to some value at large distance for p
=
2 than for lower values of p, whereas λ has an opposite effect as shown in Figure 6.Figure 5
Steady contaminant profile for concentration given in (37). Here λ
=
1.Figure 6
Steady contaminant profile for concentration given in (37). Here R
=
R
a
=
2.Figure 7
Effects ofλ and p on concentration profile for solution given in (37). Here R
=
R
a
=
2.Example
b (transient-state solution). It is quite difficult to construct exact solutions for transient contaminant transport subject to these boundary conditions (35) and (36). However, if one assumes that at an initial time, say τ
=
1, the concentration at the point source is given by a constant and that this concentration vanishes at large distances and prolonged periods, then using the symmetry combination X
1
+
a
X
2 from Table 4, the group invariant solution is given by
(39)
C
=
C
0
e
a
τ
exp
{
(
1
-
1
-
4
a
λ
4
λ
)
R
2
}
,
∀
a
<
0
,
λ
>
0
.
Solution (39) is depicted in Figures 8 and 9.Figure 8
Concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and λ
=
1Figure 9
Effects ofλ on the concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and τ
=
2.In his work, Philip [6] considered the instantaneous point source for contaminant dispersion during radial water flow in porous media. Exact solutions were constructed for the two- and three-dimensional models with dispersion coefficient depending on Péclet number. Here, we consider models in stream functions coordinates. Exact close-form (similarity) solutions are constructed using the elements of the one-dimensional optimal systems. These new solutions may be viewed as representing the continual supply of contaminant at a point (source) which are dispersed radially.


## 5.1. Construction of the Optimal System of Subalgebras

In this section we adopt the method in [23] to construct the one-dimensional system of subalgebras of the algebra spanned by the base vectors in array of (11). To construct the optimal system we first need to determine the commutators of the admitted symmetries.Definition 2 (see [26]).
Given the generators(13)
X
1
=
ξ
1
i
(
x
,
u
)
∂
∂
x
i
+
η
1
(
x
,
u
)
∂
∂
u
,
X
2
=
ξ
2
i
(
x
,
u
)
∂
∂
x
i
+
η
2
(
x
,
u
)
∂
∂
u
,
admitted by a kth-order partial differential equation
(14)
u
(
k
)
=
F
(
x
,
u
,
u
(
1
)
,
…
,
u
(
k
-
1
)
)
,
then the commutator of X
1 and X
2 is defined by
(15)
[
X
1
,
X
2
]
=
(
X
1
(
ξ
2
i
)
-
X
2
(
ξ
1
i
)
)
∂
∂
x
i
+
(
X
1
(
η
2
)
-
X
2
(
η
1
)
)
∂
∂
u
.
One may list the commutators of symmetries in (11) together with the time translation and concentration scaling in Table 1.Table 1
Commutators of the admitted symmetries.
[X
i
,
X
j]
X
1
X
2
X
3
X
4
X
5
X
6
X
1
0
0
0
2
X
1
2
X
3
4
X
4
-
2
X
2
X
2
0
0
0
0
0
0
X
3
0
0
0
X
3
-
X
2
2
X
5
X
4
-
2
X
1
0
-
X
3
0
X
5
2
X
6
X
5
-
2
X
3
0
X
2
-
X
5
0
0
X
6
2
X
2
-
4
X
4
0
-
2
X
5
-
2
X
6
0
0Furthermore we construct a set of one-dimensional subalgebras which are equivalent to a unique element of the set under some element of the adjoint representation given by(16)
Ad
(
e
ϵ
X
i
)
X
j
=
∑
n
=
0
∞
ϵ
n
n
!
(
Ad
X
i
)
n
X
j
=
X
j
-
ϵ
[
X
i
,
X
j
]
+
ϵ
2
2
!
[
X
i
,
[
X
i
,
X
j
]
]
-
⋯
,
where the commutator of X
i and X
j is defined above. The adjoint representation of base vectors in (11) is given in Table 2. Let
(17)
X
=
a
1
X
1
+
a
2
X
2
+
a
3
X
3
+
a
4
X
4
+
a
5
X
5
+
a
6
X
6
.
It remains to use the adjoint table to simplify as much as possible the constants in (17). The key here is the recognition of the function η
(
X
)
=
a
4
2
-
4
a
1
a
6, which is invariant under the full adjoint action [23]. To begin the simplification we concentrate on the constants a
1
,
a
4
,
a
6. If X is given as in (17), then
(18)
X
~
=
∑
i
=
1
6
a
i
~
X
i
=
Ad
(
exp
(
α
X
1
)
)
∘
Ad
(
exp
(
β
X
6
)
)
X
has the coefficients
(19)
a
1
~
=
a
1
-
2
α
a
4
+
4
α
2
a
6
,
a
4
~
=
a
4
-
4
α
a
6
+
4
β
(
a
1
-
2
α
a
4
+
4
α
2
a
6
)
,
a
6
~
=
a
6
+
2
β
(
a
4
-
4
α
a
6
)
+
4
β
2
(
a
1
-
2
α
a
4
+
4
α
2
a
6
)
.
Three cases arise.Table 2
Adjoint representation for the base vectors.
Ad
X
1
X
2
X
3
X
4
X
5
X
6
X
1
X
1
X
2
X
3
-
2
ϵ
X
1
+
X
4
-
2
ϵ
X
3
+
X
5
4
ϵ
2
X
1
+
2
ϵ
X
2
-
4
ϵ
X
4
+
X
6
X
2
X
1
X
2
X
3
X
4
X
5
X
6
X
3
X
1
X
2
X
3
-
ϵ
X
3
+
X
4
ϵ
X
2
+
X
5
-
ϵ
2
X
2
+
X
6
-
2
ϵ
X
5
X
4
e
2
ϵ
X
1
X
2
e
ϵ
X
3
X
4
e
-
ϵ
X
5
e
-
2
ϵ
X
6
X
5
X
1
+
2
ϵ
X
3
-
ϵ
2
X
2
X
2
-
ϵ
X
2
+
X
3
X
4
+
ϵ
X
5
X
5
X
6
X
6
X
1
-
2
ϵ
X
2
+
4
ϵ
X
4
+
4
ϵ
2
X
6
X
2
X
3
+
2
ϵ
X
5
X
4
+
2
ϵ
X
6
X
5
X
6Case 1.
Ifη
(
X
)
>
0, then we choose α to be a real root of a
1
-
2
α
a
4
+
4
α
2
a
6
=
0 and β
=
a
6
/
(
8
α
a
6
-
2
a
4
). This implies a
1
~
=
a
6
~
=
0 and a
4
~
=
η
(
X
)
≠
0. Thus, X is equivalent to the vector X
~
=
X
4
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
5 and X
3, namely, Ad
(
exp
(
-
a
5
~
X
5
)
) and Ad
(
exp
(
a
3
~
X
3
)
), as such X
5 and X
3 in X
~ vanish. No further simplifications are possible; therefore X is equivalent to a multiple of X
4
+
a
X
2 for some a
∈
ℝ provided η
(
X
)
>
0.Case 2.
Ifη
(
X
)
<
0, we set α
=
0 and β
=
-
a
4
/
4
a
1 so that a
4
~
=
0. One may then assume both the coefficients of X
1 and X
6 to be unity. Thus, X is equivalent to the vector X
~
=
X
1
+
X
6
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
5 and X
3, namely, Ad
(
exp
(
(
-
a
3
~
/
2
)
X
5
)
) and Ad
(
exp
(
(
a
5
~
/
2
)
X
3
)
), as such X
5 and X
3 in X
~ vanish. No further simplifications are possible; as such X is equivalent to X
1
+
X
6
+
a
X
2
,
a
∈
ℝ given η
(
X
)
<
0.Case 3.
Considerη
(
X
)
=
0. Two subcases arise. If not all the coefficients a
1, a
4, a
6 are zero, then we are free to choose α and β such that a
1
~
≠
0 and a
4
~
=
a
6
~
=
0. In this case X is equivalent to a multiple of X
~
=
X
1
+
a
2
~
X
2
+
a
3
~
X
3
+
a
5
~
X
5. Acting on X
~ by adjoint maps generated by X
1 and X
3, namely, Ad
(
exp
(
(
a
3
~
/
2
a
5
~
)
X
5
)
) and Ad
(
exp
(
(
-
a
2
~
/
a
5
~
)
X
3
)
)
,
a
5
~
≠
0, as such X
2 and X
3 in X
~ vanish. Therefore X is equivalent to a multiple of X
1
±
X
5. If a
5
~
=
0, then acting on X
~ by Ad
(
exp
(
-
(
a
3
~
/
2
)
X
5
)
), so X is equivalent to a multiple of X
1
+
a
X
2
,
a
∈
ℝ. If all the coefficients a
1
,
a
4
,
a
6 are zero, assuming a
3
≠
0 (say a
3
=
1), then acting on X by Ad
(
exp
(
a
2
~
X
5
)
) and Ad
(
exp
(
(
-
a
5
~
/
2
)
X
5
)
) yields X
3, a multiple of X. If a
3
=
0 but a
5
≠
0, acting by any group generated by X
1, namely, Ad
(
exp
(
ϵ
X
1
)
), gives a nonzero coefficient in front of X
3 implying that this case is similar to the case when a
3
≠
0. Thus, the only remaining vectors are multiples of X
2. The set of one-dimensional optimal system is given by
(20)
{
X
4
+
a
X
2
,
X
1
+
X
6
+
a
X
2
,
X
1
±
X
5
,
X
1
+
a
X
2
,
X
2
,
X
3
}
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
a
∈
ℝ
.Repeating these calculations, we construct the one-dimensional optimal system for the Lie point symmetries in (12) together with the time translation and the scaling of C and obtain
(21)
{
X
5
+
a
X
2
,
X
1
+
X
6
+
a
X
2
,
X
1
±
X
4
,
X
1
+
a
X
2
,
X
2
,
X
3
}
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
a
∈
ℝ
.


## 5.2.   Construction of Group-Invariant Solutions

The group-invariant solutions are constructed in this section. We further classify the group-invariant solutions according to the elements of the one-dimensional optimal systems. The reductions and group invariant solutions are listed in Tables3 and 4.Table 3
Group invariant solutions for (4) with constant D.
Symmetry
Reductions
a
∈
ℝ
Group invariant solutions
X
4
+
a
X
2
C
=
τ
a
/
2
G
(
γ
), where γ
=
τ
e
2
ϕ and G satisfies γ
3
G
′
′
+
1
4
γ
(
6
γ
-
1
)
G
′
′
-
1
8
a
G
=
0.
a
<
0
C
=
τ
a
/
2
[
1
2
1
τ
e
ϕ
k
2
1
F
1
(
a
+
1
2
,
3
2
;
-
1
4
τ
e
2
ϕ
)
+
k
1
1
F
1
(
a
2
,
1
2
;
-
1
4
τ
e
2
ϕ
)
].
a
=
0
C
=
[
k
2
-
2
π
k
1
erf
(
1
2
τ
e
2
ϕ
)
].
a
>
0
C
=
τ
a
/
2
[
1
2
1
τ
e
ϕ
k
2
1
F
1
(
1
-
a
2
,
3
2
;
-
1
4
τ
e
2
ϕ
)
+
k
1
1
F
1
(
-
a
2
,
1
2
;
-
1
4
τ
e
2
ϕ
)
].
X
1
+
a
X
2
C
=
e
a
τ
G
(
γ
), where γ
=
ϕ and G satisfies
G
′
′
+
G
′
-
a
e
-
2
ϕ
G
=
0.
a
<
0
C
=
e
a
τ
[
k
1
cos
a
e
-
ϕ
-
k
2
sin
a
e
-
ϕ
].
a
=
0
C
=
-
e
-
ϕ
k
1
+
k
2.
a
>
0
C
(
τ
,
ϕ
)
=
e
a
τ
[
k
1
cosh
(
a
e
-
ϕ
)
-
i
k
2
sinh
(
a
e
-
ϕ
)
].
X
1
+
X
5
C
=
e
(
τ
e
-
ϕ
+
(
2
/
3
)
τ
3
)
G
(
γ
), where γ
=
τ
2
+
e
-
ϕ and G satisfies
G
′
′
-
γ
G
=
0.
—
C
=
exp
(
τ
e
-
ϕ
+
2
3
τ
3
)
(
Ai
(
τ
2
+
e
-
ϕ
)
k
1
+
Bi
(
τ
2
+
e
-
ϕ
)
k
2
).
X
1
+
X
6
+
a
X
2
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
tan
-
1
(
2
τ
)
]
G
(
γ
),  where γ
=
e
ϕ
1
+
4
τ
2 and G satisfies γ
6
G
′
′
+
2
γ
5
G
′
+
(
1
-
a
γ
2
)
G
=
0.
a
<
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
tan
-
1
(
2
τ
)
]
w
w
w
w
w
w
w
w
w
w
×
[
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
2
1
F
1
(
3
+
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
w
w
w
w
w
w
w
w
w
w
w
w
w
+
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
1
U
(
3
+
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
]
.
a
=
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
]
w
w
w
w
w
w
w
w
=
×
[
1
2
e
ϕ
1
+
4
τ
2
J
-
1
/
4
(
1
2
e
2
ϕ
(
1
+
4
τ
2
)
)
k
1
Γ
(
3
4
)
w
w
w
w
w
w
w
w
e
e
w
-
+
1
2
e
ϕ
1
+
τ
2
J
1
/
4
(
1
2
e
2
ϕ
(
1
+
4
τ
2
)
)
k
2
Γ
(
5
4
)
]
.
a
>
0
C
=
exp
[
-
τ
e
2
ϕ
(
1
+
4
τ
2
)
-
1
4
ln
(
1
+
4
τ
2
)
+
a
2
arctan
(
2
τ
)
]
w
w
w
w
w
w
w
w
w
×
[
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
2
1
F
1
(
3
-
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
w
w
w
w
w
w
w
w
w
e
w
=
+
e
-
i
/
2
e
2
ϕ
(
1
+
4
τ
2
)
1
2
e
2
ϕ
(
1
+
4
τ
2
)
k
1
U
(
3
-
i
a
4
,
3
2
;
i
e
2
ϕ
(
1
+
4
τ
2
)
)
]
.Table 4
Group invariant solutions for (4) given D
(
v
)
=
λ
e
2
ϕ.
Symmetry
Reductions
a
∈
ℝ
Group invariant solutions
X
5
+
a
X
2
C
=
exp
[
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
8
λ
ln
(
τ
)
]
G
(
γ
), where γ
=
e
ϕ
τ
4 and G satisfies
λ
2
γ
6
G
′
′
-
1
4
γ
λ
(
1
-
12
γ
λ
4
)
G
′
-
a
8
G
=
0.
a
<
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
ln
(
τ
)
8
λ
)
×
[
k
2
e
-
2
ϕ
4
λ
τ
1
F
1
(
4
+
a
8
λ
,
3
2
;
-
e
-
4
ϕ
16
λ
τ
)
+
k
1
1
F
1
(
a
8
λ
,
1
2
;
-
e
-
4
ϕ
16
λ
τ
)
]
.
a
=
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
)
[
k
2
-
λ
π
k
1
erf
(
e
-
2
ϕ
4
λ
τ
)
]
.
a
>
0
C
=
exp
(
e
-
2
ϕ
4
λ
-
τ
4
λ
+
a
8
λ
ln
(
τ
)
)
W
W
W
W
W
W
W
W
W
W
W
W
W
W
W
W
×
[
k
2
e
-
2
ϕ
4
λ
τ
1
F
1
(
4
-
a
8
λ
,
3
2
;
-
e
-
4
ϕ
16
λ
τ
)
+
k
1
1
F
1
(
-
a
8
λ
,
1
2
;
-
e
-
4
ϕ
16
λ
τ
)
]
.
X
1
+
X
4
C
=
exp
[
8
λ
τ
3
3
+
τ
e
-
2
ϕ
-
τ
2
]
G
(
γ
), where γ
=
4
τ
2
+
e
-
2
ϕ and G satisfies
G
′
′
-
1
2
λ
G
′
-
γ
4
λ
G
=
0.
—
C
=
exp
(
8
λ
τ
3
3
+
τ
e
-
2
ϕ
-
τ
2
)
×
{
e
4
τ
2
+
e
-
2
ϕ
4
λ
[
k
1
Ai
(
1
+
4
(
4
τ
2
+
e
-
2
ϕ
)
(
4
λ
)
4
/
3
)
W
W
W
W
W
W
W
+
k
2
Bi
(
1
+
4
(
4
τ
2
+
e
-
2
ϕ
)
(
4
λ
)
4
/
3
)
]
}
.
X
1
+
a
X
2
C
=
e
a
τ
G
(
γ
), where γ
=
ϕ and G satisfies
G
′
′
+
(
2
+
e
-
2
ϕ
λ
)
G
′
+
a
λ
e
-
4
ϕ
G
=
0.
a
<
0
C
(
τ
,
ϕ
)
=
e
a
τ
[
exp
(
-
1
4
λ
(
-
1
-
1
-
4
a
λ
)
e
-
2
ϕ
)
k
1
+
exp
(
-
1
4
λ
(
-
1
+
1
-
4
a
λ
)
e
-
2
ϕ
)
k
2
].
a
=
0
C
=
[
k
2
+
exp
(
e
-
2
ϕ
2
)
k
1
].
a
>
0
C
=
e
a
τ
exp
[
(
1
4
λ
(
1
-
1
+
4
a
λ
)
e
-
2
ϕ
)
k
1
+
exp
(
1
4
λ
(
1
+
1
+
4
a
λ
)
e
-
2
ϕ
)
k
2
].
X
1
+
X
6
+
a
X
2
Leading to reductions which are unsolvableDefinition 3 (see [24]).
u
=
H
(
x
) is a group-invariant solution of a kth-order PDE corresponding to the appropriate admitted symmetry generator if and only if u
=
H
(
x
) satisfies
(22)
X
(
u
-
H
(
x
)
)
=
0
,
when
u
=
H
(
x
)
.
That is,
(23)
ξ
i
(
x
,
H
(
x
)
)
∂
H
∂
x
i
=
η
(
x
,
H
(
x
)
)
.
The solution of (23) is obtained from the characteristics equation
(24)
d
ξ
1
d
x
1
=
d
ξ
2
d
x
2
=
⋯
=
d
η
u
.
Classifications of the group-invariant solutions by the elements of the optimal system in (20) and (21) are listed in Tables 3 and 4, respectively. Wherever they appear, k
1 and k
2 are arbitrary constants, 1
F
1
(
b
,
c
;
z
) and U
(
b
,
c
;
z
) are the confluent hypergeometric functions, Ai
(
z
) and Bi
(
z
) are the Airy functions, while J
n
(
z
) and Γ
(
z
) represent the Bessel function of  the first kind and Euler gamma function, respectively, and erf
(
z
) represents an error function. A well-documented review of such functions is presented by Abramowitz and Stegun [27].


## 5.3.  Some Physical Examples

### 5.3.1. Given Constant Dispersion Coefficient

Example  a. Suppose a concentration C
0 of a solute is supplied to a single point in an instant of time. We require to determine the subsequent concentration of the pollution at various distances from where it was released (see, e.g., [28]). We would expect concentration to vanish at large distance; that is,
(25)
C
(
ϕ
,
τ
)
⟶
0
,
as
ϕ
⟶
-
∞
.The generatorX
4 in (11) (also, an element of the one-dimensional optimal systems given a
=
0) leads to an invariant solution in functional form given by
(26)
C
=
G
(
γ
)
,
where γ
=
τ
e
ϕ and G satisfies the equation
(27)
2
γ
3
G
′′
+
(
4
γ
2
-
1
)
G
′
=
0
,
and hence
(28)
G
=
c
1
+
c
2
erf
(
1
2
γ
)
.
Here erf
(
z
)
=
(
2
/
π
)
∫
0
z
‍
e
-
ν
2
d
ν is the error function [27]. In terms of the original variable and subject to the boundary conditions, we obtain
(29)
C
=
C
0
erfc
(
e
-
ϕ
2
τ
)
.
Here, the erfc(·) is the complement error function defined by (
1
-
erf
(
·
)
). Solution (29) is depicted in Figure 1.Figure 1
Contaminant concentration profile along the radius. HereC
0
=
5.Total flux acrossR
=
R
a is given by
(30)
C
+
1
R
∂
C
∂
R
|
R
=
R
a
=
C
0
erfc
(
R
a
2
τ
)
-
C
0
e
-
R
a
2
/
4
τ
π
τ
R
a
.
The contaminant flux (30) is depicted in Figure 2. Total flux across R
a
=
4 increases and flattens at large time.Figure 2
Solute flux across a fixed radius. HereR
a
=
4 and C
0
=
10.Example  b. The X
6-invariant solution is given in functional form as
(31)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
G
(
γ
)
,
where
(32)
γ
=
τ
e
ϕ
,
G
satisfies
the
ODE
γ
G
′′
+
G
′
=
0
.
We impose the boundary conditions
(33)
C
⟶
0
,
ϕ
⟶
-
∞
,
C
=
ω
(
τ
)
,
ϕ
⟶
∞
.
Infinite concentration at the origin implies that there is a high supply of contaminants at this point. Furthermore, contaminant concentration vanishes when time evolves. In terms of the original variables we obtain the exact (group-invariant) solution given by
(34)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
.Solution (34) is depicted in Figures 3 and 4. In Figure 3, a sharp peak of concentration is observed shortly after τ
=
0 and decreases at later stage. This may be interpreted as an injection of contaminants at a single point; that is, the concentration at a single point increases but due to diffusion at larger time it smoothes out. Note that here we have restricted our analysis using symmetry generator X
6. This symmetry generator leads to simpler and realistic exact solution. In Figure 4, we observe that concentration at the origin decreases with time. Furthermore, this concentration vanishes at large distances.Figure 3
Contaminant concentration profile as time evolves.Figure 4
Contaminant concentration profile at various fixed times.


### 5.3.2. Given Velocity-Dependent Dispersion Coefficient

Example
a (steady-state solution). The time translation leads to the analysis of the steady-state contaminant transport. Steady-state solutions may be constructed and subjected to the following imposed boundary conditions:
(35)
C
=
C
0
,
R
=
0
,
(36)
C
+
1
R
d
C
d
R
=
0
,
R
=
R
a
.
The boundary condition (35) implies that pollutants are supplied at the origin, and boundary condition (36) correspond to the assumption that pollutants are not carried though at some distance R
a, rather it accumulates here. We obtain the exact solution
(37)
C
=
C
0
{
(
1
-
1
Δ
)
+
1
Δ
exp
(
R
p
λ
p
)
}
,
where Δ is given by
(38)
Δ
=
1
-
(
1
+
R
a
p
-
2
λ
)
exp
(
R
a
p
λ
p
)
.
The solution (37) is depicted in Figures 5, 6, and 7. We observe in Figure 5 that concentration starts decreasing and converges to some value at large distance for p
=
2 than for lower values of p, whereas λ has an opposite effect as shown in Figure 6.Figure 5
Steady contaminant profile for concentration given in (37). Here λ
=
1.Figure 6
Steady contaminant profile for concentration given in (37). Here R
=
R
a
=
2.Figure 7
Effects ofλ and p on concentration profile for solution given in (37). Here R
=
R
a
=
2.Example
b (transient-state solution). It is quite difficult to construct exact solutions for transient contaminant transport subject to these boundary conditions (35) and (36). However, if one assumes that at an initial time, say τ
=
1, the concentration at the point source is given by a constant and that this concentration vanishes at large distances and prolonged periods, then using the symmetry combination X
1
+
a
X
2 from Table 4, the group invariant solution is given by
(39)
C
=
C
0
e
a
τ
exp
{
(
1
-
1
-
4
a
λ
4
λ
)
R
2
}
,
∀
a
<
0
,
λ
>
0
.
Solution (39) is depicted in Figures 8 and 9.Figure 8
Concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and λ
=
1Figure 9
Effects ofλ on the concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and τ
=
2.In his work, Philip [6] considered the instantaneous point source for contaminant dispersion during radial water flow in porous media. Exact solutions were constructed for the two- and three-dimensional models with dispersion coefficient depending on Péclet number. Here, we consider models in stream functions coordinates. Exact close-form (similarity) solutions are constructed using the elements of the one-dimensional optimal systems. These new solutions may be viewed as representing the continual supply of contaminant at a point (source) which are dispersed radially.


## 5.3.1. Given Constant Dispersion Coefficient

Example  a. Suppose a concentration C
0 of a solute is supplied to a single point in an instant of time. We require to determine the subsequent concentration of the pollution at various distances from where it was released (see, e.g., [28]). We would expect concentration to vanish at large distance; that is,
(25)
C
(
ϕ
,
τ
)
⟶
0
,
as
ϕ
⟶
-
∞
.The generatorX
4 in (11) (also, an element of the one-dimensional optimal systems given a
=
0) leads to an invariant solution in functional form given by
(26)
C
=
G
(
γ
)
,
where γ
=
τ
e
ϕ and G satisfies the equation
(27)
2
γ
3
G
′′
+
(
4
γ
2
-
1
)
G
′
=
0
,
and hence
(28)
G
=
c
1
+
c
2
erf
(
1
2
γ
)
.
Here erf
(
z
)
=
(
2
/
π
)
∫
0
z
‍
e
-
ν
2
d
ν is the error function [27]. In terms of the original variable and subject to the boundary conditions, we obtain
(29)
C
=
C
0
erfc
(
e
-
ϕ
2
τ
)
.
Here, the erfc(·) is the complement error function defined by (
1
-
erf
(
·
)
). Solution (29) is depicted in Figure 1.Figure 1
Contaminant concentration profile along the radius. HereC
0
=
5.Total flux acrossR
=
R
a is given by
(30)
C
+
1
R
∂
C
∂
R
|
R
=
R
a
=
C
0
erfc
(
R
a
2
τ
)
-
C
0
e
-
R
a
2
/
4
τ
π
τ
R
a
.
The contaminant flux (30) is depicted in Figure 2. Total flux across R
a
=
4 increases and flattens at large time.Figure 2
Solute flux across a fixed radius. HereR
a
=
4 and C
0
=
10.Example  b. The X
6-invariant solution is given in functional form as
(31)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
G
(
γ
)
,
where
(32)
γ
=
τ
e
ϕ
,
G
satisfies
the
ODE
γ
G
′′
+
G
′
=
0
.
We impose the boundary conditions
(33)
C
⟶
0
,
ϕ
⟶
-
∞
,
C
=
ω
(
τ
)
,
ϕ
⟶
∞
.
Infinite concentration at the origin implies that there is a high supply of contaminants at this point. Furthermore, contaminant concentration vanishes when time evolves. In terms of the original variables we obtain the exact (group-invariant) solution given by
(34)
C
=
1
τ
exp
(
-
e
-
2
ϕ
4
τ
)
.Solution (34) is depicted in Figures 3 and 4. In Figure 3, a sharp peak of concentration is observed shortly after τ
=
0 and decreases at later stage. This may be interpreted as an injection of contaminants at a single point; that is, the concentration at a single point increases but due to diffusion at larger time it smoothes out. Note that here we have restricted our analysis using symmetry generator X
6. This symmetry generator leads to simpler and realistic exact solution. In Figure 4, we observe that concentration at the origin decreases with time. Furthermore, this concentration vanishes at large distances.Figure 3
Contaminant concentration profile as time evolves.Figure 4
Contaminant concentration profile at various fixed times.


## 5.3.2. Given Velocity-Dependent Dispersion Coefficient

Example
a (steady-state solution). The time translation leads to the analysis of the steady-state contaminant transport. Steady-state solutions may be constructed and subjected to the following imposed boundary conditions:
(35)
C
=
C
0
,
R
=
0
,
(36)
C
+
1
R
d
C
d
R
=
0
,
R
=
R
a
.
The boundary condition (35) implies that pollutants are supplied at the origin, and boundary condition (36) correspond to the assumption that pollutants are not carried though at some distance R
a, rather it accumulates here. We obtain the exact solution
(37)
C
=
C
0
{
(
1
-
1
Δ
)
+
1
Δ
exp
(
R
p
λ
p
)
}
,
where Δ is given by
(38)
Δ
=
1
-
(
1
+
R
a
p
-
2
λ
)
exp
(
R
a
p
λ
p
)
.
The solution (37) is depicted in Figures 5, 6, and 7. We observe in Figure 5 that concentration starts decreasing and converges to some value at large distance for p
=
2 than for lower values of p, whereas λ has an opposite effect as shown in Figure 6.Figure 5
Steady contaminant profile for concentration given in (37). Here λ
=
1.Figure 6
Steady contaminant profile for concentration given in (37). Here R
=
R
a
=
2.Figure 7
Effects ofλ and p on concentration profile for solution given in (37). Here R
=
R
a
=
2.Example
b (transient-state solution). It is quite difficult to construct exact solutions for transient contaminant transport subject to these boundary conditions (35) and (36). However, if one assumes that at an initial time, say τ
=
1, the concentration at the point source is given by a constant and that this concentration vanishes at large distances and prolonged periods, then using the symmetry combination X
1
+
a
X
2 from Table 4, the group invariant solution is given by
(39)
C
=
C
0
e
a
τ
exp
{
(
1
-
1
-
4
a
λ
4
λ
)
R
2
}
,
∀
a
<
0
,
λ
>
0
.
Solution (39) is depicted in Figures 8 and 9.Figure 8
Concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and λ
=
1Figure 9
Effects ofλ on the concentration profile for solution given in (39). Here R
=
R
a
=
2
,
C
0
=
1, and τ
=
2.In his work, Philip [6] considered the instantaneous point source for contaminant dispersion during radial water flow in porous media. Exact solutions were constructed for the two- and three-dimensional models with dispersion coefficient depending on Péclet number. Here, we consider models in stream functions coordinates. Exact close-form (similarity) solutions are constructed using the elements of the one-dimensional optimal systems. These new solutions may be viewed as representing the continual supply of contaminant at a point (source) which are dispersed radially.


## 6. Some Concluding Remarks

In this paper, we have focused only on the two-dimensional solute concentration field within water from a single injection well. The considered problem is a significant improvement in the study of solute transport under radial water background since we analyze the convection-dispersion equation in stream functions. We have observed that extra Lie point symmetries are admitted when the dispersion coefficient is a constant or when it is given as a power law function of velocity, with exponent being given by two. We have classified the group invariant solutions by the elements of the optimal systems. In fact, new exact solutions are constructed. The symmetry (invariant) solution is obtainable when dispersion coefficient is a constant or is given by Taylor’s theory of mixing in soils.



---

*Source: 138289-2014-02-02.xml*