EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 11, No. 1, 2018, 110-137
ISSN 1307-5543 – www.ejpam.com
Published by New York Business Global

A new mathematical model of heat equations and its
application on the agriculture soil

Ali Allahem1,2, Salah Boulaaras2,3,∗, Khaled Zennir2,4 and Mohamed Haiour5

1 Department of Mathematics, College of Sciences, Qassim University,
Kingdom of Saudi Arabia
2 Department of Mathematics, College of Sciences and Arts, Al-Rass, Qassim University,
Kingdom of Saudi Arabia
3 Laboratory of Fundamental and Applied Mathematics of Oran (LMFAO), University of
Oran 1, Ahmed Benbella, Algeria.
4 LAMAHIS Laboratory, Department of Mathematics, University August 20, 1955- Skikda,
21000, Algeria
5 Department of Mathematics, Faculty of Sciences, Annaba University, Algeria

Abstract. The goal of this paper is to provide a mathematical model which based on numerical
simulation applied to heat equations and comparing the obtained results with virtual temperatures
results for different depths of the homogeneous agricultural soil. The model was treated through
this research that enable to know the degree of soil temperature homogeneous without recourse
to tools of measuring physical expensive. The importance of this lies in the agricultural seed
germination at temperatures suitable.

2010 Mathematics Subject Classifications: AMS 65N06, 65N12, 65F05

Key Words and Phrases: Numerical model, Finite elements, Temperature fields, Meteorology,Heat-
Soil, Agricultural

1. Introduction

In the recent years, the alternative methods are used to speed in solving partial differ-
ential equations that have been explained the natural phenomena, one of these methods
is the finite element method. In this project, we use the finite element method which
belong to the family in general methods of Galerkin. These are used instead of the style
differences limited (considered simplistic) to resolve the horizontal and vertical fields in
numerical model. Galerkin methods can be used in solving systems of partial differential
equations fixed and the problems of border free evolutionary (see reference [1], [2], [5],

∗Corresponding author.

Email addresses: a.allahem@qu.edu.sa (A. Allahem), saleh boulaares@yahoo.fr or

S.Boularas@qu.edu.sa (S. Boulaaras), khaledzennir2@yahoo.com or k.zennir@qu.edu.sa (Kh.
Zennir), haiourm@yahoo.fr (M. Haiour)

http://www.ejpam.com 110 c© 2018 EJPAM All rights reserved.



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 111

[6], [7], [12], [13] and [7]) do not use this method directly in the field of values at grid
points after the partition of scale, so as to improve the accuracy of ordinary differential
equations systems solutions. There are two methods in this process: the specific elements
of the functions of the zero-way, and the way spectral. Finite element method is one of
the ways of applied Mathematics. And we put it in space, using the principles inherited
from the transformative formulation or weak wording, separate mathematical algorithm
to find approximate solution to the problem of free boundary (see [3], [4] and [14]). It is
said Dirichlet conditions or Neumann data.

Finite element method differs from the spectral method because it is not exhaustive,
but rather to determine local values. However, it is a rough distinct function defined on
the whole region and not just separate points (see [8] and [10]).

The techniques used can solve with numerical methods, the equations describing the
behavior of the atmosphere [16] and [13], that is to say to determine future values of its
parameters characteristics starting from known initial values through meteorological ob-
servations [9], [11], [18] and [15]. The use of numerical computation is required to solve
these systems of nonlinear equations whose solutions cannot be determined analytically
in the general case. The construction of a numerical model of the atmosphere consists of
two distinct stages: the first is to establish a system of equations governing the behavior
of continuous atmosphere that called the second scan is to replace the equations on con-
tinuous variables by equations involving discrete variable whose solutions are obtained by
means of a suitable algorithm. The results of a numerical weather prediction, solutions of
discretized equations of dynamic meteorology, therefore depend on the adopted scanning.
In the late ’40s methods of grid points were used to model the universally large-scale at-
mospheric flow. During recent years ([12], [19] and [20]), alternative methods have been
used, one of these methods is the finite element method. We present here the main tools
for implementations of the finite element method which belongs to the more general family
of Galerkin methods. These are commonly used instead of the finite difference method to
treat the horizontal and vertical fields in weather prediction models. Galerkin methods,
which can solve numerically systems of partial differential equations, do not directly use
the field values at the points of a grid, but use series expansions of functions suitably cho-
sen so to be reduced to the solution of a system of ordinary differential equations. There
are two types of methods within this process: the finite element method for which the
functions are zero, except for a small area or they are equal to the low-order polynomials,
and the spectral method in which the functions are the functions of a spatial operator
defined on the whole area of the works.

The finite element method is one of the tools of applied mathematics. It is puts in
place, using principles inherited from the variational formulation or weak formulation, a
discrete mathematical algorithm for finding an approximate solution of a partial differ-
ential equation on a compact domain with boundary conditions and / or in the interior
of the compact. It speaks Dirichlet conditions (values at the edges) or Neumann (gra-
dients at the edges). The finite element method is different than the spectral because
it is not comprehensive, but rather determined by local values. However, it is distinct
approximations gridded because the function is defined over the entire region and not just



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 112

the discrete points. As in the case of other numerical methods, questions arise as to the
discretization: Existence and uniqueness of solutions, stability, convergence and of course,
the error measure between a discrete and a unique solution of the initial problem.

Our goal in this paper is to highlight the application of the method to an equation of
parabolic type. In the practical part, we applied this method to the diffusion equation
for virtual temperatures results in the ground in Qassim province in Kingdom of Saudi
Arabia.

This paper is organized as follows: In the first part, we introduce some necessary
notations of the problem , then we study the numerical analysis side of the heat equations
and we obtain the weak formulation of our model using the theta time scheme combined
with a finite element method. Then the stability analysis of the proposed schemes and
the theoretical side of the full discrete problem is established. In the second part, we
present a one-dimensional numerical prediction model that predicts the temperature field
for different depths of soil homogeneous and for different times and different days of a given
month. By using the finite element approximation of the Lagrange polynomial of order
1, determining the change in temperature in the soil over time by the knowledge of the
parameters and characteristics of the soil. During a period of 24 hours, the temperature
surface is substantially the response to radiation received diurnal. Exchanges heat does
make the firsts in soil layers that can return to the surface if necessary, part of the energy.
By cons, in a day, the heat storage in the soil (thermal inertia) concerns only the first layer
of soil not exceeding 60 cm. I.e. the temperature field in the surface layer of the soil, which
extends to a depth of 60 cm, is still sensitive to climatic variations recorded on the surface.
The objective of this work is precisely to propose a prediction model based on dimensional
numerical finite element approximation for the evolution of the temperature field in a
homogeneous soil. Finally, we perform virtual measurements of the soil temperature and
in the soil, that is to say, we use virtual data to prove the accuracy of the numerical
method by using the error estimate.

2. Temperature diffusion equation

The objective of this work is precisely to propose a one-dimensional numerical predic-
tion model based on the finite element approximation for the evolution of the temperature
field in a homogeneous soil. For the case of a homogeneous and isotropic soil, the ther-
mal balance on a volume element, gives rise according to the depth z, to the equation of
propagation of the heat with internal production of energy f which is written:

ut −Duzz = f in [a, b]× [0, tf ] ,

u(t, a) = us = α, u(t, b) = β,

u(0, z) = u0(z),

(1)



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 113

where

D =
ρ.λ

cp
, (2)

and
ρ is the density, λ is the thermal conductivity and cp is the specific heat.
According to ([17]), the thermal diffusively in the soil is 10−6m2/s.

2.1. Discrete problem

To solve this problem (1), one can proceed in two stages. In a first step, we discretizes
in space by finite elements, that is to say that we approach the space H1 by a finite
dimension space V h

1 defined as:

V h
1 =


vh ∈ C2 (Ω) ∩H1 (Ω) ,

vh |I=[a,b]∈ P1 in Ω

and u (a) = α, u (b) = β.


2.1.1. Spatial discretization

We then define the semi-discretized problem, multiplying equation (1) by a test function v ∈
H1 ([a, b]) and we integrate on [a, b] = Ω, we get∫

Ω

utvdz −D
∫
Ω

uzzvdz =

∫
Ω

fvdz. (3)

To introduce an approximation uh ∈ V h
1 ⊂ H1(Ω), we will introduce a mesh of Ω

composed of N + 1 points and if the mesh is uniform, the step of spatial discretization is

h =
b− a
N

.

But in our case the mesh is not uniform, we take:

h = maxhj .

We introduce the following piecewise linear approximation:

V h
1 =



vh ∈ C2(Ω) ∩H1(Ω) :

vh

Ij

= p1 ∈ P1

and vh(a) = α, vh(b) = 0,





A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 114

where

vh =

N−1∑
j=1

vj(z)ϕj(z), u
h =

N∑
i=0

ui(t)ϕi(z)

with {ϕi}i=1,...,N form a vector basis uh and {ϕj}j=1,...,N for vh.

Problem (3) is equivalent to∫
Ω

ut.v

︸ ︷︷ ︸
dz

A

+D(

∫
Ω

uzz.vz

︸ ︷︷ ︸
dz

B

−
∫
Γ

uzz.
−→η .v.dσ

︸ ︷︷ ︸
)

C

=
∫
Ω

fvdz.

(4)

Terms B and C are obtained using Green’s formula as:

∫
Ω

uzzvz.vdz =

∫
Ω

uzvzdz

︸ ︷︷ ︸
B

−
∫
Γ

uz.η.v.dσ

︸ ︷︷ ︸
)

C

=

∫
Ω


∂
∂z

N−1∑
i=1

ui(t)ϕi(z)

∂
∂Z

N−1∑
j=1

vj(z)ϕj(z)

 dz

︸ ︷︷ ︸
B

−
∫
Γ

uz.
−→η .v.dσ

︸ ︷︷ ︸
)

C

.

To eliminate the term C, we can change the summation of u from 1, ..., N − 1 to 0, N
(that is to say, to consider the boundaries α and β), thus

∫
Ω

uzzvdz =

∫
Ω


∂
∂z

N∑
i=0
ui(t)ϕi(z)

∂
∂z

N−1∑
j=1

vj(z)ϕj(z)

 dz



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 115

with
N∑
i=0
ui(t)ϕi(z) = αϕ0(z)

+
N−1∑
i=1

ui(t)ϕi(z) + βϕN (z).

and

A =
∂

∂t
(uh, vh), B = a(uh, vh).

So equation (4) is written in the following variational form:

∂

∂t
(uh, vh) +D.a(uh, vh) = (fh, vh). (5)

2.1.2. Time discretization

In a second step, we discretize (5) with respect to time by using theta scheme. We are
therefore looking for a series of elements uhn ∈ V h which will have to approach uh(tn) , tn =
n∆t and which are defined by the following formula: find uhn ∈ V h , for all vh ∈ V h and
θ ∈ [0, 1] 

(
uhn+1 − uhn

∆t
, vh)

+Da
(
θuhn+1 + (1− θ)uhn , vh

)
= (θfn+1

h + (1− θ) fnh , vh).

(6)

Equation (6) equivalent to

(
uhn+1 , v

h
)

+Dθ∆t.a (uhn+1 , v
h)

+(θfn+1
h , vh) =

(
uhn , v

h
)

−D (1− θ) ∆t
(
a
(
uhn , v

h
)

+ (fnh , v
h)
)
.

(7)

By recurrence, given uh0 , we have the existence and uniqueness using Lax-Milgram
theorem.

3. Stability analysis of the theta scheme

To study the stability, we use the energy increment technique, which allow to conclude
without hypothesis on V h when the scheme is unconditionally stable. For the case of
a homogeneous and isotropic soil, the thermal balance on a volume element, gives rise
according to the depth z, to the equation of propagation of the heat without internal



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 116

production of energy (the right hand side of the problem equal zero), thus we can write
equation (7) as: (

uhn+1 − uhn
∆t

, vh

)
= D (θ − 1) a

(
uhn , v

h
)

−Dθa (uhn+1 , v
h)

= −D
[
a
(
(1− θ)uhn , vh

)
+ a (θuhn+1 , v

h)
]

= −D
[
a
(
(1− θ)uhn + θuhn+1, v

h
)]

= −D
[
a
(
θ(uhn+1 − uhn) + uhn, v

h
)]
.

Then ∣∣∣∣∣∣∣∣
(

uh
n+1−uh

n

∆t , vh
)

+D
[
a
(
θ(uhn+1 − uhn) + uhn, v

h
)]
∣∣∣∣∣∣∣∣ = 0. (8)

Putting the test function

vh = θ(uhn+1 − uhn) + uhn

= ũh.

Equation (8) becomes (
uhn+1 − uhn, θ(uhn+1 − uhn) + uhn

)
︸ ︷︷ ︸

A

+D∆ta
(
ũh , ũh

)
︸ ︷︷ ︸

B

= 0.

(9)

The term A of the equation equivalent to

A =
(
uhn+1 − uhn , θ(uhn+1 − uhn)

)
+
(
uhn+1 − uhn , uhn

)
= θ

(
uhn+1 − uhn , uhn+1 − uhn

)
−1

2

(
uhn+1 − uhn , −2uhn

)
.



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 117

Thus
A = θ

(
uhn+1 − uhn , uhn+1 − uhn

)
−1

2

(
uhn+1 − uhn , −uhn − uhn + uhn+1 − uhn+1

)
= θ

(
uhn+1 − uhn , uhn+1 − uhn

)
−1

2

(
uhn+1 − uhn , uhn+1 − uhn

)
= θ

(
uhn+1 − uhn , uhn+1 − uhn

)
−1

2

(
uhn+1 − uhn , uhn+1 − uhn

)
,

so
A = θ

(
uhn+1 − uhn , uhn+1 − uhn

)
−1

2

(
uhn+1 − uhn , uhn+1 − uhn

)
= (θ − 1

2)
(
uhn+1 − uhn , uhn+1 − uhn

)
+1

2

(
uhn+1 − uhn , uhn+1 + uhn

)
.

Therefore,

A = (θ − 1
2)
∥∥(uhn+1 − uhn)

∥∥2

L2(Ω)

+
1

2
(
∫

(uhn+1)2 − (uhn)2)

= (θ − 1
2)
∥∥(uhn+1 − uhn)

∥∥2

L2(Ω)

+1
2

∥∥uhn+1

∥∥2

L2(Ω)
− 1

2

∥∥uhn∥∥2

L2(Ω)
.

By adding the term A into equation (9)

(θ − 1
2)
∥∥(uhn+1 − uhn)

∥∥2

L2(Ω)

+1
2

∥∥uhn+1

∥∥2

L2(Ω)

−1
2

∥∥uhn∥∥2

L2(Ω)
+D∆ta (ũ , ũ ) = 0.

This means



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 118

1
2

∥∥uhn+1

∥∥2

L2(Ω)
− 1

2

∥∥uhn∥∥2

L2(Ω)

+∆tDa (ũ , ũ )

= −(θ − 1
2)
∥∥(uhn+1 − uhn)

∥∥2

L2(Ω)
.

(10)

If θ ≥ 1
2 , we have

1
2

∥∥uhn+1

∥∥2

L2(Ω)
− 1

2

∥∥uhn∥∥2

L2(Ω)

+∆tDa
(
ũh , ũh

)
≤ 0,

thus ∥∥uhn+1

∥∥2

L2(Ω)
+ 2∆tDa

(
ũh, ũh

)
≤
∥∥uhn∥∥2

L2(Ω)
,

then ∥∥uhn+1

∥∥2

L2(Ω)
+ 2∆tD

∥∥∥∂ũh

∂z

∥∥∥2

L2(Ω)

≤
∥∥uhn∥∥2

L2(Ω)
.

By summation with respect to n, we find∥∥uhn∥∥2

L2(Ω)
+ 2∆tD

n∑
n=1

∥∥∥∂ũh

∂z

∥∥∥2

L2(Ω)

≤
∥∥∥uh0∥∥∥2

L2(Ω)︸ ︷︷ ︸
known

.

It is deduced that for θ ≥ 1

2
the scheme is unconditionally stable.

3.1. Matrix formulation of the problem

We can reformulate equation (7)(
uhn+1 , v

h
)

+Dθ∆ta (uhn+1 , v
h)

=
(
uhn , v

h
)
−D (1− θ) ∆ta

(
uhn , v

h
)
,

ω = θ∆tD.



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 119

We have again

∫
Ω

(
N−1∑
i=1

un+1
i (t)ϕi(z)

n−1∑
j=1

vjϕj(z)

)
dz

+ω
∫

Ω
∂
∂z

(
N−1∑
i=1

un+1
i (t)ϕi(z)

)
∂
∂z

(
n−1∑
j=1

vjϕj(z)

)
dz

=
∫

Ω

(
N−1∑
i=1

uni (t)ϕi(z)
n−1∑
j=1

vjϕj(z)

)
dz

−(1− θ)∆tD
∫

Ω
∂
∂z

(
N∑
i=0
uni (t)ϕi(z)

)
∂
∂z

(
n−1∑
j=1

vjϕj(z)

)
dz.

(11)

Equation (11) is equal to

N−1∑
i=1

un+1
i (t)

∫
sup pϕi∩sup pϕj

ϕi(z)ϕj(z)dz

+ω
N−1∑
i=1

un+1
i (t)

∫
sup pϕi∩sup pϕj

∂
∂zϕi(z)

∂
∂zϕj(z)dz

=
N−1∑
i=1

uni (t)
∫

sup pϕi∩sup pϕj

ϕi(z)ϕj(z)dz

−(1− θ)∆tD
N∑
i=0
uni (t)

∫
sup pϕi∩sup pϕj

∂
∂zϕi(z)

∂
∂zϕj(z)dz,

thus



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 120

N−1∑
i=1

un+1
i (t)


∫

suppϕi∩suppϕj

ϕi(z) · ϕj(z)dZ

+ω
∫

sup pϕi∩sup pϕj

∂
∂Zϕi(z)

∂
∂Zϕj(z)dz


= −α(1− θ)∆tD

∫
sup pϕi∩sup pϕj

∂
∂zϕ0(z) ∂

∂zϕ1(z)dz

−β(1− θ)∆tD
∫

sup pϕi∩sup pϕj

∂
∂zϕN (z) ∂

∂zϕN−1(z)dz

+
N−1∑
i=1

uni (t)


∫

sup pϕi∩sup pϕj

ϕi(z) · ϕj(z)dx

−(1− θ)∆tD
∫

sup pϕi∩sup pϕj

∂
∂zϕi(x) ∂

∂zϕj(z)dz.



(12)

Then

N−1∑
i=1

un+1
i (t)


xi+1∫
xi−1

ϕi(z) · ϕj(z)dx

+θ∆tD
xi+1∫
xi−1

∂
∂zϕi(z) · ∂

∂zϕj(z)dz


= −α(1− θ)∆tD ·

∫
sup pϕi∩sup pϕj

∂
∂zϕ0(z) · ∂

∂zϕ1(z)dz

−β(1− θ)∆tD ·
∫

sup pϕi∩sup pϕj

∂
∂zϕN (z) ∂

∂zϕN−1(z)dz

+
N−1∑
i=1

uni (t)


xi+1∫
xi−1

ϕi(z)ϕj(z)dz

−(1− θ)∆tD
xi+1∫
xi−1

∂
∂zϕi(z)

∂
∂xϕj(z) dz



(13)

with basic functions

ϕi(z) =


z−zi−1

h , z ∈ [zi−1, zi] ,

Zi+1−Z
h , z ∈ [zi, zi+1]

and



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 121

ϕi−1(z) =


z−zi−2

h , z ∈ [zi−2, zi−1] ,

zi−z
h , z ∈ [zi−1, zi]

with
h = (zi+1 − zi) = (zi − zi−1) .

To make the calculations easier, we make a change change such as{
z = as+ b

z ∈ [z, zi+1]
Linear withdrawal−→ s ∈ [0, 1] ,

thus, it can be easily found:

ϕ0(s) =


(1 + s) , s ∈ [−1, 0] ,

(1− s) , s ∈ [0, 1] ,

ϕ−1(s) =


(1 + s) , s ∈ [−2,−1] ,

−s , s ∈ [−1, 0]

and

ϕ1(s) =


(1− s) , s ∈ [0, 1] ,

−s, s ∈ [1, 2] .

Equation (13) becomes



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 122

N−1∑
i=1

un+1
i (t)



zi+1∫
zi−1

ϕi(z) · ϕj(z)dz

︸ ︷︷ ︸
M

+θ∆tD

zi+1∫
zi−1

∂

∂z
ϕi(z) ·

∂

∂z
ϕj(z) dz

︸ ︷︷ ︸
A


= −α(1− θ)∆tD

∫
sup pϕi∩sup pϕj

∂
∂zϕ0(z) · ∂

∂zϕ1(z)dz

−β(1− θ)∆tD
∫

sup pϕi∩sup pϕj

∂
∂zϕN (z) · ∂

∂zϕN−1(z)dz

+
N−1∑
i=1

uni (t)


zi+1∫
zi−1

ϕi(z) · ϕj(z)dz

−(1− θ)∆tD
zi+1∫
zi−1

∂
∂Zϕi(z) · ∂

∂zϕj(z)dz.


Therefore, we get the following linear system

(M + θ∆tDA)un+1 = (M − (1− θ)∆tDA)un

−(1− θ)∆tDB,
(14)

where M and A are symmetric tridiagonal matrix, their dimension is (N − 1)×(N − 1)
and

M = m (ϕi, ϕj) =

zi+1∫
zi−1

ϕi(z) · ϕj(z)dz,

A = a (ϕi, ϕj) =

zi+1∫
zi−1

ϕi(z) · ϕj(z)dz

with

B= (αa(ϕ0, ϕ1), ..., βa(ϕN , ϕN−1))t .



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 123

In addition, we have un+1 =
N−1∑
i=1

un+1
i (t) and un =

N−1∑
i=1

Tn
i (t). Now, it is easy to

calculate M(ϕi, ϕj), A(ϕi, ϕj) as following:

M(ϕi, ϕj) =


2h
3 , i = j,

h
6 , |i− j| = 1

and

A(ϕi, ϕj) =


2
h , i = j,

− 1
h , |i− j| = 1.

4. Application

4.1. Error evaluation

To put our work into practice, measurements of the temperature on the ground and in
the soil are made, that is, virtual data was used to highlight the accuracy of the method
using the error assessment according to weather temperatures data at Qassim province,
Qassim province is one of the thirteen administrative regions of Kingdom of Saudi Arabia.
Located at the heart of the country, and almost in the center of the Arabian Peninsula, it
has a population of 1, 370, 727 and an area of 58, 046 km2. It is known to be the ”alimental
basket” of the country, for its agricultural assets.

4.2. Obtaining initial data

1-The data are considered according to the weather data: the temperature at agricul-
tural ground level, the maximum and minimum temperature on the ground, the tempera-
tures of four levels in the soil (10, 20, 30 and 50cm). We consider the data of March 2017
whose selected days are from 7 to 10.

2-As with all numerical forecast models, data at the boundaries of the domain is
needed. The boundary α is the ground temperature that is calculated every 10 minutes
(using a linear interpolation between the minimum and the maximum of a given day, while
considering the maximum and minimum tim) and the boundary β is considered according
to a previous study in [17], during a day, the storage of heat in the soil concerns only the
first layer of soil not exceeding 60cm. The results obtained show that the temperature
field in the upper soil layer, which extends over a depth of 60cm, remains sensitive to
climatic variations recorded at the surface.

Since the data are not available at 60cm in the soil, the temperature was estimated at
60cm according to the variation of the thermal gradient observed in the soil. Moreover,
after this temperature has been fixed during the whole experiment (this hypothesis is
proved by a previous study, the results are in Fig 1).



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 124

Note that the temperature at 60 cm remains virtuality constant during the 24 hours.
3-Taking ∆t = 10 min and h = maxhj = 20cm .
4-We can turn the model every 10 minutes until while displaying every 12 hours the

results obtained and compared with the actual data to evaluate the mean squared error.

EQM =
n=N∑
n=1

√
1

N
(Texpected − Treal)2

6- The calculations can be done with two methods: the first with initialization, where
we can turn the model and initialize the data at 6h, 12h and 18h and the second is without
initialization: the model is rotated during the entire deadline without resetting the data.

4.3. Virtual experience

The code used in this work is based on the finite element technique to predict the
ground temperature at different levels and for different integration deadlines. For that,
we realized two simulations. The first is to start the calculation without initialization and
retrieve the forecasts at different times, while for the second, the calculation is done with
initialization every twelve hours. We compared the distribution of observed and simulated
temperatures. We have also calculated for each deadline the quadratic error and then the
mean squared error for the whole simulation and the results of the two simulations are as
follows:

4.3.1. Without initialization

According to the results of Fig 2, we note that the model is in good agreement with
the observations for shallow depths for which it gives fairly virtual temperatures, but
differences despite which are weak begin to appear when the depth increases, it is reported
that the model returns exactly the same temperature for the depth of 25 cm.

According to the results of Fig 3, we note that the model always starts with a very
good behavior for the shallow depths and is able to calculate the same temperature as that
measured for an interval between 25 and 30 cm. According to the results in the figures 2
until 8, we note for deadlines that are too high, namely 36 hours, 48 hours and 60 hours,
we notice that there is not a very good agreement between the expected temperatures
and those observed for the different depths. The figures show that the model slightly
overestimates the temperature for shallow depths and underestimates them for significant
depths. It arrives at a good forecast for an average depth between 25 and 30 cm, but
remains close to observations and for 72 hours, 84 hours and 96 hours, the comparison
between the results of the simulation and the observed data shows that the temperatures
are not very well reproduced but remains close to the measurements in the limit where
the difference between the two temperatures do not exceed 0.6 ◦ C.

According to the results of Fig 9, from a statistical point of view, the calculation of
the quadratic error allows a more detailed view of the behavior of the model. It is clear
that the smallest errors were observed for the lowest deadlines this is logical as the model



A. Allahem, S.Boulaaras, K. Zennir and M. Haiour / Eur. J. Pure Appl. Math, 11 (1) (2018), 110-137 125

starts to move away from the value when the maturities increase, but for our case we
notice that the model returns at a fairly good estimate for the deadlines of 84, 96 and 108
hours. We note that the lowest error (0.19) corresponds to the best behavior of the model
at the expiry of 24 hours.

4.3.2. With initialization

According the results of the 12-hours (see Fig 12) simulation with initialization do not
show significant differences between the two simulations, where we find the same behavior
of the model with a slight over estimation for shallow depths and a slight underestimate
for deeper depths. This can be verified by the value of the squared error which gives
exactly the same value for both simulations.

The figures results for the 24 h and 36 h (see Fig 13 and Fig 14) maturities allow the
results of the simulation to be compared to the measurements taken. If the model starts
with a very good behavior for the 24 h deadline, it moves away slightly from the obser-
vations thereafter. A distance materialized by a slight overestimation of the temperature.
This does not say that the model diverges completely because the differences are quite
small. The opposite is observed for the 36 hour deadline for which the model starts an
over-estimation and ends with a good behavior. For the deadlines 96 and 108 h (see Fig
15 up to Fig 20), we notice that the model is in good agreement with the observations for
the important depths for which it gives quite virtual temperatures, but differences despite
which are weak not exceeding 0.6 cm, begin to appear when the depth decreases for the
96 hour deadline. An unverified aspect for 108 hours of expiry where the model has an
underestimation of the temperature for different depths. It is reported that the statistical
results show a positive evolution between the two simulations by the comparison of the
mean squared error which goes from 0.46 for the simulation without initialization to 0.34
for the simulation with initialization. However, Note that there is not a very good agree-
ment between the expected temperatures and those observed for different depths. The
figure shows that the model slightly overestimates the temperature for the weak and the
deep depths and arrives at a good prediction for a depth of 10 cm, but remains close to
observations. However, the results of the simulation do not indicate a fairly clear behavior
of the model with respect to deadlines and depths. In some cases, there is an overes-
timation. in other cases an underestimation and a good reproduction in the remaining
cases. Although the terms underestimation and overestimation are used, the general be-
havior remains fairly virtual because of the small difference between measurement and
calculation.

5. Conclusion

The study showed that the finite element method satisfactorily reflects the variation
of temperature fields at different depths of the soil. According to the results obtained, the
method gives a good simulation and explains to a precision that does not exceed 0.6 degree
maximum, that the model did not diverge even after 108 hours. Logically, the model will



REFERENCES 126

start to move away from the measure as maturities increase, but this is not the case for our
present model and this is due to the assumptions we made: Linear interpolation between
minimum and maximum ground temperature to obtain temperatures every 10 minutes;
the presence of clouds and dew each morning during the entire deadline, the presence of
rains, the presence of plants. We noticed after the simulations that the finite element
method estimates with good accuracy the temperature in the soil. The results are quite
conclusive as well with an initialization as without. On the other hand, we would like to
point out that the results of the simulation could be better if the measurements were made
under a greenhouse and thus protected from the hydrometeor which clearly influence the
variation of the surface temperature.

Acknowledgements

The first and second authors gratefully acknowledge Qassim University, represented
by the Deanship of Scientific Research, on the material support for this research under
number (1460-alrasscac-2016-1-12-S) during the academic year 1437 AH/ 2016 AD. Also
the authors would like to thank the anonymous referees and the handling editor for their
careful reading and for relevant remarks/suggestions which helped them to improve the
paper.

References

[1] M. Ainsworth, Dispersive and dissipative behavior of high order Discontinuous
Galerkin finite element methods. Journal of Computational Physics, 198(1):106–130,
2004.

[2] A. Arakawa, V.R. Lamb, Computational design of the basic dynamical processes of
the UCLA general circulation model Methods Comput. Phys., 17:173–265.

[3] S. Boulaaras and M. Haiour. “L∞-asymptotic behavior for a finite element approxi-
mation in parabolic quasi-variational inequalities related to impulse control problem”,
Applied Mathematics and Computation, Vol. 217, No. 3, 2011, pp. 6443–6450,

[4] 15- S. Boulaaras , M. Haiour: A General Case for the Maximum Norm Analysis of an
Overlapping Schwarz Methods of Evolutionary HJB Equation with Nonlinear Source
Terms with the Mixed Boundary Conditions, Applied Mathematics & Information
Sciences 9(3): (2015)1247-1257 ·

[5] J.H. Atkinson, J.J. Westerink, R.A. Luettich, Two dimensional dispersion analyses of
finite element approximations to the shallow water equations International Journal
for Numerical Methods in Fluids, 45:715–749, 2004.

[6] M.L. Batteen, Y.-J. Han, On the computational noise of finite-difference schemes used
in ocean models Tellus, 33:387–396, 1981.



[7] J.-M. Beckers and E. Deleersnijder, Stability of a fbtcs scheme applied to the propa-
gation of shallow-water inertia-gravity waves on various space grids. Journal of Com-
putational Physics, 108:95–104, 1993.

[8] P.-E. Bernard, J.-F. Remacle, V. Legat and E. Deleersnijder, Dispersion analysis
of Discontinuous Galerkin Schemes applied to Poincaré, Kelvin and Rossby waves
Journal of Scientific Computing, 34:26–47.

[9] B. Cushman-Roisin, Introduction to Geophysical Fluid Dynamics Prentice-Hall, 1994.

[10] P.G. Ciarlet, Introduction à l’analyse numérique matricielle et à l’optimisation Mas-
son, paris, 1990.

[11] J.K. Dukowicz Mesh effects for Rossby waves Journal of Computational Physics, 119:
188–194, 1995.

[12] R. Glowinski, Numerecal methods for Nolinear variational problems. Springer -Verlag,
1982.

[13] G. Guyot, Cours de Bioclimatologie, Chapitre II : Echanges de Chaleur et de Masse
par Conduction et Convection’, INRA Bioclimatologie, BP. 91, Montfavet cedex,
1992.

[14] M. Boulbrachene, M. Haiour, The finite element approximation of Hamilton Jacobi
Bellman equations, Comput. Math. Appl. (2001) 41, 993-1007.

[15] R. Haxaire, ‘Caractérisation et Modélisation des Ecoulements d’Air dans une Serre ’,
Thèse de Doctorat de l’Université de Nice, France, 1999.

[16] A. Majda. Introduction to PDE’s and waves for the atmosphere and ocean . American
Mathematical Society, 2003.

[17] R.Nebbali et S. Makhlouf, Dtermination de la distribution du champ de tempratures
dans le sol, par un modle semi-analytique.Conditions aux limites pour les besoins de
simulation d.une serre de culture, Universit Mouloud Mammeri, 15000 Tizi-Ouzou,
Algrie. 2000.

[18] D.A. Randall Geostrophic Adjustment and the Finite-Difference Shallow-Water Equa-
tions Monthly Weather Review, 122:1371–1377, 1994.

[19] V. Rostand, D.Y. Le Roux Raviart-Thomas and Brezzi-Douglas-Marini finite element
approximations of the shallow water equations International Journal for Numerical
Methods in Fluids, (published online), DOI:10.1002/fld.1668, 2007.

[20] L. White and E. Deleersnijder. Diagnoses of vertical transport in a three-dimensional
finiteelement model of the tidal circulation around an island. Estuarine Coastal and
Shelf Science, 74:655–669, 2007.

127



REFERENCES 128

6. ANNEX

6.1. Figures of without initialization results

Figure 1: Temperature evolution in the soil for different depths.

Figure 2: Comparison after 12h.

6.2. Figures of with initialization results



REFERENCES 129

Figure 3: Comparison after 24h.

Figure 4: Comparison after 36h.



REFERENCES 130

Figure 5: Comparison after 48h.

Figure 6: Comparison after 60h.



REFERENCES 131

Figure 7: Comparison after 72h.

Figure 8: Comparison after 84h.



REFERENCES 132

Figure 9: Comparison after 96h.

Figure 10: Comparison after 108h.



REFERENCES 133

Figure 11: Comparison after 12h.

Figure 12: Comparison after 24h.



REFERENCES 134

Figure 13: Comparison after 36h.

Figure 14: Comparison after 48h.



REFERENCES 135

Figure 15: Comparison after 60h.

Figure 16: Comparison after 72h.



REFERENCES 136

Figure 17: Comparison after 84h.

Figure 18: Comparison after 96h.



REFERENCES 137

Figure 19: Comparison after 108h.