1NUMERICAL ANALYSIS OF NON-LINEAR 
LARGE-STRAIN CONSOLIDATION AND FILLING 
A. HUERTA and A. RODRIGUEZ 
Departamento de Matematica Aplicada m, Escuela Tecnica Superior de Ingenieros de Caminos, 
Canales y Puertos, Universitat Politecnica de Catalunya, Jordi Girona Salgado 31, 
08034 Barcelona, Spain 
Abstract-Finite strain consolidation and filling of soft sediments at high water level is a challenging 
problem because of its highly non-linear physical and mathematical aspects. Several numerical schemes 
designed for this problem are presented as well as simple numerical improvements for a better handling 
·.of the extremely high variations of the material properties with depth . The numerical algorithms developed 
., are robust and verify convergence of the iterative schemes instead of the more classical approaches based 
t,-00 choosing time increments 'sufficiently' small and assuming convergence at every step. 
· · · A set of computer programs has been developed to predict magnitude and rate oflarge-strain self-weight 
:. ·one-dimensional and pseudo bi-dimensional (i.e. one-dimensional deformation, bi-dimensional flux) 
Jt·consolidation during and after deposition, that is, coupling filling and consolidation phenomena. The 
".·actual life of the deposit can be numerically simulated combining filling periods and quiescent periods 
',_ where surcharges ( or capping) can exist. Consequently, they are a basic technique for the design of disposal 
' ponds. 
INTRODUCTION 
"'rediction of densification magnitudes and rates 
ajor concern in the disposal of fine grained 
H mineral wastes, such as mine tailings, 
a red mud, oil sand sludge, dredged materials, 
. 'e primary issues center around predicting the 
ependent capacity of a given disposal area-as 
s.ultimate storage capacity and useful life-and 
· e-rate improvement of material properties for 
·tion purposes. Intrinsic to both issues is the 
of consolidation and the associated material 
. .. · of computer programs has been developed 
i-e'dict magnitude and rate of large-strain self-
) consolidation during and after deposition, 
is; coupling filling and consolidation phenom-
T:_he actual life of the deposit can be numerically 
... ~ted combining filling periods and quiescent 
ods. 
· .~e to the loose nature of the sedimented )it (solid contents ranging from less than 10% 
p· 30% ), the ensuing vertical settlements are 
,Ilrextremely large and beyond the range that can 
~,l!dled by classical small-strain consolidation 
. ?~s; accordingly non-linear finite deformation 
~1.s must be developed. Since the physical and 
. _h~tnatical aspects of this problem are highly 
~g-1_near, the numerical algorithms developed 
t1~Y. convergence of the iterative schemes instead 
. .. more classical approaches based on choosing 
increments 'sufficiently' small and assuming 
r~_ence at every step. This procedure precludes 
,,,);:.~real aberrations that appear with other 
o·aches and increases accuracy in the first half 
of the consolidation process. It should be noted that 
due to the alternation between periods of filling and 
quiescent consolidation good accuracy at the first or 
middle stages of consolidation is needed, because the 
results obtained after a cycle are the initial conditions 
for the following one. This is not the case of classical 
consolidation analyses concerned mostly with the · 
final stages of consolidation. 
The next section contains a brief background 
of finite-strain theory and a discussion of the 
material relationships employed. This is · followed 
by a statement of the numerical problem together 
with the pertinent algorithms studied and developed. 
Next a couple of one-dimensional examples and 
the bi-dimensional extension of the problem are 
presented. 
THEORETICAL CONSIDERATIONS 
The theory presented here is patterned after 
the early work by Gibson et al. [l],. and several 
modifications that are described in (2-7] among 
others. Two standard formulations for the consoli-
dation equation are employed in finite-strain theory, 
depending upon whether the void ratio, e, see for 
instance [I, 2, 4, 5, 8), or the excess pore-water press-
ure, u, see (3, 9], is selected as the dependent variable. 
The approach adopted here uses excess over hydro-
static pore-water pressure, u, as the dependent vari-
able. The one-dimensional consolidation equation 
obtained is [3]: 
(1) 
Huerta, A. and Rodríguez-Ferran, A., Numerical Analysis of Nonlinear Large-Strain 
Consolidation and Filling, Computers and Structures, Vol. 44, Issue 5, pp. 905-919, 1992
2in which the material coordinate z is given by 
r· da' 
z(a) = Jo I+ e
0
(a ')' (2) 
where a is the initial vertical (Lagrangian) coordinate, 
e0 is the initial void ratio, k is the coefficient of 
permeability, y •. is the unit weight of water, t is time, 
and u' and <Jb are the effective stress and buoyant 
stress·; respectively. The material coordinate, z, can 
be interpreted as the height of solid particles between 
the planes determined by a = 0 and a, which is 
constant with time. 
The buoyant stress is defined as the total stress 
minus the hydrostatic pressure, that is 
<Jb(z, t) = (y, - y .. )(z,op(t)- z) + q,un (3) 
where y, is the unit weight of the solids, q,ur is the 
surcharge or capping load, and z10P is the material 
coordinate corresponding to the top layer of the 
deposit. Since the surcharge, q,ur, is usually assumed 
to be deposited instantaneously, the buoyant stress 
depends only on time through z10P. That is, the 
variations of the buoyant stress are directly depen-
dent on the increase of the total height of solid 
particles, i.e. any additional deposition of solid par-
ticles, namely the filling rate. Thus it is precisely the 
term oub/ot, in eqn (1), which differentiates periods 
of filling and quiescent consolidation. 
- The resolution of eqn (I) requires the speci-
fication of boundary and initial conditions. Three 
different boundary conditions have been con-
sidered: drained, undrained or prescribed piezometric 
head at any boundary . The initial condition is 
the distribution of excess pore-water pressure with 
depth. 
Appropriate relationships between void ratio and 
effective stress . and between the coefficient of per-
meability and the void ratio must be implemented. 
The · computer programs developed include several 
specific types of equations proposed for these two 
relationships, see [6, 10). Available test data on fine 
grained materials at low solid contents suggest highly 
non-linear e-u' and-k-e behaviors which are usually 
written as power Jaws. 
One important drawback of the usual compress-
ibility relations is that zero effective stn,ss implies 
infinite void ratio. To avoid this problem without 
changing the compressibility relationship and the fact 
that at the drained surface the buoyant stress and the 
excess pore-water pressure are zero, the effective 
stress is defined as: 
(4) 
where u0 is obtained by introducing in th~_compress-
ibility relation the initial void ratio, e0 • The validity 
of this assumption is discussed in (11) for steady-
state seepage-induced consolidation; moreover, a 
sensititivity analysis conducted with the developed' 
programs concluded that the influence of llo in 
settlement is negligible-· and stays under 5% for-the 
material property distributions, as advanced in [I I]. 
Physically, this should be expected since the initial 
void ratio which arbitrarily represents the lirnir 
between sedimentation and consolidation is usually 
very large [IO] and this implies a small u0. Other 
definitions, such as, 
u' = <Jo 
induce unrealistic situations in the upper part of the· 
deposit; namely, constant distributions of void rati~ 
and permeability, thus no consolidation. 
It should be noted that sand mixing is easily 
introduced by assuming that the clay phase governs: 
both the mechanical and flux behavior of the mixture, 
That is, in the material relationships a clay void ratio' 
(volume of voids divided by volume of clay) is use_d 
instead of the classical void ratio (volume of voirlt 
divided by volume of solids). 
Finally, it should be noted that the coefficien(s'! 
k/(1 +e) and de/du' in eqn (I) depend ultimatelyoWi 
the excess pore-water pressure, u. Both perrneabili~ 
and compressibility have such large variations durinf 
the consolidation process, that these coefficient{ 
cannot be considered constants, or even 'almost'; 
constants. Consequently, the parabolic eqn (lfi il 
clearly highly non-linear. 
NUMERICAL SCHEMES 
A fully implict finite difference method was chosert 
for the solution of eqn (1) because of its stabilit 
However, most of the conclusions observed are ind~F 
pendent of the particular method employed and ~j 
be generalized to other techniques such as the finite 
element method . 
The dimensionless form of eqn (1) is given by 
p[ou - O<Jb  = !... [q OU J. (6,J 
or ot oz oz 
where the same symbols are used for dimensionle~r 
variables in order to simplify the expressions. T~t 
coefficients p and q in eqn (6) can be interpreted_f"" 
the dimensionless compressibility and perrneabih_tr,7 
therefore, they are functions of u. _ ., 
After discretization, the resulting system of eqU~'1 
tions is non-linear and given by 
fl 
where u" + 1 and u" are the unknown nodal veeto d 
of excess pore-water pressure at time tn+·' an 
t", respectively, A is the tridiagonal matrix wit~~· 
coefficients are function of p and q (this functto 
3z Initial 
filling 
Quiescent 
consolidation 
Second 
filling 
T 
Ll.T1 Ll.T2 Ll.T3 
Fig. 1. Domain discretization ( one node added per time-step). 
"nds· on the discretization scheme implemented), 
fisthe nodal vector of increments in the buoyant 
·•· for filling periods (f = 0 in the quiescent con-
!ltion periods). 
f:;t, the domain discretization is discussed 
use· of the difficulties encountered during the 
gperiods particularly in the upper areas of the 
oslt Then, several approaches for solving eqn (7) 
:,presented. Due to the highly non-linear behavior 
the problem, finding the solution of eqn (7) is not 
'reasy task. 
];)omain discretization 
Finite differences are employed; consequently, a 
discretization in both independent variables, time (t) 
lij!d space (z ), is required. During filling periods 
new, material is deposited and, since Lagrangian 
coordinates are used (the spatial mesh follows the 
particles), new nodes are introduced in the spatial 
discretization. 
Figure I shows a classical [3] domain discretization 
where· one node is added per time-step during filling 
periods. This is a usual and simple technique where 
the filling rate is the spatial increment divided by the 
time-step (recall that z is the height of solid particles). 
However, the strong dependency between both incre-
ments is sometimes a serious drawback. 
z • Real new nodes 
o Ficticious nodes 
T 
Fi~. 2. Multiple node technique per 'deposited' layer. 
Self-weight consolidation problcu1s present ex-
tremely high variations of permeability and com-
pressibility with depth in the upper area of the 
deposit. This requires a very fine spatial mesh able to 
capture the changes in material properties; recall that 
linear interpolation (two nodes per layer) is employed 
in every 'deposited' layer. The ensuing time-steps are 
usually very small and the computer cost clearly 
uneconomical. 
If several nodes are introduced per 'deposited ' 
layer (time-step), as shown in Fig. 2, the high vari-
ations of compressibility and permeability are better 
captured and both discretizations (t and z) are inde-
pendent of each other. To illustrate it, Fig. 3 shows 
the void ratio distributions after the first time-step 
of an initial filling period over an undrained bed. 
As expected, the five nodes per time increment com-
pares much better than the one node technique with 
an 'exact' (extremely fine mesh) solution. 
The area under the void ratio distribution is 
directly related to the depth of the deposit. The 
over-estimation of e with the one node per time-step 
discretization induces an under-estimation of settle-
ments, while the multiple node technique developed 
compensates the areas and consequently predicts 
settlements much more accurately. This is an 
N 
0 
/ 
I 
I 
I 
I 
I 
/ 
/ 
/ 
/ 
/ 
/ 
/ 
I 
/ ---- I node for 6 T 
I 
/ -- 5 nodes for Ll.T 
· / -·-· Very fine mesh // 
. / 
5 10 15 20 25 30 
Void ratio ec 
Fig. 3. Comparison between one node and multiple node 
techniques. 
4important feature that will be clearly observed in the 
presented examples. 
Linearization of the-problem 
In order to solve the non-linear set ,of algebraic 
eqns (7), the classical approaches [3] linearize the 
problem replacing in the matrix A the solution of the 
last time increment, i.e. they assume A(u•+ 1) = A(un). 
Physically, this means that the material properties lag 
one time increment behind the solution; recall that 
A iS"·a function of the-dimensionless-·compressibility 
and permeability. This interpretation advances the 
problems to be found , since high variat ions of 
material properties occur in the first stages of con-
solidation. Equation (7) is then reduced to a linear 
system of algebraic equations 
A(u")u"+ 1 = u" + f (8) 
which is easily solved at each time increment. 
. The .main advantage of this technique is that its 
implementation is extremely simple, but it does, 
however, present serious drawbacks. The basic hy-
potheses are that the time increment is sufficiently 
small and that small increments of time induce negli-
gible modifications in A. Therefore, only under these 
circumstances can both matrices A(u" + 1) and A(u") 
be assumed · ·approximately equal. Obviously, this 
method may need uneconomical time-steps. More-
over, it should be emphasized that with the lineariz-
ation , convergence of the solution is not imposed. 
Thus, this method requires experienced operators 
who know from the final results (settlements vs time, 
void ratio or pore-water variations with depth, etc.) 
when the time increment is small enough and conver-
gence was attained. In conclusion, the linearization 
of the problem is, in spite of its computational 
simplicity, a dangerous technique for solving eqn (7). 
Iterative techniques 
In order to solve the non-linear set of equations, 
standard ·non-linear numerical techniques should be 
implemented. There are a large number of methods 
for the resolution ofeqn (7) which are based on the 
general scheme defined by 
u;;t l = u;;+ i - [C(u;;+ i )J-1 
where the subscript k denotes the iteration counter 
and C is an arbitrarily chosen matrix. Convergence 
of the scheme is assumed when the difference between 
two successive approximations of un+ 1 and the resi-
dual vector, [A(u;; + 1 )uZ + 1 - (DZ+ f)J, are smaller than 
a predetermined tolerance. 
In fact, the choice of C classifies the numerical 
technique employed.and induces .the. order .of conver-
gence of the method. Therefore, the choice of C is 
important and must be done carefully depending on 
the computational cost (storage and CPU time). 
Most of the usual techniques take C as the Jacobian 
matrix, J, or an approximation to it. The Jacobian 
matrix for this problem may be written as 
(10) 
The first method implemented consists in the 
following approximation of the Jacobian matrix 
C = A. This method of order one is in fact. th; 
well-known 'fixed-point iterative method', which 
may be written as 
(Ill 
and was chosen because its implementation in .a 
linearized computer code is extremely simple. The. 
computer cost per time-step increases because egn 
(11) must be constructed and solved, at every instant 
t" + 1, for each iteration, (k + I), up to convergence . 
However, the overall computer cost decreases com-
paratively with the linearized technique because the 
time-step may be larger. 
In some problems, mostly during filling period$: 
and because of the highly non-linear behavior, con-
vergence is too slow or never obtained. This result. 
may seem negative but it actually gives valuable 
information: the time-step is too large for the pre-· 
scribed tolerance. Obviously, the same time-step fo.,a 
linearized code would give incorrect results. 
In order to solve this problem several improve· 
ments where tried: Aitken acceleration, time-st~.P 
splitting, etc. None of them was satisfactory enougb, 
the main reason being that any acceleration technique 
for convergence was not sufficient because of the 
extreme non-linear behavior of the problem. 
Consequently, a full Newton- Raphson method 
is used to avoid the previously cited difficulties fan 
some of the studied problems. This is a second order 
method where the iteration matrix , C, is taken eqµa.1 
to the Jacobian matrix, J . As shown in eqn (10),Jlfe 
derivatives of the coefficients of A with respect to l!Zt 
must be computed to evaluate J. The dependence· 
of these coefficients in u is as follows: 
A=A(p , q) (12@ 
. P.= p(u ') oJ@ 
q=q(e,k) (!Zi 
k=k(e) (I~ 
e =e(u ') (Ji 
u' = u '(u). 02n 
A modifi~~tion of the compu~er ~ode is requi~ , 
a modular and structured orgamzat10n must t,e.u 
5. sed on the above presented equations. That is, 
~!,, q. uation in (12) represents a necessary module 
everye .. dh'd '. fi • ' i re -rhe ·desff,ed ..funct1ons .an t eir envat1ves..are 
rJ3 ted At the end, the coefficients of A and the 
,.·.JllpU . W bian matrix are obtained. In this manner, modifi-
Jaco d 1·1 · 1 
,\ •.. s of the computer co e are re at1ve y s1mp e. 
eatJOil . 
·;; -;-.·nstance, if the numerical scheme changes , only 
i;Of I . d (12 ) d . d . . tfiiJllodule ass?c1~te to eqn . ~ an 1t~ en~at1ves 
,'.,: h nged; while, 1f the const1tut1ve relat10nsh1ps are 
!SC a 2 ) d h . d . . iii~dified, eqns (12d) and (1 e an t e1r envatives 
/ ' t-be updated without perturbation of the general f u\t of ,the code. Knowledge of the analytical 
·JY;essions of eqn (12) and their derivatives renders 
.exp · I d 'hth fuU Newton-Raphson optlffia compare wit o er 
,fil h order convergence schemes that do not evaluate 
t ·!xplicitly, such as the quasi-Newton methods. . 
-···This second order scheme reaches convergence m 
£few iterations and for very small tolerances, on the 
'Jfder of 10-6. Obviously, the computer cost per 
m;ration is larger because of the evaluation of the 
dedvatives. In the linear or first order schemes, 
th~ construction of A takes approximately 70% of 
.,(he CPU time while its resolution only uses the 30% 
7f~ft;' here, two matrices are computed in every iter-
"dtion, thus the cost per iteration almost doubles. 
However, in general it is less expensive than other 
t~ hniques because of the important reduction in the 
'Aumber of iterations needed for convergence (three or 
four iterations during the first stage of consolidation 
'f6r tolerances of 10- 4), It should be noted that, as 
·expected, the solutions obtained after convergence 
3yith any method are identical. This proves the 
Wfisistency of the techniques employed and the 
,lth.iqueness of the solution . 
·Finally, a simple but illustrative example of 
Ufo extreme non-linear behavior of eqn (7) and the 
J#iportance of adequate initial approximation is 
·~hown. Similarly as before, the first time increment of 
;1::. initial filling analysis over an undrained bed is 
~thdied. The top surface is drained and consequently 
if one node is added per time-step, there is only one 
_unlcnown: the excess pore-water pressure, u, at the 
:b_ottom. The system of non-linear eqns (7) is now 
feduced to one equation and the residue can be 
plotted vs u, see Fig. 4. As shown in the figure, only 
one solution exists for the physical range of u. 
tiowever, it can also be observed that considerable 
difficulties in the resolution of this problem are 
J~duced by a second solution exterior to the domain 
ilnd the shape of the curve. The developed technique 
wli!ch uses as an initial approximation the solution in 
t,he previous step, does the first iteration with the 
?rder one scheme and then uses the second order 
method, converges in a few steps (less than five for 
tolerances under 10- 6) for all the studied cases. 
ONE-DIMENSIONAL EXAJ'1PLES 
To show the applicability of the developed com-
puter codes and the advantages of the non-linear 
Excess pore - water pressure 
Fig. 4. Residue vs excess pore-water pressure for a one 
degree of freedom filling analysis. 
numerical schemes, two illustrative examples are 
presented. ·' 
The first example is the classical [12] self-weight 
consolidation of a soil column, starting from uni-
form conditions (i.e. constant void ratio with depth). 
This condition is never actually encountered, but 
it represents laboratory tank tests where filling 
is 'instantaneous' compared with pore pressure 
dissipation. 
The constitutive relations are 
e = 4.674(u ')- 0·22 
(e )S.511 
k = 3.358 x 10-1 -1-· 
+e 
(13) 
where u' is in tonnes per square meter (102 kPa) and 
k in meters per day. The operating parameters in 
eqn (13) are those of one test case considered in [12]. 
The spatial mesh is composed of 101 nodes and 
consolidation is simulated over a year. 
Figure 5 shows a . comparison . between .. con-
solidation height-time curves. Three analyses are 
presented: linear model with 20 time-steps (discon-
tinuous line); iterative model with 20 time-steps 
(continuous line); and iterative model with 1000 
---- Linearization method 
-- Iterative method 
6 
-- Exact SOlution at initial 
instants 
----
2 
0 0 . 1 0 . .8 09 1.0 
Fig. 5. Height of the deposit vs time. Quiescent consoli-
dation . 
6N 
c: 
0 
l 
.3 
0 
I 
I 
I 
--- Linearization method I 
-- Iterative technique I 
,,.~· 
,..., 
/ 
/ 
/ 
I 
5 15 
Void ratio ec 
I 
I 
I 
I 
I 
,/ 
_,... 
20 25 
Fig. 6. Comparison between linear and non-linear void ratio 
distributions after 3 months of quiescent consolidation. 
time-steps (continuous line). The linear model clearly 
underestimates settlement and has an unrealistic 
change in curvature, while the non-linear model 
gives similar results to those of a 1000 time-steps 
imafysis. The straight line which represents the 
theoretical exact initial behavior (constant consoli-
dation rate, see [8]) is in good agreement with both 
non-linear solutions in spite of the large time incre-
ments taken here (20 time-steps per year implies an 
increment almost 350 times larger than the one used 
in [8]). 
A comparison between the void ratio distribution 
obtained with both techniques after 3 months is 
illustrated in Fig. 6. The linear solution over-
estimates the void ratio in the upper area where the 
initial conditions still persist. One iteration, with 
the implemented time-step, is not enough to model 
the important changes in the void ratio. 
It should be noted that, as expected, both method-
ologies converge to the final conditions which can be 
evaluated analytically. However, the next example 
shows that due to the alternation between periods of 
filling and quiescent consolidation, good accuracy at 
the first or middle stages of consolidation is needed, 
because the results obtained after a cycle are the 
initial conditions for the following one . 
By proper manipulation, virtually any sequence of 
days. It is clearly observed that the linear approach 
overestimates the needed depth of the deposit , while 
the non-linear scheme with five nodes per time-step 
converges clearly to the dash-dotted line computed 
with a non-linear scheme and time-step of I day. 
It is obvious that the computer cost per time-step 
is larger in the non-linear codes. It would therefore 
seem adequate to compare linear and non-linear 
results using different time-steps so that both pro-
cedures require similar amount of computer time. 
However, predicting the proper time-step reduction 
at every instant for the linear analysis is not an easy 
problem and may be as costly as solving the non-
linear problem itself. This is shown by the fact that 
the number of iterations varies with time, for a given 
time increment, in a non-linear simulation . 
PSEUDO BI-DIMENSIONAL EXTENSION 
One-dimensional large-strain consolidation is an 
acceptable model for the overall behavior of a usual 
waste slurry containment. However , when consoli~ 
dation is accelerated by vertical drains, or when 
consolidation occurs in narrow trenches, bi-dimen-
sional effects clearly influence the final solution,. 
Similarly to classical small-strain consolidation, an 
extension of eqn (7) to bi-dimensional flux a11:d, 
400000~---------------------------, 
0 
0 2 3 4 5 
Time ( years ) 
(a) 
30~---------------------i 
filling and quiescent settling, with or without sur- 20 , 
charge, may be simulated. To illustrate the versatility ] /\ 
,\ 
I' I 
I 
I 
of the computer programs, and their applicability to t 1s I \ 
the prediction of the containment requirements and .; / '",, 
useful life, a problem composed of four filling and J: 10 1 ' 
· I -----Linearized, t:,.T • 10 days 
qutescent consolidation phases is shown . The filling I\ I -Iterative, "T. 1o days 
periods span over half a year , while the quiescent 
5 
1
/ ~---- - ·- ·-Iterative,~T· I day 
consolidation takes a year. Three different filling rates L--...l.----.J. 1----+--~_._~-~ E> h O 2 3 4 5 
ave been used: 100 x 106 kg/year (initial filling), Time (years) 
3_00 x 106 kg/year (second and third cycles) and (b) 
200 _x 106 kg/year (last cycle); (see Fig. 7a). · .. . . . . ,cons0li:. 
Figure 7 presents the height variations with time Fig. 7. (a) Sequence of fillmg periods and quiescent es. dation phases. (b) Comparison between different s~heJllof 
for different numerical schemes. It compares the during the simulation of the prediction of the containJlle 
linear and non-linear schemes with a time-step of 10 requirements and useful life. 
70.1 0.2 0 .3 0 .4 0.5 0 .6 0 7 0 .8 0.9 1.0, _ 1.00 0 . 1 0 .2 0 .3 0.4 0 .5 .0.6· 0.7 0.8 
I 
~ 0 . 1 0 . '==10 .~ 0.9 _--,0 . 1 0 . 1 0 .9 
,. --- ~ -
i0.7~l'/~o, o?Jo, :: r~: :: :: 
0, 06~1/;,,, . . ~ A 06 0.6 1(9 ~ 0.6 
h. 
1i o.5 ~~ 0 _5 o.5 o.511111/ / J' . ,/"" o.5 
o '<I- o!'> i 0.4 / ~ -6 -10 .4 0 .4UIIIIYf I /' n!'>_. -10 .4 
'' ''~[11f  ~"' 
. 
0
'.~~ I I:'. ::_j? B,, :'. 
o . 2 o .3 o .4 o.5 ·o .6 o . 1 o.8 0.0 1.c? o ·0 .1 0 .2 · o .3 o .4 0 .11 o .6 o. 1 o .8 0 .0 1.0 
: .:
0
~:;=• 1' '( 
0
~' 'i' j;:~ 
' --0 .2 -0 .7 
O' ---. o.Z __-
____. _ --Jo .s 
~ -3 
::i II! f/~ ---io _ ~= 
::11111~ f/ ,L tCl 
0 0 .1 0.2 0 .3 0 .4 0 .5 0.6 0.7 10.8 I0.91 1.0 
t = 3 months t = 6 months t = 12 months 
,,o ,, ,, ,, o, ,., ,, o, ,., o.s ",, ,,~ ~ o_c:..·-' ."""-o"°'.2 _=--0 .~-~ o-.4 __ 0::-.~~ o~._s ~o_ .1 __ 0_.8 
- - - r.,c- ~ ~-=--=--==- ----==-
= 09 09 
. "// 
0 .8 0.8 ~ 0.8 
0.7 
1.00 0 . 1 0.2 0 .3 0 .4 0.5 0.6 0 .7 0 .8 0.9 1.01.0 
----------------=I0 .9 
9
~9-,.____,9 
0 .8 
0 .7 
~05~\\\  ~ ~/.~~ - Jo5 o.J\ ~. ~ J:: 0.6 
- --,Q .!I 
:g 
g ol\ \ \ ~--:- IIJ:: :: \-\. ' ~ -~~::: 
~ ........_ 0 .2 0 .21-\ \. ------ -10.2 
~9 9 
-=:IQ . I 0.1 I- \ "-. --1 0. I 
9.2 
0.4 
>~?----.7 
?·3 
0 .1 p.1 
0 0 . 0 0 . Q 
0 .1 0 .2 o .3 0.4 o .5 o .6 o .7 o .e 0.9 1.0 o 0. 1 0 .2 o .3 o .4 o .5 o .e o.7 o .e o .9 1:0 o 0 .1 0 .2 0 . 3 o .4 o .s o .6 0.1 o .s o.9 1.a 
t =c;?, mon\l,s t = 6 months t = 12 monttis 
\l;, g~-'2 ~ ~~-~m~te ~w_a~t !~~µt·f( fiiA~}hVJ~J~ i@.Ub?'tdiitfiD\iti); ,l_iikti;i.JttJ,) ( =.~ 1,'.i t~ ~J ~ JM1-~~ ,~;{JP.QPt;ijS:I·_AA~iii;::~ {-1:2. mp1;ithS: ~f- -r~dial _ @rgc-sthlill <~~blls,;,Jidaiion .. 
8The extremely high variations of the void ratio 
-&~)coefficient of permeability in the upper area of 
, ,,,._·.d osit must be accurately captured; the pro-
. ffie ep d · h ' · h 
',t~, ed multiple no e per time-step tee mque 1s muc 
P-~ ffici·ent for filling analysis than the classical one 
:moree . 
:Ed(ie teehruque. 
n.'(J) Since the physic~) and mat~ematical as~ects 
( the problem are highly non-Jmear, non-lmear 
R./ . h ds rather than the classical linearization 
:met o 
;''" hod should be employed; the latter does not 
met . 1 . l . ;'i-· re convergence, may imp y uneconom1ca t1me-
assu 1 . h "d . 
r.;• 5 and in genera , overestimates t e vo1 ratio step ' · I d h ( h' thus underestimatmg set~ ement an strengt t 1s 
·::-,· ·mportant for reclamation purposes). 
15 t4) A second order non-linea~ method (Newto°:-
Raph~on) is necessary to ~scertam ~onvergence, this 
:'"tovides a robust and efficient algonthm; however , a 
~ood initial approximation of the solution and a first 
·[~ration with the 'fixed-point iterative method' are 
~so recommended. 
'.(~) The bi-dimensional extension ampl ifies the 
non~Iinear aspects, however the 'split-step method ' 
allows a ready generalization of the one-dimens ional 
techniques, and successful results are also obtained 
~ith'·the developed non-linear algorithm. 
REFERENCES 
J. R. E. Gibson, G . L. England and M . H. L. Husse y, The 
· theory of one-dimensional consolidation of saturated 
·clays, I. Finite nonlinear consolidation of thin homo-
. geneous layers. Geotechnique 17, 261-273 (1967). 
'2:, J. L. Monte and R . J . Krizek , One-dimensional · 
. · ~thematical model for large-strain consolidation . 
Geotechnique 26, 495-510 (1976). 
3. F . Somogyi , Analysis and prediction of phosphatic 
clay consolidation: implementation package . Technical 
report ; Florida Phosphatic ··cla y Research Project , 
Lakeland, FL (1979) . 
4. R. L. Schifffman, Finite iind infinitesimal strain consol-
idation . J. Geotech. Engng Div., ASCE 106, 203-20 7 
(1980). 
5. R . E. Gibs on, R . L. Schiffman and K. W. Cargill , The 
theory of one-dimensi onal consolidation of satura ted 
clays, II . Finite nonlinear consolidation of thick homo-
geneous layers. Can. Geotech. J. 18, 280- 293 (1981). 
6. R. J. Krizek and F . Somogyi, Perspectives on modelling 
consolidation of dredged materials . Proceedings of the 
Symposium on Sedimenta tion/ Consolidatiori Models , 
ASCE, San Francisco , CA , pp . 296-332 (1984). 
7. R. L. Schiffman, V. Pane and V. Sunara, Sedimentation 
and consolidation. Proceedings of the Symposium 
on Sedimentation /Consolidation Models , ASCE , San 
Francisco , CA, pp. 57-121 (1984). 
8. J. R . Feldkamp, Numerical analysis ofone-dimensional 
nonlinear large-strain consolidation by the finite ele-
ment method . J. Porous Media 4, 239-257 (1989). 
9. S. D . Koppula and Morgens tern, On the consolidation 
of sedimenting clays . Can. Geotech. J . 19, 260-268 
(1982). 
10. W. D . Ca rrier, L. G. Bromwell and F. Somogyi , Design 
capacity of slurried mineral waste ponds. J. Geotech. 
Engng Div., ASCE 109, 699- 716 (1983). 
11. A. Huerta , G . A. Kri egsmann and R . J. Krizek, Per-
meabil ity and compressibility of slurries from seepage-
induced consolidation. J. Geotech. Engng, ASCE 114, 
614-{,27 (1988). 
12. Symposium on Consolidation and Disposal of 
Phosphat ic and Other Was te Clays . Department of 
Civil Engineering , University of Flo_rida , Gainsville , 
Florida . Sponsored by the Florida Institute of Phos-
phate Research , Bartow , FL (1987). 
13. F . Somogyi, W. D. Carrier, J. E. Lawyer and J. F . 
Beckman, Waste phosphatic clay disposal in mine 
cuts . Proceedings of the Symposium on Sedimentation / 
Consolidat ion Models, ASCE , San Francisco , CA 
(1984).