# Studying grout propagation in granular soils by image processing techniques

A. Aït Alaiwa, Nadia Saiyouri, Pierre-Yves Hicher

► To cite this version:

A. Aït Alaiwa, Nadia Saiyouri, Pierre-Yves Hicher. Studying grout propagation in granular soils by image processing techniques. Strain, 2011, 47 (2), pp.194-202. 10.1111/j.1475-1305.2008.00503.x. hal-01007348

HAL Id: hal-01007348

https://hal.science/hal-01007348

Submitted on 15 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Studying Grout Propagation in Granular Soils by Image Processing Techniques

A. Ait Alaiwa, N. Saiyouri and P.-Y. Hicher

Institut de Recherche en Génie Civil et Mécanique (GeM), UMR CNRS 6183, Ecole Centrale de Nantes, 1 rue de la Noë, BP 92101, F-44321 Nantes Cedex 3, France

This article presents two image-processing techniques which allow studying grout propagation in granular soils. The first one deals with the evaluation of the cement grout concentration using UV spectrophotometry. The second is based on image segmentation processing: active contour model. We propose to develop the active contour model image analysis technique. This process, currently applied in the field of medicine, can be transposed to grout flow detection during soil injection. They are powerful image segmentation techniques that combine geometry, physics and approximation theory. Two distinct formulations exist to employ these techniques, parametric or geometric curves. These models have proven to be effective in segmenting and tracking non-rigid structures. They exploit features derived from the image data together with a prior knowledge about the location, size and shape of these structures. We focus on a parametric approach for deformable models. They support highly intuitive interaction mechanisms. Those mechanisms allow, when necessary, researchers to bring their expertise to bear on the model-based image interpretation task. In this paper, we show the development of the two different image analysis methods employed in our experiments for grout detection and grout concentration measurement.

**KEY WORDS:** grout, image processing, porous media, segmentation

## Introduction

Deformable models or ‘snakes’ are object-delineating curves that move within two-dimensional (2D) digital images under the influence of both internal and external forces and user-defined constraints. Since their introduction by Kass et al. [1] or Malik and Perona [2], these algorithms have been at the heart of one of the most active and successful research areas in edge detection, image segmentation, shape modelling and visual tracking. There are two general types of deformable models in the literature today: parametric deformable models and geometric deformable models. In particular, we focus on parametric deformable models which are represented explicitly as parameterised contours (i.e. curves) in a Lagrangian framework.

Segmenting structures from images and reconstructing a compact geometric representation of these structures are difficult because of the small size of the data sets and the complexity and variability of the properties of medium studied, i.e. soil. Furthermore, the typical shortcomings of sampled data such as sampling artefacts, spatial aliasing and noise, may cause the boundaries of structures to be indistinct

and disconnected. The challenge is to extract boundary elements belonging to the same structure and integrate these elements into a coherent and consistent model of the structure.

This article presents a particular deformable model, a spectral method based on computer-assisted image analysis. The other deformable model techniques stem from their ability to segment, match and track images of various structures. Deformable models are capable of accommodating the often significant variability of soil structures over time. Furthermore, deformable models support highly intuitive interaction mechanisms. We will review the basic formulation of deformable models and their application to fundamental image analysis problems, including segmentation and matching.

As this approach does not have a high convergence rate, spectral methods were introduced. They allow application of real-time detection. The cement grout transport evolution will be shown as an application of this technique. Indeed, fluid displacement studies are important tools in soil physics for understanding the transport of adsorbed and non-adsorbed solutes through soil. Such experiments provide valuable insight into the porous media, the behaviour of
chemicals, and associated processes such as diffusion, dispersion and sorption during transport [3]. The mathematical conceptualisation of these processes at an individual pore scale is available (e.g. Navier-Stokes equation). However, it is not possible to completely characterise the variability or heterogeneity of the pore structure (micro-scale). Therefore, most models provide the description of solute transport at a macroscopic scale. Later, we focus on this last approach.

## Experimental Setup

The convective dispersion transport equation remains the foundation on which most analyses of solute transport in porous media have been based [4]. In this model, fluid velocities and concentration must be correctly evaluated. Therefore, we developed an application which allows tracking fluid motion and determining fluid concentration with respect to time. Fluid motion provides fluid velocity, and concentration must be evaluated to perform solute transport models. Figure 1 shows the scheme of the 1D laboratory for cement grout injection experiment.

For this study, dry Loire river sand was packed into polymethyl methacrylate plastic cylinders having a length of 0.9 m and a diameter of 0.10 m at different defined densities. Care was taken to follow exactly the same procedure for packing all the soil columns. Each soil column was first saturated by water from the bottom using a point source. Cement grout transport evolution was controlled by monitoring the weight of the column, the pressure at different heights and effluent volumes with respect to time (Figure 1). We used a charge-coupled device (CCD) camera to analyse grout concentration during the test. Later, cement grout was injected through the column. Cement grout position was tracked by using the active contour model. Cement grout concentra-

<figure><img src="image_1.png" /><figcaption>Figure 1: Experimental setup</figcaption></figure>

tion along the column was measured and correlated with a convection-dispersion model to identify the dispersive parameters. We present the classical solute transport model mostly used to study fluid propagation in the porous media.

## Transport and Detection Theories

### Convection-dispersion model

Complex behaviours were frequently observed in transport process because of heterogeneities in the porous media (grain size). Therefore, developing and using alternative ways to quantify such complex transport process in porous media have become one of the major research topics for transport modelling in the last several decades.

The hydrodynamic dispersion theory [3] is one of the most commonly used for describing the solute transport in porous media. The convection-dispersion equation (CDE) deduced from the theory is based on two major assumptions. The first one is that the mass centre of the solute plume travels with the velocity of the average fluid. The second assumption is that the dispersion behaves macroscopically as a Fickian diffusive process, with the dispersivity being assumed constant in space and time. For 1D flow, the CDE for the non-reactive solute is given by:

$$ \frac{\partial}{\partial t} c(x, t) = \frac{\partial}{\partial x} \left( D \frac{\partial}{\partial x} c(x, t) \right) - \frac{\partial}{\partial x} (v_f c(x, t)) \quad (1) $$

where c is the solute concentration; D the hydrodynamic dispersion coefficient; v_f the average pore velocity of solute transport equal to q/n, where q is the volumetric flux per unit area of porous medium, n the porosity; x the distance with respect to the injection point; and t the time. The CDE has been found with reasonable accuracy in describing the solute migration in aquifers with relatively small degrees of heterogeneity [5, 6] and in small-scale experiments [3].

### Absorption and reflection

Absorbance and reflectance spectroscopy commonly referred to spectrophotometry are analytical techniques based on measuring the amount of light absorbed and reflected by a sample. Particularly, in the visible and UV portions of the electromagnetic spectrum, this technique is one of the most widely used in chemistry.

Molecular absorption and reflectance spectroscopy in the visible spectrum (VS) is concerned with the
measured absorption of radiation in its passage through a liquid or a solid. The wavelength region generally used is from 190 to about 1000 nm, and the absorbing medium is at room temperature (20 °C).

In our case, we try to evaluate the concentration of a homogenous suspension of cement particles, which percolate in a homogenous sand column, by using the reflectance technique. As only the transmission coefficient of the suspension varies in the column, it is possible to relate cement particle concentration and the reflected radiant power.

## Beer-Lambert law

The Beer-Lambert law states that the concentration of a substance in solution is directly proportional to the 'absorbance A' of the solution as:

$$A = -\log\left(\frac{I_t}{I_0}\right) = \alpha \cdot L \cdot C \quad (2)$$

where $I_t$ is the radiant power transmitted by the absorbing medium, $I_0$ the radiant power incident on the medium, $\alpha$ the extinction coefficient – a constant dependent only on the nature of the solute, L the length of the radiation path through the sample, A the (decadic) absorbance and C the solute concentration.

The Beer-Lambert law holds only if the absorbing species behave independently of each other and if the absorption occurs in a uniform medium. Indeed, during the grout injection process, there is a non-optical interaction between the grout and soil. Photons, after crossing to interstitial fluid, are attenuated and reflected from grain surfaces. The measure of the power spectrum of these reflected photons depends directly on the concentration of cement particles if we consider the Snell-Descartes optical relations. We obtain the equivalent equation:

$$A' = -\log\left(\frac{I_r}{I_0}\right) = \alpha' \cdot L \cdot C \quad (3)$$

where $I_r$ is the radiant power reflected by the absorbing medium, $\alpha'$ is the reflection coefficient always dependent on the solute and A' is the 'reflection'.

As CCD cameras using International Telecommunications Union (ITU) standard (analogic capture) detect radiant power spectrum, we used the differential expression of Equation (3) and calibration of reflected power spectrum for several cement grout concentrations in order to get the concentration curve.

## Active contour models

The basic premise of the deformable contours is to find a parameterised curve that minimises the weighted sum of internal energy and potential energy contained in an image. The internal energy specifies the tension or the smoothness of the contour. The potential energy is defined over the image domain and typically possesses local minima at the image intensity edges occurring at object boundaries. Minimising the total energy yields internal forces and potential forces. Internal forces hold the curve together (elasticity forces) and keep it from bending too much (bending forces). External forces attract the curve towards the desired object boundaries. To find the object boundary, parametric curves are initialised within the image domain and are forced to move towards the potential energy minima under the influence of both these forces.

We first recall some definitions and formulate the mathematical problem. Mathematically, a deformable contour is a curve $v(s) = (x(s); y(s))$, $s \in [0, 1]$ in which $s$ is a parameter that increases as the curve is traversed. $x$ and $y$ are particular functions of $s$, which describe Cartesian coordinates of each point of the curve, which moves through the spatial domain of an image to minimise the following energy functional:

$$E(v) = S(v) + P(v) \quad (4)$$

The first term is the internal energy functional [7, 8] and is defined as

$$S(v) = \int_{0}^{1} w_{1}(s) \left| \frac{\partial v}{\partial s} \right|^{2} + w_{2}(s) \left| \frac{\partial^{2} v}{\partial s^{2}} \right|^{2} ds \quad (5)$$

The first-order derivative discourages stretching and makes the model behave like an elastic string. The second-order derivative discourages bending and makes the model behave like a rigid rod. The weighting parameters $w_1$ and $w_2$ can be used to control the strength of the model's tension and rigidity, respectively. In practice, $w_1$ and $w_2$ are often chosen as the constants.

The second term is the potential energy functional and is computed by integrating a potential energy function $P(x,y)$ along the contour $v(s)$:

$$P(v) = \int_{0}^{1} P(v(s)) ds \quad (6)$$

The potential energy function $P(x,y)$ is derived from the image data and takes smaller values at object boundaries as well as other features of interest.
Given a grey-level image $I(x,y)$ viewed as a function of continuous position variables $(x,y)$, a typical potential energy function designed to lead a deformable contour towards step edges is

$$P(x, y) = w_{\sigma}|\mathbf{G}_{\sigma}(x, y) * \nabla \mathbf{I}(x, y)|^2 \quad (7)$$

where $w_σ$ is a positive weighting parameter, $G_σ(x,y)$ is a 2D Gaussian function with standard deviation $σ$, $∇$ is the gradient operator, and '*' is the 2D image convolution operator. If the desired image features are lines, then the appropriate potential energy function can be defined as follows:

$$P(x, y) = w_l[\mathbf{G}_{\sigma}(x, y) * \mathbf{I}(x, y)], \quad (8)$$

where $w_l$ is a weighting parameter. Positive $w_l$ is used to find black lines on a white background, while negative $w_l$ is used to find white lines on a black background.

For both edge and line potential energies, increasing $σ$ can broaden its attraction range. However, larger $σ$ can also cause a shift in the boundary location, resulting in a less accurate result. Regardless of the selection of the exact potential energy function, the procedure for minimising the energy functional is the same. The problem of finding a curve $v(s)$ that minimises the energy functional $E$ is known as a variational problem [9]. It has been shown that the curve that minimises $E$ must satisfy the following Euler-Lagrange equation [1, 10]:

$$-\frac{\partial}{\partial s}\left(w_1(s) \frac{\partial \mathbf{v}}{\partial s}\right)+\frac{\partial^2}{\partial s^2}\left(w_2(s) \frac{\partial^2 \mathbf{v}}{\partial s^2}\right)+\nabla P(\mathbf{v}(s, t))=0 \quad (9)$$

To gain some insight into the physical behaviour of deformable contours, we can view Equation (9) as a force balance equation

$$\mathbf{F}_{\text{int}}(\mathbf{v}) + \mathbf{F}_{\text{pot}}(\mathbf{v}) = 0 \quad (10)$$

where the internal force is given by

$$\mathbf{F}_{\text{int}}(\mathbf{v}) = \frac{\partial}{\partial s} \left( w_1 \frac{\partial \mathbf{v}}{\partial s} \right) - \frac{\partial^2}{\partial s^2} \left( w_2 \frac{\partial^2 \mathbf{v}}{\partial s^2} \right) \quad (11)$$

and the potential force is given by

$$\mathbf{F}_{\text{pot}} = -\nabla P(\mathbf{v}) \quad (12)$$

The internal force $F_{\text{int}}$ discourages stretching and bending while the potential force $F_{\text{pot}}$ pulls the contour towards the desired object boundaries. In this chapter, we define the forces, derived from the

potential energy function $P(x,y)$ given in either Equation (7) or Equation (8), as Gaussian potential forces.

To find a solution to Equation (9), the deformable contour is made dynamic by treating $v(s)$ as a function of time $t$ as well as $s$, i.e. $v(s,t)$. The partial derivative of $v$ with respect to $t$ is then set equal to the left-hand side of Equation (9) as follows:

$$\gamma \frac{\partial \mathbf{v}}{\partial t} = \frac{\partial}{\partial s} \left( w_1 \frac{\partial \mathbf{v}}{\partial s} \right) - \frac{\partial^2}{\partial s^2} \left( w_2 \frac{\partial^2 \mathbf{v}}{\partial s^2} \right) - \nabla P(\mathbf{v}) \quad (13)$$

The coefficient $γ$ is introduced to make the units on the left side consistent with the right side. When the solution $v(s,t)$ stabilises, the left side vanishes and we achieve a solution of Equation (9).

We note that this approach of making the time derivative term vanish is equivalent to applying a gradient descent algorithm to find the local minimum of Equation (4) [11]. Thus, the minimisation is solved by placing an initial contour on the image domain and allowing it to deform according to Equation (13).

## Fourier spectral algorithm for the active contour model equation

In this section, the principles behind the Fourier spectral method will be discussed [12, 13] and an algorithm for implementation of the snake evolution equation using Fourier spectral method will be presented.

The basic idea of spectral methods for solving partial differential equations is to assume that $u(x)$, the equation solution can be approximated by a sum of $N + 1$ basic functions $φ_n(x)$

$$u(x) \approx u_N(x) = \sum_{n=0}^{N} a_n \phi_n(n) \quad (14)$$

When this series is substituted into the differential equation

$$Lu = f(x) \quad (15)$$

where $L$ is the differential operator, the result is the so-called residual function defined by

$$R(x_i, a_0, a_1, \dots, a_N) = Lu_N - f \quad (16)$$

If the basic functions individually satisfy the homogeneous boundary conditions on $u(x)$, then their sum will too. Therefore, the only error in $u_N(x)$ is that it does not exactly satisfy the differential equation. Since the residual function is identically 0 for the
exact solution, the challenge is to invent systematic methods for choosing the series coefficients $a_N$ to obtain residual function as small as possible. The different spectral methods differ in their way of minimising the residual.

When the boundary conditions require the solution to be spatially periodic, the solution generation procedure is normally much simpler than for non-periodic boundaries. In the periodic case, the sines and cosines of a Fourier series automatically and individually satisfy the boundary conditions. Consequently, our only remaining task is to choose the coefficients of the Fourier series to minimise the residual function.

Consider a set of points

$$x_j = \frac{2j\pi}{N}, \quad j \in \mathbb{N} \cap [0, N-1] \qquad (17)$$

which are called collocation points. The discrete Fourier coefficients are

$$\tilde{u}_k = \frac{1}{N} \sum_{j=0}^{N-1} u(x_j) e^{-ikx}, \quad -\frac{N}{2} \le k \le \frac{N}{2} - 1 \qquad (18)$$

Because of the orthogonality relation

$$\frac{1}{N} \sum_{j=0}^{N-1} e^{ipx_j} = \begin{cases} 1 & \text{if } p = Nm, m \in \mathbb{Z} \\ 0 & \text{otherwise} \end{cases} \qquad (19)$$

we have the inversion formula

$$u(x_j) = \sum_{k=0}^{N-1} \tilde{u}_k e^{-ikx_j}, \quad j \in \mathbb{N} \cap [0, N-1] \qquad (20)$$

So, the polynomial

$$I_N u(x) = \sum_{k=0}^{N-1} \tilde{u}_k e^{-ikx_j} \qquad (21)$$

is the N-degree trigonometric interpolant of $u$ at the collocation points.

## Fast Fourier Transformation

Transforming the internal energy terms should therefore not be a problem. The image energy, on the other hand is a nonlinear function. Transforming such a function and integrating it into the resulting differential equation in the Fourier space may cause some computational problems. An alternative algorithm was therefore developed. Each iteration of the snake evolution equation was split into two steps.

The internal energy term was iterated in the Fourier space, transformed back to the physical space where the image energy term was added.

We start by transforming the Equation (9) into Fourier domain

$$F\left[\frac{\partial v}{\partial t}\right] = w_1 F\left[\frac{\partial^2 v}{\partial s^2}\right] + w_2 F\left[\frac{\partial^4 v}{\partial s^4}\right] \qquad (22)$$

which gives

$$\begin{aligned} \frac{\partial \tilde{v}(\omega, t)}{\partial t} &= w_1(i\omega)^2 \tilde{v}(\omega, t) - w_2(i\omega)^4 \tilde{v}(\omega, t) \\ \frac{\partial \tilde{v}(\omega, t)}{\partial t} &= -w_1\omega^2 \tilde{v}(\omega, t) - w_2\omega^4 \tilde{v}(\omega, t) \\ \frac{\partial \tilde{v}(\omega, t)}{\partial t} &= -(w_1\omega^2 + w_2\omega^4)\tilde{v}(\omega, t) \end{aligned} \qquad (23)$$

This is a first-order differential equation with respect to $t$ and its general solution is

$$\tilde{v}(\omega, t) = F[v(s, 0)]e^{-i(w_1\omega^2 + w_2\omega^4)t} \qquad (24)$$

where $F[v(s, 0)]$ is evaluated using FFT.

We iterate this explicit solution with a time step and after each step transforms it back to the physical space where the image energy term is added. In the kth step, the iteration equation is

$$v(s, t_k) = F^{-1}[\tilde{v}(\omega, t_k)] - \nabla E_{\text{image}}(v(s, t_k)) \qquad (25)$$

## Results and Discussion

The experiments with synthetic and preprocessed grout captured images are performed using Matlab® based functions for its discrete implementation. The active contour was initialised on the image extracted from captured video in case of grout injection and generated by user in the other cases. The initialization was followed by temporal and spatial snake development using an iterative procedure. With these two examples (Figure 2A–H), we can examine the effect of the image force $\nabla P$. The direction of this force implies the steepest descent in $P$, which is natural because we want to get a minimum of $P$. Equilibrium is achieved when $P$ is minimum in the direction normal to the curve. Although the initial estimated contour is close to an edge, instability can occur because of the discrimination of the evolution problem (Figure 2E–H). All the tests were performed using a PC Pentium III 866 MHz. We used $w_1 = 0.6$, $w_2 = 0.0$ and $w_\sigma = 0.0$ for these snakes. For the examples presented in Figure 2, time consumption is
<figure><img src="image_2.png" /><figcaption>**Figure 2A–H:** Synthetic images (64 × 64 pixel) proposed by Xu and Prince [14] to test the algorithm</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 3A–C:** Deformable models to extract object boundaries from SEM images</figcaption></figure>

less than 60 s. Fortunately, this numerical phenomenon does not occur when we realise a grout injection. To extend this application to microscopy analysis, we test an SEM image of grouted sand.

The initial snake contour is shown in Figure 3A. The Figure 3B shows a developed contour at the next level after 50 iterations. It illustrates that the gradient forces ‘pulled’ the snake towards the pore. Figure 3C gives the final result of the snake, fully enclosing the pore after 100 iterations.

The active contour performance is still dependent on the local maxima in the image, if we only use this feature as potential function. The high gradient local values yielded from reflection, although smoothed from the filtering, still have important contribution as a local potential energy maximum. The performance of the active contour depends on the image preprocessing, and edge map [15] or gradient vector flow (GVF) [14] computation. Besides, small contours which are not previously expected appear because of singularities in the gradient intensity field captured by the spectral active contour. To solve Equation (13), numerical explicit schemes are often employed. This type of scheme prevents from matrix inversion (great time consumption). Moreover, it is sensitive to the edge map singularities.

To reduce time consumption because of edge map and/or GVF field, we implemented FFT formulation to transform 2D computation to 1D computation.

This latter procedure allows converging four times faster than the initial algorithm presented in Ref. [16]. Indeed, by reducing the data storage dimension of the system to solve, we increase the convergence rate. This can also be done by users having knowledge of image processing analysis to set relevant parameters of $w_1$ and $w_2$. However, it does not exist fast techniques to evaluate them. In order to get a first estimation of these parameters, we can employ an ad hoc algorithm [1] or a decomposition of the solution in the hermite finite element basis [17]. Thus, for our SEM images, the time consumption varies from 45 min for initial algorithm to 12 min using the FFT procedure. We analyse the 1D grout injection process to target different possible shapes of the grout front. The FFT active contour improves the performance of the ordinary contour by minimising memory consumption. This technique is relevant when the size of the analysed images is significant (more than 256 × 256 pixels).

Figure 4(A–C) presents a first set of frames acquired during the water saturation of a dry Loire river sand column. It constitutes the reference to compare slightly the algorithm performances. This case is equivalent to a threshold segmentation which can be reduced to a binary image. Then it is easy to extract front and compare it.

The contour detected by the FFT active contour model is shown in Figure 4(D–F) for water saturation
<figure><img src="image_4.png" /><figcaption>**Figure 4:** (A–C) Initial frame from water saturation; (D–F) contour from water saturation; (G) transition zones; (H–L) contour from grout injection</figcaption></figure>

<table><caption>**Table 1:** Parameters of different grout injections of Loire river sand columns</caption><thead><tr><th>Soil column</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>d (m)</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>Q (m³ s⁻¹)</td><td>5.6e-6</td><td>5.6e-6</td><td>5.6e-6</td></tr><tr><td>ν<sub>f</sub> (m s⁻¹)</td><td>2.3e-3</td><td>2.8e-3</td><td>2.8e-3</td></tr><tr><td>D (m² s⁻¹)</td><td>7e-6</td><td>7e-6</td><td>1e-5</td></tr><tr><td>n</td><td>0.32</td><td>0.32</td><td>0.303</td></tr><tr><td>C/W</td><td>0.42</td><td>0.42</td><td>0.36</td></tr><tr><td>C₀ (kg m⁻³)</td><td>1410</td><td>1410</td><td>1380</td></tr></tbody></table>

process. Figure 4G shows an entire cycle of grout injection. Figure 4(H–I) presents relative grout position during the test. In all cases, the contour encloses the correct front position. Less than 10 iterations at each step are needed to converge towards the final shape. This type of algorithm uses spectral method and it is well recommended to close contour. Fortunately, this technique still gives best results in grout injection image analysis. Results show that active contour model and in particular parametric deformable model can be easily employed to track grout propagation in laboratory columns. However, other kinds of active contour models may be used to realise image segmentation and motion tracking.

Table 1 summarises the column diameter (d), sample porosity (n), injection flow rate (Q), cement to water ratio (C/W), and initial grout concentration (C₀) used for each test. It also shows the superficial velocity (νf) and the dispersion coefficient (D) obtained by CDE formulation.

<figure><img src="image_5.png" /><figcaption>**Figure 5:** Typical optical breakthrough curve for Loire river sand</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6:** Relative grout concentration estimated (using CDE) from BTC measured in Loire river sand column 1 for different instants (50, 100 and 150 s)</figcaption></figure>

Figures 5–8 present relative concentration C/C₀ measured by the Beer–Lambert theory (data) and calculated by CDE formulation (CDE) with respect to the distance x to the injection source point. Figure 5 shows that grout concentration can be well estimated in the injected zone, i.e. before the sharp concentration variation. It is coherent with classical concentration versus time breakthrough curves (BTC) commonly observed in porous media.

However, because of image noise effect caused by intensity illumination variation along the column, removing this noise which occurs in the area where the column is only saturated by water is not simple without altering the linear CCD response. Fortunately, grout concentrations in this zone are known and equal 0. So this artefact does not disturb the region of interest.

After all these considerations, some results obtained from three soil columns are presented in Figures 5–8.
<figure><img src="image_7.png" /><figcaption>**Figure 7:** Relative grout concentration estimated (using CDE) from BTC measured in Loire river sand column 2 for different instants (15, 120 and 275 s)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8:** Relative grout concentration estimated (using CDE) from BTC measured in Loire river sand column 3 for different instants (20, 70, 140 and 250 s)</figcaption></figure>

Samples 1 and 2 have the same porosity *n*, same cement to water ratio *C*/*W* and same flow rate *Q*. Sample 3 has the same characteristics except its porosity. After optical calibration of CCD camera we show identical BTC curves and identical noise effect.

Besides, we reach the limit on the CCD camera sensitivity. Close to these intensity rates, we may have deviation between Beer-Lambert formulation and video-captured data. We can also notice a similar noise effect on various experiments. Thus, further image analysis of this phenomenon may allow attenuating this offset value. CDE formulation estimated nearly equivalent superficial velocities and dispersion coefficients. Besides, superficial velocities are close to the commonly used value $4Q/(\phi \cdot \pi \cdot d^2)$.

A part of the BTC curve fits the ADE model curve and reveals that the superficial velocity and the

dispersion coefficient are constant during the first part of the injection, i.e. until grout does not reach the outlet of the column. Furthermore, these figures show that the cement particle concentration decreases before sharp concentration variation, which is linked to the transition zone, i.e. area highly affected by the miscibility of the grout. It is characteristic of particle retention. Incorporating the filtration mechanism in the CDE formulation may allow capturing this filtration phenomenon.

## Conclusions

Column and field systems are usually treated separately because the column has finite length and identifiable outlet while the field usually extends to infinite length in all directions. Using a simple numerical model to predict breakthrough curves and optical concentration measurements, we were able to test the effects of different conditions on the predicted grout propagation. The use of the generally accepted CDE formulation allowed us to show that we can estimate the superficial velocity and the dispersion coefficient during grout propagation.

It was also demonstrated that convection and dispersion mechanisms must be completed by a filtration mechanism. This explains why we observed that the component concentration decreased before the transition zone because of the miscibility of the grout. It is concluded that the CDE formulation can be used to estimate the dispersion coefficient. Optical concentration measurements using the Beer-Lambert theory can be used within their region of validity near the transition zone.

## REFERENCES

1. Kass, M., Witkin, A. and Terzopoulos, D. (1988) Snakes: active contour models. *Int. J. Computer Vision* **1**, 321-331.

2. Malik, J. and Perona, P. (1987) A scale space and edge detection using anisotropic diffusion. *IEEE Comp. Soc. Workshop Comp. Vision*, Miami, FL, USA.

3. Bear, J. (1972) *Dynamics or Fluids in Porous Media*. American Elsevier Publ., Co., New York.

4. Ogata, A. and Banks, R. B. (1961) A solution of the differential equation of longitudinal dispersion in porous media. U.S. Geol. Surv. Prof. Pap. No. 411A.

5. Harvey, R. W. and Garabedian, S. P. (1991) Use of colloid filtration theory in modeling movement of bacteria through a contaminated sandy aquifer. *Env. Sci. Technol.* **25**, 178-185.

6. Pachepsky, Y., Benson, D. and Rawls, W. (2000) Simulating scale-dependent solute transport in soils with the fractional advective-dispersive equation. *Soil Sci. Soc. Am. J.* **64**, 1234-1243.
7. Terzopoulos, D. (1986) Regularization of inverse problems involving discontinuities. *IEEE Trans. PAMI* **8**, 4.

8. Thikonov, A. and Arsenine, V. (1974) *Méthodes de résolution de problèmes mal posés*. Mir, Moscou.

9. Courant, R. and Hilbert, D. (1953) *Methods of Mathematical Physics*, vol. 1. Interscience, New York.

10. Cohen, L. D. (1991) On active contour models and balloons. *CVGIP: Imag. Under.* **53**, 211–218.

11. Cohen, I., Cohen, L. D. and Ayache, N. (1992) Using deformable surfaces to segment 3-D images and infer differential structures. *CVGIP: Imag. Under.* **56**, 242–263.

12. Canuto, C. Hussaini, M., Quarteroni, A. and Zang, T. (1988) *Spectral Methods in Fluid Dynamics*. Springer-Verlag, Berlin, Heidelberg.

13. Boyd, J. P. (1985) *Chebyshev and Fourier Spectral Methods*. Springer-Verlag, Berlin, Heidelberg.

14. Xu, C. and Prince, J. L. (1998) Snakes, shapes, and gradient vector flow. *IEEE Trans. Image Processing* **7**, 359–369.

15. Deriche, R. (1987) Using Canny's criteria to derive recursively implemented optimal edge detector. *Int. J. Computer Vision* **1**, 167–187.

16. Ait Alaïwa, A., Saiyouri, N. and Hicher, P.-Y. (2006) Experimental detection of one dimensional grout spread using image processing active contour model, 2nd Int. Conf. Problematic Soils, Selangor, Malaysia.

17. Cohen, L. D. and Cohen, I. (1993) Finite-element methods for active contour models and ballons for 2-d and 3-d images. *IEEE Trans. PAMI* **15** 1131–1147.