Complexity Analysis of Unsaturated Flow in Heterogeneous Media
Using a Complex Network Approach
Hamed.O.Ghaffari, Mamdou Fall and Erman.Evgin
Dept. Civil Engineering, University of Ottawa, Ottawa, Ontario, Canada

Abstract: In this study, we investigate the complexity of two-phase flow (air/water) in a heterogeneous
soil sample by using a complex network theory .Based on the different similarity measurements (i.e.,
correlation, Euclidean metrics) over the emerged patterns from the evolution of saturation of non-wetting
phase of a multi-heterogeneous soil sample, the emerged complex networks are recognized.
Understanding of the properties of complex networks (such degree distribution, mean path length,
clustering coefficient) can be supposed as a way to analysis of variation of saturation profiles structures
(as the solution of finite element method on the coupled PDEs) where complexity is coming from the
changeable connection and links between assumed nodes. Also, the path of evolution of the supposed
system will be illustrated on the state space of networks either in correlation and Euclidean
measurements. The results of analysis showed in a closed system the designed complex networks
approach to a small world network where the mean path length and clustering coefficient are low and
high, respectively. As another result, the evolution of macro -states of system (such mean velocity of air
or pressure) can be scaled with characteristics of structure complexity of saturation. In other part, we tried
to find a phase transition criterion based on the variation of non-wetting phase velocity profiles over a
network which had been constructed over correlation distance.

Keywords: Unsaturated Flow, Complex Networks, heterogeneity.

1

1. Introduction

Complexity of flow in porous media–especially in multiphase state- is the result of multi
dimensional (geological) heterogeneity of particles where each element of system (or group of
elements) based on nonlinear mutual relation with gravity, capillarity and pressure gradient
shows its behaviour. In analysis of such rich complex behaviour apart of physical simplification
to solve the problem, the possible uncertainties in different format (spontaneously or nonspontaneously) will induce other facet of complexity. This type of complexity has allocated
considerable literature where uncertainty (i.e. stochastic or fuzzy) in direct form construct
stochastic differential equations or fuzzy differential equation while in indirect branch
uncertainty acts basically over collected information from system and gives inferences associated
with in(i.e. fuzzy) or on(i.e. rough set) boundaries of system [1],[17],[19].
Along this complexity, the interaction of particles which their dynamic directly are
governed by forces (agents) exhibits another complexity (absolutely coupling with uncertainty
and controlling equations): structure and topology based complexity. Undoubtedly, the later
complexity has an effective role in recognizing of state of each particle where the interactions of
the (local) states with each other and external forces give practical measurement of the system. In
other world, dynamical superposition of complexity in external forces, uncertainty in behavior
and interaction of particles may give a general evolution path of a system [1, 4]. The point that
must be considered is the role of local interactions of particles (in spatial or /and temporal form)
on the evolution of system so that such locality in different computational based methods (i.e.,
game theory, cellular automata, lattice Boltzmann methods, finite difference solutions, etc) has
been taken in to account. However, the idea behind the mentioned methods (as the basic and
initial idea) is considering the function of immediate neighborhoods[2-3].
Understanding of structures Complexity in variation and dynamic of a physical system
can be evaluated by constructing links (may be with weights) over nodes (as particles or agents).
This interpretation of topology complexity during the last years has opened a new perspective to
traditional graph theory, called complex networks [2]. Picturing, modeling and evaluation in a
simple and intuitive way are some of the discriminated features of complex networks [3, 4].
Based on the mentioned description of the topology structure over space and/or time, possible
2

explicit or implicit of agents /governing equations within or behind of complexity can lead to
networks on differential equations (NODs-NPDs) or linguistic rules such fuzzy rules (NFRs) [4]
.Complex networks have been developed in the several fields of science and engineering for
example social, information, technological, biological and earthquake networks are the main
distinguished networks [5-12]. Regard to the geosciences field, during past decade, some
considerable efforts have been done on the complex earthquake networks (generally recursive
events [7, 10] and climate dynamic networks. The dynamical growth of a thin film surface (and
in a nano-structural way) and evolution of a rock joint (complex aperture network-millimeter
scale) by this approach have been devised [8, 10].
The aim of this study is to analysis of complexity evolution of saturation (of air) profiles in a
closed heterogeneous soil sample (under mildly linear evolution of gas) either in spatial or
temporal way. In this way based on two metric spaces, firstly along each sampled time steps,
networks will be constructed. Characterization of the appeared networks and comparison of
network properties with each other and other distinguished mean properties of system (such
pressure and velocity) will be the next step of this study. In fact we are searching a suitable space
in where emerged properties of networks over one parameter can scaled with other dimension
(attribute) of system. In other word, the question “is there any hidden space so that structure
evolution of one property over that space is related to another attribute of same system”. Also,
we will make state spaces of complex networks where show the path of evolution of system
where are dissimilar for evaluated distances. The organization of this paper is as follow: in the
second part we will introduce the model parameters and accompanied governing equations; the
next section will cover a brief on complex networks properties and the results of the designated
networks on saturation and velocity profiles will be completed this section.
2.

Advection flow in porous media

Generally, the process of displacement of wetting/non-wetting phase is mainly affected by the
properties of the permeable medium, and fluids in single-phase both in homogeneous and in
heterogeneous media. In two-phase immiscible flows, interactions between the permeable
medium and the fluids also affect the fluid flow paths. Because flow dynamics depend on a
combination of conditions such as heterogeneity, moisture content, and chemistry, the resulting
3

transient flow and transport are usually complex. A large number of numerical methods have
been developed to model two-phase flow in heterogeneous media. The finite difference and finite
volume methods are the general frameworks for numerical simulation for the study of fluid flow
in very large problems. It has been proved that the finite difference method, however, is strongly
influenced by the mesh quality and orientation, which makes the method unattractive for
unstructured gridding [13-14].
Finite element methods (in different forms), have been used to model single-phase and
two-phase flow [13-17] in fractured media and heterogonous permeable media with different
capillarity pressures. In this part of study, based on the transition modeling of immiscible twophase flow, the role of permeability, porosity and the width of pore-size distribution in a small
sample of soil using FEMLAB [18] will be evaluated, while non-wetting phase (air) changes
with time, mildly. In all of this study, we assume that the porous medium is a non-deformable
(constant porosity) and that cross-product permeability terms associated with the viscous drag
tensor can be neglected. The general form of the two-fluid flow equations (without source–sink
terms) is described by the two-fluid, volume-averaged momentum and continuity equations [1517]:
qw = −

φ

kw

µw

∂ ( ρ w Sw )
+ ∇.( ρ wq w ) = 0
∂t

q nw = −

φ

[∇Pw + ρ wg ]

k nw

µnw

.[∇pnw + ρ nwg ]

∂ ( ρ nw S nw )
+ ∇.( ρ nw q nw ) = 0
∂t

(1)

(2)

(3)

(4)

In Eqs (1-4), the subscripts w and nw denote the wetting and non-wetting fluids,
respectively; Pi (i = w, nw) , Si , q i ,

g , µ i and k i denote pressure ,degree of saturation relative to

the porosity φ , the flux density vector , the gravitational acceleration vector, dynamic viscosity
, density and the effective permeability tensor, respectively. The effective permeability can be
defined as the relation between intrinsic permeability ( k ) and relative permeability ( k ri ):

ki = kri k .
4

With definition of volumetric fluid content as

θi = φSi we can write:

Sw + Snw = 1

(5)

θw + θnw = φ
Assuming one-dimensional vertical flow and that the wetting fluid is incompressible,
substitution of Eq. (1) into Eq. (2) gives:


∂θw ∂  kw ∂Pw
=  (
+ ρw g )
∂t ∂z  µw ∂z


(6)

When replacing fluid pressures and capillary pressure ( Pc = Pnw − Pw ) with pressure head

hi =

Pi

K =
ρH2O .g and defining the hydraulic conductivity of fluid i by i

ki ρ H2O . g

µi

, transient flow

of the wetting fluid by considering the capillary capacity ( C w = dθ w
) is described by :
dh
c

Cw (

∂hnw ∂hw
∂ 
∂h

−
) =  K w ( nw + 1) 
∂t
∂t
∂z 
∂z


(7)

Substitution of Eq.(3) in Eq.(4) yields :


∂( ρnwθnw ) ∂  ρ nwknw ∂Pnw
= 
(
+ ρ nw g )
∂t
∂z  µnw
∂z


(8)

For air, the density of non –wetting fluid has a dependency to the pressure head:
ρ nw = ρ 0 ,nw + (

ρ 0 ,nw
h0

(9)

) hnw

where ρ0,nw , ρ0,nw are the reference density and pressure head (at atmospheric pressure),
respectively. The ratio of

ρ 0,nw
h0

is defined as the compressibility λ (λ = 1.24 × 10−6 g

cm4

) . Then,

the flow of non-wetting fluid can be concluded:

{(φ − θ w )λ − ρ nwCw }

∂hnw
∂h
∂ 
∂h
ρ 
+ ρ nwC w w =  ρ nw K nw ( nw + nw 
ρ H 2O 
∂t
∂t
∂z 
∂z

5

(10)

Equations (7) and (10) for the pre-defined boundary conditions can be solved
simultaneously (coupled PDES). By functional description of the capillary pressure–saturation,
hc(Sw), and permeability functions, ki(Sw), the evolution of wetting and non-wetting phases
distribution can be estimated.. In this study, we use the van Genuchten equation (VG) [21]:

Sew = [1 + (αhc )n ]−m
where Sew denotes the effective saturation of the wetting fluid, Sew = (θ w − θ wr )
where

θwr and θwr

(11)

(θ ws − θ wr )

are the saturated and residual wetting fluid saturation, respectively;

α

,

and n

are fitting parameters, that are inversely proportional to the non-wetting fluid entry pressure
value and the width of pore-size distribution, respectively. We assume that m = 1 − 1 , and that
n
the effective saturation of the non-wetting fluid (Sen) is derived from Sen =1-Sew.
The capillary pressure–saturation function can be considered a static soil property, while
the permeability function is a hydrodynamic property describing the ability of the soil to conduct
a fluid. The basic assumption behind Capillary pressure–permeability prediction models are from
conceptual models of flow in capillary tubes combined with pore-size distribution knowledge
which are derived from the capillary pressure–saturation relationship. A typical representation of
this type of model follows Mualem formulation [22]:

 Se dSe 
∫

η  0 hc 
krw = Sew 1
 dSe 
∫

 0 hc 

2

 1 dS e
∫
hc
η  Se
k rn = (1 − S ew ) 1
 dS
∫ e
 0 hc

6

(12a)








2

(12b)

Combining the van Genuchten capillary pressure–saturation Eq. (11) with the Mualem
(VGM) model-with introducing new parameter the tortuosity parameter ( η ) gives permeability
functions as follows by [15,22]:

k rw

k rnw

1

k
η 
m m
= w = S ew
1
−
(1
−
S
ew ) 

k



2

1


k
m
= nw = (1 − S ew )η 1 − S ew

k



(13a)

2m

(13b)

Combination Eqs. (8) and (10), the boundary and initial conditions (Fig. 1), and the
constitutive relationships in Eq.(12-13) constitute the mathematical model of the assumed
system. In the upper bound of the geometry (fig. 1), a slight growth of the air pressure in 0 to
1hour (see Fig. 2) will be increased where the lower bound of the sample is in-transferable to
non-wetting phase (Fig. 1). The experiment injects air over the surface of a laboratory column
filled with water and sand. The incoming air (the nonwetting phase for this fluid pair) forces the
water (the wetting phase) toward the outlet at the base of the column. At the inlet, air pressure
increases by steps in time, and no water exits through the column top. In moving to the outlet,
the water passes through a disc that is impermeable to air flow. Neither the air nor the water can
pass through the vertical column walls. The water pressure at the outlet, which changes in time,
corresponds to the height of fluid rise in a receiving buret. The column has a total length of 8.34
cm, a 6-cm radius, and the disk is 0.74 cm thick [15]. Initial non-wetting phase pressure head is
assumed 0.2 (m water) and at the water outlet, the fluid level increases linearly in time from 0 m
to 0.1 m. The employed parameters for the sample come from Lincoln sand (88.6% sand, 9.4%
silt and 2.0% clay). The constant parameters for our model are as follow (see figure 3 for
uncorrelated normal distribution of random parameters):

θr (cm3cm−3 ) = 0.0210


θ s = [θ s ];



2
−9
k (cm ×10 ) = [kij ] 


−1
α (cm ) = 0.0189

n = [nij ]




7

in which θs , k and n will follow an uncorrelated Gaussian distribution drawn in Figure 3.
y

qw = 0, pnw = pnw (t )

7.6 cm
top
pnw = pnw (t0 ) = pnw
+ (8.34(cm) − y ) ρ nw g

pw = pw (t0 ) = − ρ w gy

Impermeable layer for air
0.74 cm

qnw = 0; pw = pw (t )
6cm
Figure 1. Schematic representation of boundary and initial conditions [15,18].

8

x

6.59
6.58

Gas Pressure (H)

6.57
6.56
6.55
6.54
6.53
6.52
6.51
6.5
0

0.2

0.4

0.6

0.8

1

Time(h)

Figure 2. Gas evolution transition on a supposed sample

m2

Figure 3. Intrinsic Permeability variation based on a normal distribution

9

250

250

t=0

250

t=0.5 h

250

t=0.8h

200

200

200

200

150

150

150

150

100

100

100

100

50

50

50

50

50

100

150

200

250

50

100

150

200

250

50

100

150

200

250

t=1 h

50

100

150

200

S_nw
1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

Figure 4. Effective saturation of non-wetting phase evolution

(h)
Figure 5. Absolute air-velocity fluctuation (inset: effective saturation of non-wetting phase evolution along time)

10

250

(h)
Figure 6. Air pressure changes vs. time for four randomly selected points

Figures 4-6 represent the result of

simulation with 1857 triangular elements for

saturation of non-wetting phase ,absolute value of air-velocity and pressure of non-wetting
phase evolution, respectively. Based on the later mentioned figure it seems in association with
initial boundary condition a phase transition before 0.1 h can be observed. The point of phase
transition is known as a critical point to include some basic terminology. Usually two types of
phase transitions are distinguished, first order and second order. First order designates phase
transition where the macroscopic states change in a discontinuous way upon passage through the
critical point, and second order indicates phase transitions where the states change in a
continuous way [23]. When a system passes through a second order transition; it may be left at
that point. At this point the system is like a pencil balanced on its end (only small variation may
be observed). We cannot tell which direction it will fall, but a small perturbation can send it
falling in a certain direction. In our case, the phase transition can be pronounced around 0.074 h
where overall signals of the variables states show an indented behaviour change.

11

Frequency of normalized velocities

Normalized total velocity
Figure 7. frequency distribution of the total velocity over time in log-log coordinate

Depicting of normalized (between 0 and 1) velocity frequency (Fig.7) shows the general overall
probability patterns of locally scaled- velocity. Roughly speaking, It can be followed that the
frequency of displacements rates of non-wetting phase is well approximated by a truncated
power-law [23, 24]: N ( v ) = ( v + v0 ) − β exp( − v ) ; where v0 = 0.05; β =1.5;κ = 0.85 . This equation
κ
suggests that air particles motion follows approximately a truncated Levy flight (Levy
processes), indicating that, despite the diversity of their travel history rate, the scaled rate of
displacement of non-wetting phase particles follow reproducible patterns. The main property of
the levy walk procedures is obeying from the infinite variance despite of the Gaussian process
which has a finite variance. However, extension of Levy flight to truncated Levy flight transfers
the infinite variance to finite variance [25-26].

12

3. Complex networks

As we mentioned in the previous section, the emerged patterns are following relatively
complex forms. Such complex forms come out in the signals and time series formats. Recently,
one of the methods for analysis of signals (information patterns) which has been considered is
employing of complex networks over the time-space (spatial/temporal). Different methodologies
have been proposed to capture the signal evolution based on the networks formation [31].
A network (graph) consists of nodes and edges connecting them. To set up a network, we
consider each profile (in y-direction- saturation or pressure of non wetting phase) as a node. To
make edge between two nodes, a relation should be defined. Several similarity or metric spaces
has been proposed for a construction of a proper network. The main point in selection of each
space is to explore the explicit or implicit hidden relation among different distributed elements of
a system. In this study we will use p-value criterion (a matrix for testing the hypothesis of no
correlation) related to correlation measurement and Euclidean distance over the non-wetting
saturation profiles. For each pair of signals (profiles) Vi and V j , containing L elements (pixels) the
correlation coefficient can be written as [30-31]:
L

∑ [V (k )− ≺ V  ]. V ( k )− ≺ V
i

Cij =

i

j

j

 

(14)

k =1
L

L

∑ [V ( k )− ≺ V  ] . ∑ [V (k )− ≺ V  ]
2

i

i

2

i

k =1

i

k =1

L

∑V (k )
i

where ≺ Vi  =

k =1

L

. Obviously,

Cij is restricted to the −1 ≤ Cij ≤ 1 , where Cij = 1,0 and −1

are related to perfect correlations, no correlations and perfect anti-correlations, respectively.
Each p-value is the probability of getting a correlation as large as the observed value by random
chance, when the true correlation is zero. The p-value is computed by transforming the
correlation to create a t- statistic having N-2 degrees of freedom, where N shows number of
profiles (samples) [27]. If p is small-here less than p ≤ ξ =0.05- then the correlation Cij is
significant. Euclidean distance is given as:

13

L

d ij =

∑ (V (k ) − V (k ))
i

2

(15)

j

k =1

Selection of threshold ( ξ ) is a challengeable discussion that can be seen from different
view. Choosing of such constant value may be associated with the current accuracy at data
accumulation where after a maximum threshold the system loses its dominant order. In fact,
there is not any unique way in selection of constant value, however, preserving of general
patterns of evolution must be considered while the hidden patterns (in our study come out from
the solution of two-PDEs based on a two-phase flow constitutive model) can be related to the
several characters of a networks . These characters can express different facets of the relations,
connectivity, assortivity (hubness), centrality, grouping and other properties of nodes and/or
edges. Generally, it seems obtaining stable patterns of evolution (not absolute) over a variation of

ξ can give a suitable and reasonable formed network. For the Euclidean distance we set
d ij ≥ ξ = 0.7 d ijmax . Consider with this definition we are investigating dissimilarity of profiles over
the metric space while for former metric the variations of profiles related to each other are
measured.
Let us introduce some properties of the networks: clustering coefficient (C) , the degree
distribution ( P (k ) ) and average path length ( L ). The clustering coefficient describes the
degree to which k neighbors of a particular node are connected to each other. Our mean about
neighbors is the connected nodes to the particular node. To better understanding of this
concept the question “are my friends also friends of each other?” can be used. In fact
clustering coefficient shows the collaboration between the connected nodes to one. Assume
the i

th

node to have k i neighboring nodes. There can exist at most k i (k i −1) / 2 edges

between the neighbors (local complete graph). Define c i as the ratio

ci =

actual

number

of

edges

between the neighbors
k i ( k i − 1) / 2

of

the

i th

node

(16)

Then, the clustering coefficient is given by the average of c i over all the nodes in the
network [7,10]:

14

C =

For k i ≤ 1 we define C

1
N

N

∑c .
i

(17)

i =1

≡ 0 . The closer C is to one the larger is the interconnectedness of the

network. The connectivity distribution (or degree distribution), P (k ) is the probability of finding
nodes with k edges in a network. In large networks, there will always be some fluctuations in the
degree distribution. The large fluctuations from the average value ( < k > ) refers to the highly
heterogeneous networks while homogeneous networks display low fluctuations [6]. The average
(characteristic) path length L is the mean length of the shortest paths connecting any two nodes
on the graph. The shortest path between a pair (i, j) of nodes in a network can be assumed as
their geodesic distance gij , with a mean geodesic distance L given as below [11]:
L=

2
∑ gij ,
N ( N − 1) i < j

(18)

where gij is the geodesic distance (shortest distance) between node i and j, and N is the number
of nodes. We will use a well known algorithm in finding the shortest paths presented by Dijkstra
[28]. Based on the mentioned characteristics of networks two lower and upper boards of
networks can be recognized: regular networks and random networks (or Erd˝os-Renyi networks
[5-6]). Regular networks have a high clustering coefficient (C ≈ 3/4) and a long average path
length. Random networks (construction based on random connection of nodes) have a low
clustering coefficient and the shortest possible average path length. However Watts and Strogatz
[30] introduced a new type of networks with high clustering coefficient and small (much smaller
than the regular ones) average path length (this is called SW property).
To extraction of saturation networks we set up X=273*Y=273 points on the constrained
area (only upper part-see figure 1). During the evolution of system, in 10 time-points using
aforementioned method, complex networks along Y-direction and on the non-wetting phase were
elicited. Consider that the numbers of evaluated points for construction of networks are smaller
and bigger than the number of (finite) elements and heterogeneous elements which give the
solution. This implies a limited correlation around each granule (node) with neighbourhood
nodes. This can be followed by the evolution of adjacency matrix visualization where for
15

correlation measurement and Euclidean distance have been depicted at figure 8. Appearance of a
thin layer of connections (disconnections for Euclidean distance) around diagonal of matrix
along successive time steps suggests a correlation around each node. It seems that the variation
of the emerged network over the p-value measurement is insignificant whereas networks based
on dissimilarity of Euclidean metric shows much higher variations. Evolution of number of
edges, mean clustering coefficient and average path length can be followed in Figures 9 -11,
respectively. As one can ensue the pattern and general trend of clustering coefficient and edges
are similar for each case and average path length shows an opposite variation with the later
properties. Then it can be interpreted, in mean view, increasing or decreasing of edges rate is
near to forming or decaying of triangles (clustering coefficient). Also, high clustering coefficient
and low average path length (especially for Euclidean measurement at last time steps) can be
distinguished (property of small-world networks) [6, 27]. Notice the fluctuation of properties of
saturation networks are not as well as the change of saturation (Fig.5), but is so similar to
absolute value of air velocity. Roughly speaking, one can infer from comparison of velocity and
L , < ci > , K

that

maximum

velocity

due

to

numerical

solution

coincides

with

=0.005−0.01
t =0.005−0.01
< ci >tcorrelation
, Kcorrelation
. It is noteworthy before t<0.02 the properties of two distances are

reveres each other but after this step the general trends are same so that we can write:

max{K , C , L}Euclidean = t{0.08, 0.08, 0.02}
min{K , C , L}Euclidean = t{0.02, 0.02,0.08}
max{K , C , L}p-value = t{0.01,0.01, 0.02}

.

min{K , C , L}p-value = t{0.02, 0.02, 0.01}
Consider that maximum (minimum) value of K and C are adapted to minimum
(maximum) L and vice versa. If we consider t=0.02 as transition point as has been depicted in
figure 5; this gives the point (or stage) with lowest stability which may indicate lowest
correlation (t=0.02). The temporal evolution of clustering coefficient and edges over Euclidean
distance so approve this stage is a much more chaotic step where after this point the patterns of
spatial and temporal characters of structure complexity approximately reach a quasi stationary
step. To understand the path evolution of saturation profiles, we construct state diagram of
k i − ci for each case (Fig.14). For correlation distance (Fig.14 a) nodes generally (without

16

considering time effects) lies from higher clustering coefficient to a congested area where in this
area (can be called equilibrium place) we have 100 < kicore < 120;0.4 < cicore < 0.45 .
1

1
250

0.9

t=0

250

0.9

t=0.08 h

200

200

0.7

0.7
0.6

150

0.5

0.4

0.3

0.3
0.2

50

100

150

200

250

Correlation Distance

250

0.7

0.7

0.6

0.6

150

150

0.5

0.5
0.4

100

0.4

100

0.3

0.3
0.2

50

0.2

50

0.2

50

0.1

0.1
50

0.8
200

200

0.4

100

0.9

t=0.08 h

0.8

0.6
150

0.5

100

250
0.9

t=0

0.8

0.8

1

1
250

0.1

0.1
50

0

100

150

200

250

0

50

100

150

200

250

t=1 h

250

0

50

t=1h

100
1

150

200

250

0

0.9

Euclidean Distance

0.8

200

200

150

150

0.7
0.6
0.5

a)

b)

100

0.4

100

0.3
0.2

50

50

0.1
50

100

150

200

250

50

100

150

200

250

0

Figure 8. Evolution of adjacency matrix visualization of non-wetting phase saturation along time: based on a) p-value of
Correlation measurements and b) Euclidean Distance

(h)

Figure 9. Variations of number of edges in complex saturation networks vs. time: based on a) Euclidean Distance and b) pvalue of Correlation measurements

17

0.7

0.45

0.65

0.6

b) Correlation Distance

0.445

a) Euclidean Distance

0.55

<c>

< C>

0.44
0.5

0.45

0.435

0.4

0.43

0.35

0

0.1

0.2

0.3

0.4

0.5
t(h)

0.6

0.7

0.8

0.9

1

0

0.2

0.4

0.6

0.8

1

t

Figure 10. Variations of mean clustering coefficient in complex saturation networks vs. time: based on a) Euclidean Distance
and b) p-value of Correlation measurements
1.625

1.95
1.9

b) Correlation Distance

1.62

a) Euclidean Distance

1.615

1.8

mean geodesic distance

m ean geodesic distance

1.85

1.75
1.7
1.65
1.6

1.61

1.605

1.6

1.55

1.595

1.5

0

0.1

0.2

0.3

0.4

0.5
t(h)

0.6

0.7

0.8

0.9

1.59

1

0

0.1

0.2

0.3

0.4

0.5
t(h)

0.6

0.7

0.8

0.9

1

Figure 11. Evolution of mean geodesic distance (mean length of the shortest paths) of complex saturation networks along
time: a) Euclidean Distance and b) p-value of Correlation measurements

18

c

i

1

t=0

0.5
0
1

t=0.005

0.5

t=0.01

0.5

t=0.02

0.5

t=0.04

0.5

0
1
0
1
0
1
0
1

t=0.08

0.5

t=0.1

0.5

0
1
0
1

t=0.5

0.5

t=0.8

0.5

t=1

0.5

0
1
0
1
0
0

0.5

1

0

t=0

0.5

1

0

t=0.005

Figure 12. Temporal change of

0.5

1

0

t=0.01

0.5

1

0

t=0.02

0.5

1

0

t=0.04

0.5

1

0

t=0.08

0.5

1

0

t=0.1

0.5

1

0

t=0.5

0.5

1

t=0.8

0

0.5

1

t=1

c i related to different reference point s in association with Euclidean Distance

ki
400

t=0

200
0
400

t=0.005

200

t=0.01

200

0
400
0
200

t=0.02

100
0
400

t=0.04

200
0
400

t=0.08

200
0
400

t=0.1

200
0
400

t=0.5
t=0.8

200
0
400
200
0
400
200
0
0

200

t=0

400 0

200

400 0

t=0.005

Figure 13. Temporal change of

200

400 0

t=0.01

100

200 0

t=0.02

200

400 0

t=0.04

200

400 0

t=0.08

200

400 0

t=0.1

200

t=0.5

400 0

200

400 0

t=0.8

200

t=1

k i related to different reference points (time) in association with Euclidean Distance

19

400

data1
data2
data3
data4
data5
data6
data7
data8
data9
data10

-0.3

ci

ci
ci

c(i)

10

-1

10

-0.4

10

a) Correlation Distance

data1
data2
data3
data4
data5
data6
data7
data8
data9
data10

2

1

10
k(i)

k

b) Euclidean Distance

10

i

Figure 14. Diagram of structure complexity space over two variables of networks

2

k

10
kii

k i and c i : a) p-value of Correlation, b)

Euclidean Distance-data 1 to 10 are the sampling values within the time steps.

Then, roughly speaking, the path evolution of saturation profiles in correlation
measurement shows a narrow variation over the complexity structure space so that at stationary
state the core of space can be imagined as a point of absorption. It must be considered that the
discrimination of time steps over this evolutionary grape shape path is difficult. Despite of the
former metric, Euclidean distance can give irregularities over time and space (Fig14b). Figure
14.b shows in t=0.02 there is distinguished chaotic separately different from overall way of
change of patterns which is matched with pervious results. Also, the way of evolution can be
estimated by: ci ≈ 2.61ki−0.32 . One may modify this relation with a truncated power law. In another
view and with refer to our results drawn in Figure 14; we can conclude that the correlated and
dissimilar patterns of saturation profiles (at a closed system) in a proper space of networks may
be assumed as synchronizing oscillators. In other word, based on anharmonic oscillators’ theory
and the governing non linear differential equations we can write (for Figure 14):

d 2 ci
dc
= −ω 2 sin ci (t ) − α i + f cos ωt ,
dt
dt
in which ω , α and f are constant values.
Notice that with increasing of edges the communicability of system on this space is
decreasing which shows the dissimilarity patterns tend to growth over all on nodes rather that the
20

neighbourhoods of nodes. So, comparison with the distribution of edges indicates that nodes with
high correlation and dissimilarities are decaying as well as triangles densities around each node
(Figs.15&16). This implies that low dissimilar patterns (small k i ) have higher communicability.
Let us investigate the temporal structure complexity. Two immediate solution can be presented
in first situation matrix of connections along evolution of time is extracted. In this case spatial
connection of nodes is replaced with temporal evolution of same node with same node at next (or
pervious) time steps or other nodes. Consider a reference time step ( t ref . ) then we can make the
following matrix among t ref . , t ref . + ∆ t :
ref .

M

t ref .,t ref . +∆t

t +∆t


. . .
ref . 
= t . mi .
; mi ∈ {0,1}; i = 1,.., N
. . . N × N

(19)

where mi is the connection state of element i (position i) due to one relation along a time
interval. So, a multiple array over different time references can cover completely both spatial (
∆t = 0 ) and temporal evolution of structure complexity such:

 M (t0 ,t0 +∆t )

.
B = 
.

 M (t0 ,t0 + n×∆t )



.
 ; n = 1 + t f − t0

∆t
.

( t ,t + n×∆t )
 n×n
. . M f f
. .
. .
. .

M

( t f ,t f +∆t )

(20)

in which n is the time granulation. Consider for ∆t = 0 ,we can see only evolution of
spatial complexity on

the assumed time and in collective way as is usual in networks

terminology: [ B ]1× n . The mentioned method can give general complexity patterns of time - space
and interwoven structures, explicitly. Analysis of the attributes and properties of the B is out of
the scope of this study. Second case is an implicit (more qualitatively) view to time-based
complexity. In this state, the spatial characteristics ( k i , c i ) of the system at several time steps
will be compared (see figures 12 and 13).
In another procedure we construct air-velocity networks. The procedure is as well as
correlation distance, only the number of nodes has decreased to 93. Three distinguished
21

behaviours due to the nodes and edges variation are recognized, approximately (Fig.17 a,
Fig.18a): 1) decreasing of nodes where the numbers of edges are approximately fixed (t= [0.01,
0.02]) compatibility in variations and 3) vertices are constant while sensitive variation in links
can be perused (t=[0.3,1]). The behaviour of emerged network especially after (t>0.2) can be
estimated by the properties of small-world networks. The shape of the degree distribution is
similar to that of a random graph. It has a pronounced peak at ≺ k t >0.2 ≈ 8 while the topology of
the network is relatively homogeneous and all nodes having approximately the same number of
edges (Fig.18 b).
Another point can be inferred from the plotting of the reverse of mean clustering coefficient so
that in semi-logarithmic coordinate (Fig17.c) it follows a phase-transition (first order) function
that can be followed by a sigmoid function [22]:

≺ 1/ c ∝ (1 + e
where β and

Lnt
− β (1+1+δ Lnt )

)

(21)

δ are the regulator parameters which determine the declining rate. Then, the reverse

of mean clustering coefficient for velocity networks gives an acceptable phase change criteria.
-3

x 10

0.035

z data
fit 3

a) Correlation
Distance

0.03

z data
fit 1

b) Euclidean
Distance

7

PDF
6

0.025
5

i

Density

PDF 0.02
0.015

4

3

0.01

2

0.005

1

0

0

60

70

80

90

100

110
Data

120

130

140

150

0

50

100

150

ki

ki

200

250

k

Figure15. Edge distribution of complex saturation networks: based on a) p-value of Correlation measurements and b)
Euclidean Distance

22

N (ci )

Figure 16. Evolution of mean clustering coefficient distribution along time: based on a) p-value of Correlation measurements
and b) Euclidean Distance (inset: probability density function at 0.6<t <1 h)
1600

0.7

a)

1400

b)

0.6

1200

0.5
1000

< C>

edges

0.4
800

0.3
600

0.2
400

0.1

200

0

0

0.1

0.2

0.3

0.4

0.5
time

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0.5
time

0.6

0.7

0.8

0.9

1

45

c)

40

35

30

1./C

25

20

15

10

5

0
-3
10

-2

-1

10

10

0

10

time

Figure 17. a) Evolution of total number of edges for air velocity networks , b) Mean cluster coefficient variation and c) inverse
of Mean cluster coefficient variation

23

Sampled time steps [0,1]

N(k)

Number of Nodes

b)

a)

Sampled time steps [0,1]

k

Figure 18.a) Number of nodes-time ;b) frequency of edges evolution along time for velocity (of air) networks

4. Conclusion

Analysis of two phase flow especially in heterogeneous media, has allocated a noticeable
research literature in porous media area. To take in to account the complexity due to
heterogeneity in different parameters which are characterizing the general reaction (s) of the field
(s), different methodologies have been proposed. The aim of this study was to investigate the
appeared complexity of two-phase flow (air/water) in a heterogeneous soil. The supposed filed
was under the small range of air pressure variation. By considering the capillary pressuresaturation, permeability functions, as constitutive relations and governing equations (obtained
separately for evolution of the model’s wetting and non-wetting phase parameters) the unknown
parameters can be estimated. In this way, using finite element method the role of multiheterogeneity the successive change of different variables (such relative permeability, saturation,
capillary pressure etc) was analysed. Based on the different similarity measurements (i.e., correlation,
Euclidean metrics) over the properties of the emerged patterns of the complex networks were recognized.
In this way, the path of evolution of the supposed system was illustrated based on the state space of
networks either in correlation and Euclidean measurements. The results of analysis showed in a closed
system the designed complex networks approach to small world network property where the mean path
length and clustering coefficient are low and high. As another result, we tried to scale the evolution of
macro -states of system (such mean velocity of air or pressure) with characteristics of structure
complexity of saturation.

24

5. References
1.

Doering R.C.1991.Modeling Complex Systems: Stochastic processes, Stochastic Differential Equations, and Fokker-Planck
Equations; 1990 Lectures in complex systems, SFI Studies in the science of complexity.

2.

Arenas A., Diaz-Guilera A., Kurths J., Moreno Y., Zhou C.2008. Synchronization in complex networks .Physics Reports
469,93-153

3.

Dorogovetev S.N., Goltsev A.V.2008.Crtitical phenomena in complex networks .Review Modern Physics.Vol.80.OctDec.2008

4.

Boccaletti S.,Latora V., Moreno Y.,Chavez M.,Hwang D-U.2006.Complex Networks : Structure and Dynamics. Physics
Reports 424,175-308

5.

Newman M. E. J. 2003.The structure and function of complex networks, SIAM Review; 45(2): 167- 256.

6.

Albert R., Barabasi A.-L.2002. Statistical mechanics of complex networks. Review of Modern Physics; 74, 47–97.

7.

Abe S, Suzuki N.2006. Complex network description of seismicity. Nonlin process Geophys;13:145-150

8.

Ghaffari OH., Sharifzadeh M.,Fall M.2009. Analysis of Aperture Evolution in a Rock Joint Using a Complex Network
Approach. International Journal of Rock Mechanics and Mining Sciences, In Press, Corrected Proof, Available online 5
April 2009.

9.

Valentini L Perugini D, Poli G.2007. The small-world topology of rock fracture networks. Physica A; 377:323–328.

10. Baiesi M, Paczuski M. 2004.Scale free networks of earthquakes and aftershocks. Physical review E 2004; 69 (2): 066106.1066106.8
11. Tsonis, A. A., and P. J. Roebber, 2003: The architecture of the climate network. Physica, 333A, 497–504.
12. Yuksel M,Karabacak T, Guclu H. 2007.Networking behavior in thin film and nanostructure growth dynamics . 2nd IEEE
International Conference on Nano-Networks . Available at <http://arxiv.org/abs/0711.4745v1>
13. Hoteit H, Firoozabadi A. 2008. Numerical modeling of two-phase flow in heterogeneous permeable media with different
capillarity pressures; Advances in Water Resources 31, 56–73.
14. Lunati I, Jenny P. 2006. Multiscale finite-volume method for compressible multiphase flow in porous media. J Comput
Phys;216(2): 616–36
15. Chen J., Hopmans J.W. & Grismer M.E. 1999. Parameter estimation of two-fluid capillary pressure–saturation and
permeability functions; Advances in Water Resources Vol. 22, No. 5, 479–493.
16.
Kueper B. H., Frind E. O. 1991.Two-Phase Flow in Heterogeneous Porous Media Model Development; Water Resources
Research, vol. 27, 6, p 1049-1057.
17. Zhang D .2002. Stochastic methods for flow in porous media: coping with uncertainties. Academic, San Diego, ISBN 012–
7796215, p 350.
18. COMSOL Multiphysics. 2008b. COMSOL Multiphysics Modeling Guide: Version 3.5 Stockholm, Sweden: COMSOL AB.

19. Ghanem G., Dham S. Stochastic finite element analysis for multiphase flow in heterogeneous media. Transport Porous
Media 32:239–262, 1998.
20. Parker, J.C., R.J. Lenhard, and T. Kuppusany. 1987. A parametric model for constitutive properties governing multiphase
flow in porous media. Water Resour. Research 23:618:624

25

21. Van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci.
Soc. Amer. J. 44:892-898.

22. Mualem, Y. 1976. A new model for predicting the hydraulic permeability of unsaturated porous media. Water Resources
Research 12:513-522.
23. Millonas M.M. Swarms, phase transitions, and collective intelligence, 1993: http://arxiv.org/abs/adap-org/9306002.
24. González M. C., Hidalgo C. A., Barabási A.-L. 2008. Understanding individual human mobility patterns, Nature 453, 779782.
25. Klafter J., Zumofen G., Shlesinger M.F. 1995. Levy Flights and Related Topics in Physics. Lecture Notes in Physics, vol.
450, Springer, Berlin, p. 196.
26. Mantegna, R. N. & Stanley, H. E. 1994. Stochastic process with ultraslow convergence to a Gaussian: the truncated Levy
flight. Phys. Rev. Lett. 73, 2946–2949.
27. Watts DJ, Strogatz SH. Collective dynamics of small-world networks. Nature 1998; 393:440-442.
28. E.W. Dijkstra, Numerische Math. 1, 269 (1959)
29. Massey, F. J.1951. The Kolmogorov-Smirnov Test for Goodness of Fit, Journal of the American Statistical Association,
Vol. 46, No. 253, pp. 68–78.
30. Zhongke Gao and Ningde Jin.2009. Flow-pattern identification and nonlinear dynamics of gas-liquid two-phase flow in
complex networks; Physical Review E 79, 066303.
31. Zhang, J., Small, M.2006. Complex network from pseudoperiodic time series: Topology versus dynamics. Physical Review
Letters 96(238701).

26

