Exact solution of Terzaghi's consolidation equation and extension to two/three-dimensional cases (3th versione)

Terzaghi's one-dimensional consolidation equation simulates the visco-elastic behaviour of soils depending on the loads applied as it happens, for example, when foundation are laid and start carrying the weight of the structure. Its application is traditionally based on Taylor's solution that approximates experimental results by introducing non-theoretical variables that, however, contradict the actual behaviou of soils. After careful examination of the theoretical and experimental aspects connected with consolidation, the proposal of this research is a solution consisting in a non-linear equation that can be considered correct as it meets both mathematical and experimental requirements. The solution proposed is extended to include differential equations relating to two/three dimensional consolidation by adopting a transversally isotropic model more consistent with the inner structure of soils. Finally, this essay is complete with application examples that give more reliable results than the traditional solution. Future developments are also highlighted considering that the uniqueness theorem has not been proven yet.


INTRODUCTION TO THE HYDRAULIC AND MECHANICAL BEHAVIOUR OF SOILS
The mechanical behaviour of soils was coded only after the introduction of the "concept of effective stresses" [1] which marked the birth of Soil Mechanics starting from the general structure of Continuum Mechanics from which it derives. This concept is based on the inner structure of soils that are composed of a solid skeleton and inter-particle gaps. These pores are more or less interconnected and through them run fluids of different nature. Therefore, in view of a necessary simplification of the mathematics of associated phenomena, the concept of effective stresses requires the soils to be assimilated to bi-phasic systems composed of a solid skeleton saturated with water, i.e. two continuous means that act in parallel and share the stress status [4]: In equation (1) there is the tensor of the total stresses exerted by the solid skeleton (σij), the hydrostatic pressure exerted by the fluid (u0, known as interstitial pressure) and Kronecker's delta (δij); furthermore, from a merely phenomenological point of view, equation (1) attributes the soil shear resistance only to effective stress, independent of the presence of the fluid. It should be highlighted that the concept of effective stresses is valid only in stable conditions, when the fluid is in balance with the solid skeleton. In these conditions you can calculate the hydrostatic component in all points of the underground and in all moments through the application of the laws of balance; vice-versa, in transient conditions, it is necessary to introduce other elements capable of accounting for the variation of the component u induced by stresses of various nature. At this point, the problem focus is on the permeability coefficient (K = m/sec) that, by expressing the capacity of a soil to transmit a fluid, takes on the character of a velocity and varies approximately in the range 10 -1 ÷10 -10 m/sec depending on the inner structure of the solid skeleton. As a direct consequence of this extreme variability, the soils with high permeability (such as sands) behave as open hydraulic systems where compression induced, for example, by the load of a foundation, causes simultaneous drainage of the fluid from the pores. In practice, the fluid does not take part in the mechanical response and the stress induced weighs only on the solid skeleton that, in turn, subsides in association with reduced porosity. On the contrary, soils with very low permeability (such as clays) exhibit hydraulic delay in reacting to stresses, with consequent development of an initial interstitial overpressure (u ≠ 0) that contradicts equation (1), and participate in the mechanical response with the solid skeleton. A transient filtering motion follows that comes to end only when the initial value of the interstitial pressure is reset (u = 0). In practice, given that a deformation of the solid skeleton occurs together with the expulsion of water, sands develop elasto-plastic settlements synchronous with load application while clays exhibit time-dependent consolidation settlements typically characterized as reverse hyperbolic functions. Figure 1. The load transmitted by a foundation always causes interstitial overpressures whose dissipation depends on soil permeability that in turn depends on porosity. Consolidation may last from some minutes (loose sands) to tens of years (very stiff clays).

INTRODUCTION TO THE TERZAGHI'S THEORY OF ONE-DIMENSIONAL CONSOLIDATION
The one-dimensional consolidation equation [1,2] describes the hydraulic behaviour of soils in transient conditions by making it possible to simulate the variation in time of interstitial overpressures (u), generated -for example -by the load induced by a foundation or by a road embankment (figure 1), with consequent viscoelastic settlements to which corresponds a structural reorganisation of the solid skeleton, with reduction of porosity and, concurrently, of the degrees of freedom. This formulation can be inferred from applying the continuity equation to supposedly saturated soils leading [5], with a few mathematical manipulations, to the following relationship that demonstrates how the transient filtering motion depends on the vertical permeability factor (Kz), on the compressibility factor (mv) and on the weight of the volume of water (γw = 10 kN/m 3 ) that in turn identifies the fluid: The next step consists in introducing the hypothesis (in contrast with experimental results) that the permeability factor does not change during consolidation: Finally, denoting the consolidation coefficient as cv: you come to write the classical one-dimensional consolidation equation: It should be noticed that equation (5) is analogous to Fourier's law on heat propagation to the point that you can define the theory of consolidation as the simulation of the propagation of stress-induced interstitial pressures in the subsoil.

TAYLOR'S SOLUTION TO THE ONE-DIMENSIONAL CONSOLIDATION THEORY
The solution of equation (5) may be found introducing two non-dimensional variables: Both expressed as a function of H that is the maximum drainage path of the fluid. In this manner, equation (5) becomes: whose analytical solution is [3]: valid for m ∈ Ν e based on the assumption that u remains constant with depth. According to this solution, the application of a load causes an interstitial overpressure in the soil (initial conditions: t = 0, u(0) = u0 with 0 ≤ Z ≤ 2), that rapidly drops to zero at the drainage surfaces (contour conditions: t>0, u(Z=0) = u(Z=2) = 0), generating an isochrone (regarding the concept of effective stresses in a transient phase) whose bottom coincides with the central line of the consolidating layer that, in turn, can be considered as a waterproof surface (figure 2a). On the basis of the above reasoning, the solution can be extended also to the cases with drainage at one edge, provided that attention is paid at identifying the correct value of H (figure 2b). Since drainage resulting from consolidation leads to yielding, it is also possible to define the average consolidation degree: expressed as a function of yielding in an instant (s (t)) relative to final yielded (sf). Furthermore, also the following relationship is demonstrated to be valid: whose solution is shown in figure 3 as a function of Tv. Alternatively, you can use the approximate solutions by Sivaram and Swamee [7]: To conclude, Taylor's solution corresponds to a development in series of Fourier truncated at the top, which describes consolidation through a reverse hyperbolic function coded through fixed values of Tv expressed as a function of Um.

EXPERIMENTAL DEFINITION OF THE THEORY OF CONSOLIDATION
The application of equation (5) is linked to the execution of oedometric tests (figure 4), whose first prototype dates back to Terzaghi [8], Casagrande [9], Gilboy [10] and Rutledge [11]. Only vertical stresses and deformations (σ'z, εz) can be assessed in these tests because side deformations (εr) are actually impossible and radial stresses (σ'r) cannot be measured. In this manner, its constitutive equation can be written and expressed in differential terms, in the following simplified form: where appears the compliance matrix [C], valid for oedometric conditions. See details in [12,13]. Equation (13) can also be written in the extended form: Therefore, the dσ'r expression can be inferred from equation (14b), as follows: that in turn can be introduced into equation (14a): If you analyse equations (15) and (16), you will find out that the former provides the formal relationship of the vertical and radial stresses while the latter demonstrates that vertical deformations depend on the combination of axial and radial visco-elastic constants.
In brief, as far as the phenomenological meaning of the above formulations is concerned, oedometric tests can be fruitfully used to analyse the one-dimensional effects induced by an unspecific stress assuming that the conditions of geometrical symmetry illustrated in figures 1 and 3 are valid, and taking into account that equation (16) can be also written as follows: In this manner, you work out the formulation of the elastic module in oedometric conditions (Eed) that, in turn, is in inverse proportion to the compressibility factor (mv) shown in equation (4). Finally, it should be noticed that the oedometric condition is a special case of flat deformation, in the absence of side deformations and corresponds to a null value of Poisson's factor (ν = 0), associated only to the existence of geometrical symmetries of the type shown in figure 1.

EXPERIMENTAL DETERMINATION OF THE CONSOLIDATION COEFFICIENT
The application of the theory of consolidation is strictly connected to the determination of cv, coded according to two methods, both connected with the oedometric tests [14]. These tests consist in subjecting a round test piece [15], confined into an indeformable ring (figure 4), to axial loads according to a loading schedule (Fa) following a geometric progression q like: F a n ( ) = q ! F a n"1 ( ) Each load is held for a time that guarantees the passage from load-dependent primary consolidation to secondary consolidation that depends on the merely rheological behaviour of soils.
The most reputed methodology [16] requires the identification of the following elements on the experimental curves (figure 5a): 1) t100: time necessary for conventional completion of primary consolidation; 2) t50: time to attain the condition Um = 50%, equation (10), in relation with instrumental double drainage of the test-piece (figure 4); 3) H100: height of the test-piece at the end of conventional primary consolidation, calculated starting from initial height H0 = 20 millimeters from which the value ΔH100 read in the diagram is subtracted. In this manner, the second equation (6) can be reversed: as a function of both Tv(50%) = 0.197 (figure 3) and semi-height H100/2 of the test piece in each loading step applied (figure 5b).
Known cv for the design load (meaning that it should be selected depending on the stress level expected), the next step is the calculation of the time of completion of consolidation: being Tv(95%) = 1.129 ( figure 4) and H the height of the layer according to that highlighted in figure 2.
Finally, known the consolidation settlements (that can be calculated using Geotechnics methods based on the interpretation of the oedometric tests), you can build the time-settlement curve for the different structure/foundation nodes from which to obtain the differential settlement variation field (figure 6).

LIMITS OF TAYLOR'S SOLUTION
If you analyse Taylor's solution [3], you will find out that its application relies on using the non-dimensional variable Tv, whose values are coded independent of the soil analysed, and that is used both in experimental determination of the consolidation coefficient (cv) and in the calculation of the consolidation completion time. At this point, a problem arises from experimental observation, namely the mechanical behaviour of a soil depending on its geological history (that is, a tenso-deformative history), it cannot be referred to the fixed values of a nondimensional variable. On the other hand, it suffices to reflect on clay formation environments, which depend on hydrolysis which in turn is influenced by temperature and water, i.e. by latitude, altitude and the climate system, to understand the intrinsic limits of using the parameter Tv. Also the construction of the time-settlement curve is quite difficult as it requires the following passages: -calculation of time (ti) for each value of Um, by introducing the corresponding values of Tv from figure 3 into equation (20); -product of the final settlement (sf) and the single percentage values of Um from which to obtain the corresponding values of settlement si; -assembly of the diagram si-ti ( figure 6). In the end, being not possible to simulate the real mechanical and hydraulic behaviour of soils correctly, we consider the Taylor's solution as experimentally inadequate. In other words, the exact solution of equation (5) would depend exclusively on the consolidation coefficient cv that, expressed as a function of Kz and of mv, is a real constitutive link. To conclude, it should also be considered that, once determined cv for each load step applied, the experimental practice foresees the calculation of the corresponding Kz, taken from equation (4), being known the value of mv which is obtained by applying equation (17) to the experimental curves (figure 5a). However, using this procedure includes that, due to the remarkable non-linearity of the compressibility factor, Kz markedly depends on the stress level and, consequently, on the deformative level. Furthermore, since the method is connected with the consolidation velocity through a direct relationship with cv, rather inaccurate measures are obtained due to the physical magnitudes that play a role in the process as well as to their non-linearity.

EXACT SOLUTION OF TERZAGHI'S CONSOLIDATION EQUATION
Let's assume that u, cv, kz > 0 are three positive constants assigned and that: Now, you can easily notice that u(z,t) solves the differential equation (5); indeed, given the validity of the following: you obtain: which is the differential equation given. If you analyse equation (22), you will find out that it simulates the time variation of interstitial overpressures in the subsoil through the consolidation constant kz, -starting from the point where they are triggered -by dampening their width as depth increases through a reverse hyperbolic function. Now, since equation (22) is the solution of equation (5), it should be necessarily extended also to experimental methods -oedometrically considered -in order to determine the parameters that govern it correctly. In this sense, it may be useful to analyse an important property of the function u.
Given H > 0 (figure 4), it can be useful to identify the instant tH > 0 to which the following conditions apply: which has infinite solutions like: or like: The next passage consists in selecting the positive value t closest to zero from among those determined by setting the condition t > 0 that provides: from which we obtain: The last passage includes that, if you set also the following condition: you obtain the formulation of the consolidation completion time: Finally, from equation (34) you can extract the consolidation coefficient: Now it should be noticed that for equation (35) to allow the determination of cv starting from experimental data, the additional parameter kz introduced in the initial hypotheses must be defined correctly. To this regard, considering that equation (5) after setting: as a function of the angular frequency ω = (2π)/t. The parameter sought can be taken from equation (37): In accordance with the work hypotheses, it depends exclusively on the consolidation coefficient as it should not be determined or added as an external condition to differential equation (5). Now equation (38) can be introduced into equation (35) and with a few mathematical manipulations you obtain: What is left to do now is just to assess correspondence of equation (39) and the experimental data that can be obtained from oedometric tests by setting the following conditions as a function of figure 4 and of the diagram in figure 5a: These conditions, once introduced into equation (39), give the parameter sought without adding any further hypotheses or additional variables: As an example, with H0 = 0.02 meters, t50 = 912 seconds and H100 = 0.01668 meters, obtained from a real-life test, you obtain:  equation (19) → cv = 1.5×10 -8 m 2 /s  equation (41) → cv = 3.7×10 -8 m 2 /s To conclude, the solution proposed with equation (22) can be considered correct even though the uniqueness theorem has not been proven yet, since it perfectly corresponds both mathematically and experimentally. Obviously, by analogy, the same solution can be extended also to Fourier's equation on heat transmission, after replacing the consolidation coefficient with the thermal diffusivity coefficient (cv ≡ α) and temperature (u ≡ T) to interstitial overpressures.

EXTENSION TO THE TWO-DIMENSIONAL AND THREE-DIMENSIONAL CASES
Let's put u, cx, cz, kx, kz > 0 to be positive constants assigned and be: Going through the same passages as in equations (23) to (26), you can easily notice that u(x,z,t) solves the differential equation: Similarly, if u, cx, cy, cz, kx, ky, kz > 0 are positive constants assigned and: (45) the regular function given, then: u x, z, t ( ) = u ! e " k x x"k y y"k z z cos 2 c x k x 2 + c y k y 2 + c z k z solves the differential equation: Considering the soils' markedly anisotropic hydraulic and mechanical behaviour that can be attributed to the structure of the solid skeleton that in turn depends on accumulation characteristics (figure 7), the application of the solution to the 2D case -equation (43) -does not necessarily requires the preliminary determination of all parameters introduced. Mathematically speaking, the geological development of sedimentary beds leads to assimilating soils to transversally isotropic means to which the following conditions apply: Using equation (49) requires the experimental determination of permeability in the horizontal plane (Kx), such as in the three-axial cells or permeameters, by applying the rotation of the axes to the test-pieces. It follows that the relationship in equation (4) can be extended also to the consolidation coefficients: and to consolidation constants: Finally, by repeating the same procedure as in equations (36) and (37), you obtain: where to introduce equations (50) and (51) that lead initially to: and then, by introducing equation (53) into equation (51), to the determination of the parameter sought: Notice that with the condition m = 1, equation (54) reduces to: The procedure outlined can be fruitfully extended even to the 3D case by writing in sequence: Furthermore, since the horizontal plane symmetry properties of transversal isotropic means Kx = Ky and kx = ky apply, equations (60), (61) and (62) can be simplified: To conclude, with the condition m = n = 1, the equations above reduce to: i.e. to:

PRACTICAL APPLICATIONS
At the time being, the applications of the solutions provided by equations (22), (43) and (46) led to the development of two methods to analyse consolidation development over the time.  Then, the calculation is achieved by building the curve Um÷t, noticing that each dissipation curve has a matching depth zi that voids the corresponding u (figure 9b), as it happens with the consolidation completion curve tz=10 = 17.5 years. Therefore, considering the following should apply: finally you can determine Um values for each dissipation curve, which leads to the diagram in figure 9a. However, if the application of the one-dimensional solution is immediate, with 2D consolidation it is necessary to identify the relationships expressed by equations (48) to (51) beforehand. Consequently, in order to evaluate how those parameters influence interstitial overpressure dissipation times, analyses were conducted presuming:  m = 1, 2, 5, 10  n = 0,1H, 0,5H, H, 2H, as the horizontal dimension was not known a priori. The results (figure 10a) demonstrate that as m increases, consolidation times decrease for each value of n taken, in line with the behaviour expected for a two-dimensional means. At the same time, as n increases, you observe that consolidation completion times increase and this can be attributed to a longer drainage path. If you work out a diagram with the results in the m÷log(t) plane (figure 10b), you perfectly realise the influence exerted by both m and n, considering the solution of the one-dimensional problem as a upper limit. This condition, matched with the experimental value of m, allows defining the range of action of the most probable value of n: see the grey area in figure 10 for a hypothetical m = 3. In other words, for a horizontal layer infinitely extended in that direction, the value of n to be introduced for the correct solution, depend on foundation size and consequent stress level induced in the subsoil.

CONCLUSIONS
With in mind the entire Mechanics of Soils, and the study of the soil visco-elastic behaviour [1,2] in particular, the application of Terzaghi's differential equation, is historically based on Taylor's solution [3] that approximates experimental results -limited to the one-dimensional case only -through the introduction of arbitrary and fixed non-dimensional variables, independent of the geological history of the means. After accurate examination of the theory and of the historical solution, and interpretation of lab data, this research work makes a proposal for a solution that can be considered correct as it solves the differential equation and, at the same time, allow correct interpretation of experimental data. Then, solution has been fruitfully extended to the two-and three-dimensional cases and finally tested in a real-life case. To conclude, results even more satisfactory have come to light from the analysis of data. At the same time, they have opened additional research channels, considering that the uniqueness theorem has not been proved yet and that the influence of underground geometry should evaluated in the two-and three-dimensional cases. Monitoring of future geotechnical structures and/or retrospective analysis of real-life cases may clarify these residual aspects at a later time.

CREDITS
The author wishes to thank Luca Lussari, Department of Mathematics and Physics of Università Sacro Cuore di Brescia -Italy, for his suggestions and observations.