Powerful Software to Simulate Soil Consolidation Problems with Prefabricated Vertical Drains

: The present work describes the program Simulation of Consolidation with Vertical Drains (SICOMED_2018), a tool for the solution of consolidation processes in heterogeneous soils, with totally or partially penetrating prefabricated vertical drains (PVD) and considering both the effects of the smear zone, generated when introducing the drain into the ground, and the limitation in the discharge capacity of the drain. In order to provide a completely free program, the code Next-Generation Simulation Program with Integrated Circuit Emphasis (Ngspice) has been used as a numerical tool while the Matrix Laboratory (MATLAB) code was used to program and create an interface with the user through interactive screens. In this way, SICOMED_2018 is presented as an easy-to-use and intuitive program, with a simple graphical interface that allows the user to enter all the soil properties and geometry of the problem without having to resort to a complex software package that requires programming. Illustrative applications describe both the versatility of the program and the reliability of its numerical solutions.


Introduction
The use of specific software for the understanding of complex physical phenomena involved in engineering problems has become an indispensable tool in modern teaching, mainly due to the high costs involved in laboratory tests. These programs, in spite of the effort and knowledge-especially of programming and calculation-that their development require, are generally designed by the university departments themselves as an interdisciplinary work among several areas of knowledge. The basic objective is to create an interactive computer-to-user communication program that, with minimal help, enables the learning objectives to be met in a shorter time than would be necessary with laboratory tests. In this sense, the design of these programs-avoiding the concept of the 'black box'-allows the user to access the choice of a large number of parameters associated with programming and design of the model, with the purpose of determining the potential influence of them in numerical results.
The 'network simulation' research group of the UPCT (Technical University of Cartagena) has significant experience in software development, in the fields of mechanical engineering and coupled flow and transport processes, such as: (i) Program for Designing Simple Fins (PRODASIM) and Program for Heat Conduction (PROCCA-09), for the assembling of transient heat transfer in fins [1]; (ii) Flow and Transport Simulator (FATSIM-A), for the simulation of flow and transport processes [2]; and (iii) Coupled Ordinary Differential Equations by Network Simulation (CODENS_13), for the simulation of mechanical models governed by coupled ordinary differential equations [3]. All of them are available and have been used as tools for numerous research works [4,5], but require a license to use the Personal Computer Simulation Program with Integrated Circuit Emphasis (Pspice), which is responsible for performing the numerical calculations.
In this work we present the free program SICOMED_2018 [6], a code for the simulation of transient problems of soil consolidation with the use of prefabricated vertical drains to accelerate the settlement process, easy to use but with numerous applications within the field of civil engineering. The possibility of selecting not only the value of any of the physical parameters involved, with heterogeneous soils of up to three layers, but also those associated with the grid size and the numerical computation parameters allows the user to address specific consolidation scenarios-determining the separate influence of each of the parameters involved-as well as to solve optimization problems such as the relationship between the consolidation time and the depth of the drain. On the other hand, the powerful graphic output environment of the program allows us to visualize and interpret results in a direct and simple way. For these reasons, SICOMED_2018 is presented as a complete software package, with novel contributions, that joins other specific programs for learning the theory and practice of soil consolidation such as Consolidation Testing and Multimedia (CTM) [7].
The performance of SICOMED_2018, which can be acquired through UPCT (Technical University of Cartagena), requires the installation of the free code Ngspice [8]. The graphical communication environment has been developed with MATLAB (under a license that allows its distribution), a very widespread code at a scientific and professional level that has highly developed tools to create user interfaces, while the powerful computational algorithms of Ngspice provide the nearly exact solution of the network model of the consolidation process, whose design is based on the analogy between electric and physical quantities. The design of the numerical model follows the rules of the network simulation method [9], a tool widely known in the scientific world for the implementation of network models. Each term in the equation is an electric current that is balanced with the other terms in a common node of the volume element whose voltage is the solution of the potential quantity (pressure) of the governing equation. A mesh size of the order of 50 volume elements, in one-dimensional (1D) problems, reduces errors with respect to analytical solutions to below 1% [10], which is acceptable for this type of engineering problem. Illustrative applications are presented to demonstrate the efficiency and reliability of the program.

Fundamentals of Soil Consolidation with Prefabricated Vertical Drains and Network Model
Prefabricated vertical drains (PVD) are a system of soil drainage commonly used in civil engineering in problems of soil consolidation in which the goal is to accelerate final settlement. When a fully saturated clay soil is subjected to the action of loads on its surface, it begins to experience a decrease in volume or settlement over a generally long period, in a transitory process commonly known as consolidation. This volume reduction is due to the gradual expulsion of part of the water contained in the voids between the soil particles, through one or more draining boundaries (generally across the ground surface), as the excess pore pressure initially generated by the applied load is transmitted to the soil, increasing the effective pressure on it.
In real scenarios, with ground of considerable thickness and low permeability, the main problem facing the civil engineer is the long consolidation times needed to reach final settlement (on the order of years). For this reason, the installation of PVD is a good solution to reduce the time of consolidation of a soil, as we are significantly shortening the drainage length of the water.
Assuming the hypotheses postulated by Terzaghi [11]-(i) both water and soil particles are assumed incompressible; (ii) the ground self-weight is neglected; (iii) the excess pore pressure (u) is caused by the application of the external load (q), constant over time; (iv) the fluid movement in the porous medium is assumed to obey the law of Darcy; (v) the soil skeleton does not creep under the action of a constant effective pressure; (vi) the change in soil thickness in relation to its total volume (hypothesis of 1 + e constant) is considered negligible; (vii) there is uniform applied load, so that the soil deformation only occurs in the vertical direction; (viii) the geotechnical parameters of the soil, included in the consolidation coefficient c v , remain constant; and (ix) the soil rests on an impermeable stratum-so that water drains only by the ground surface (water flow only in the vertical direction), the governing equation and boundary conditions for the 1D linear soil consolidation problem are given by The analytical solution for a single layer of soil can be found in several works [12,13], providing solutions for the local degree of settlement (U) for any depth and instant of time, as well as the average degree of settlement (U) at any time.
The installation of PVD, made of geosynthetic material, allowing drainage through them, involves the transformation of the 1D problem into a three-dimensional (3D) problem, since the direction or flow of water occurs in any of the three spatial directions (x, y, z). By adopting a rectangular layout (Figure 1), which allows us to obtain volume elements of rectangular geometry, the mathematical model for the typical section ( Figure 2) is governed by the following equations: u (x,y,z=0,t) = u (0≤x≤c,y=0,0≤z≤d,t) = 0 Soil surface and drainage area of the PVD (6) ∂u ∂n rest of the boundaries = 0 Impermeable edges (7) u (x,y,z,t=0) = u o Initial condition, where ∂u ∂n represents the partial derivative in a direction normal to the boundary surface.
Water 2018, 10, x FOR PEER REVIEW 3 of 20 stratum-so that water drains only by the ground surface (water flow only in the vertical direction), the governing equation and boundary conditions for the 1D linear soil consolidation problem are given by The analytical solution for a single layer of soil can be found in several works [12,13], providing solutions for the local degree of settlement (U) for any depth and instant of time, as well as the average degree of settlement (U ) at any time.
The installation of PVD, made of geosynthetic material, allowing drainage through them, involves the transformation of the 1D problem into a three-dimensional (3D) problem, since the direction or flow of water occurs in any of the three spatial directions (x, y, z). By adopting a rectangular layout (Figure 1), which allows us to obtain volume elements of rectangular geometry, the mathematical model for the typical section ( Figure 2) is governed by the following equations: Governing equation u ( , , , ) = u Initial condition, where represents the partial derivative in a direction normal to the boundary surface. Assuming isotropy in the (x, y) horizontal directions (a common feature of soil behavior), Equation (5) can be simplified to = c , + c , + . Assuming isotropy in the (x, y) horizontal directions (a common feature of soil behavior), Equation (5) can be simplified to At this point, it is important to note that, as in the linear theory of Terzaghi [11], the consolidation coefficients c v,z and c v,h of the governing Equation (9) have been assumed constant. As is well known, this is not entirely true, since in practice the soil properties (k v , k h , e and a v ) that are grouped by these coefficients ( are not constant, but they vary during the consolidation process, a phenomenon that began to be analyzed many years ago by Davis and Raymond [14] and many other authors of the time. At present there are numerous publications and research that incorporate the non-linearity of these coefficients in the mathematical models, although it is also true that, for the usual practical values of the surface applied load (q), the changes in properties that give rise the consolidation coefficients (c v,z and c v,h ) are small (and the variations in these properties may even compensate for each other), so that the final result is practically constant consolidation coefficients. In effect, the non-linearity of these coefficients grows the higher the σ f /σ o quotient, a ratio that, though it theoretically can reach any value, in practice very rarely exceeds the value of 2 (or, at most, 2.5); a value below which the difference between the solutions provided by the linear and non-linear models is practically insignificant [15]. In this way, the smearing effects produced by the introduction of the vertical drain are reflected in the mathematical model in a very simple way, since, in practice, the soil property that is affected by the smear, for the purposes of the consolidation coefficient, is the hydraulic conductivity [16][17][18]. So, and for the area of soil that is considered to be affected by the smear effects, the hydraulic conductivity values of the smear zone (k v,s and k h,s ) must simply be assigned. At this point, it is important to note that, as in the linear theory of Terzaghi [11], the consolidation coefficients cv,z and cv,h of the governing Equation (9) have been assumed constant. As is well known, this is not entirely true, since in practice the soil properties (kv, kh, e and av) that are grouped by these coefficients (c , = ( ) , c , = ( ) ) are not constant, but they vary during the consolidation process, a phenomenon that began to be analyzed many years ago by Davis and Raymond [14] and many other authors of the time. At present there are numerous publications and research that incorporate the non-linearity of these coefficients in the mathematical models, although it is also true that, for the usual practical values of the surface applied load (q), the changes in properties that give rise the consolidation coefficients (cv,z and cv,h) are small (and the variations in these properties may even compensate for each other), so that the final result is practically constant consolidation coefficients. In effect, the non-linearity of these coefficients grows the higher the σ′f/σ′o quotient, a ratio that, though it theoretically can reach any value, in practice very rarely exceeds the value of 2 (or, at most, 2.5); a value below which the difference between the solutions provided by the linear and non-linear models is practically insignificant [15]. In this way, the smearing effects produced by the introduction of the vertical drain are reflected in the mathematical model in a very simple way, since, in practice, the soil property that is affected by the smear, for the purposes of the consolidation coefficient, is the hydraulic conductivity [16][17][18]. So, and for the area of soil that is considered to be affected by the smear effects, the hydraulic conductivity values of the smear zone (kv,s and kh,s) must simply be assigned. The set of Equations (6)-(9) has a complex analytical solution (from the point of view of the management for the civil engineer), so it is presumed much more appropriate to solve them by numerical methods.

Network Model
The network method is a technique for the study and numerical solution of any physical processes that can be defined by a mathematical model. Starting from the latter, the procedure consists of two stages: (i) elaborate a "network model" or electrical circuit equivalent to the process, and (ii) numerically simulate the model by means of an adequate electrical circuits resolution code. Its applications are multiple; González-Fernández [9] uses the network method in solving transport problems through membranes, heat transfer and electrochemical systems. In recent years it has been successfully applied in other engineering fields such as elastic waves [19], material resistance [20], corrosion [21], magneto-hydrodynamics [22], flow and transport [4], tribology [23], mechanical problems and chaos [24]. On the other hand, among the future research lines in which the network method can be used we find: seepage [25], real-life case models in hydrologic engineering [26][27][28][29], The set of Equations (6)-(9) has a complex analytical solution (from the point of view of the management for the civil engineer), so it is presumed much more appropriate to solve them by numerical methods.

Network Model
The network method is a technique for the study and numerical solution of any physical processes that can be defined by a mathematical model. Starting from the latter, the procedure consists of two stages: (i) elaborate a "network model" or electrical circuit equivalent to the process, and (ii) numerically simulate the model by means of an adequate electrical circuits resolution code. Its applications are multiple; González-Fernández [9] uses the network method in solving transport problems through membranes, heat transfer and electrochemical systems. In recent years it has been successfully applied in other engineering fields such as elastic waves [19], material resistance [20], corrosion [21], magneto-hydrodynamics [22], flow and transport [4], tribology [23], mechanical problems and chaos [24]. On the other hand, among the future research lines in which the network method can be used we find: seepage [25], real-life case models in hydrologic engineering [26][27][28][29], and atmospheric dispersion of pollutants [30,31], some of them already in the development phase.
The formal equivalence between the network model and the physical process is that both are governed by the same equations (discretized in space, but retaining time as a continuous variable), establishing a correspondence between the dependent variables of the problem, pressure and water flow in the physical model, and the electrical variables of the circuit, voltages, and currents. That is, both models are governed by the same equations referred to a volume (or cell) element and by the same discrete equations for the boundary conditions. Consequently, the errors from the simulation are only attributable to the geometric mesh size, so that for an acceptable number of cells (of the order of 50 in 1D scenarios, 4000 for 3D problems) errors are much less than 1% in linear problems [10]. The powerful computation codes for the numerical resolution of circuits provide the exact solution of these thanks to an optimal selection of the calculation time steps, imposed by the Ngspice code or by the user.
In terms of the first spatial derivative, Equation (9) is written in the form where z + and z − denote the cell output and input locations in the z-direction, and the same for x + and x − in the x-direction and for y + and y − in the y-direction. Dividing the medium into N x × N y × N z volume elements of size ∆x × ∆y × ∆z, and using the nomenclature of Figure 3, the above equation can be written in the form Water 2018, 10, x FOR PEER REVIEW 5 of 20 establishing a correspondence between the dependent variables of the problem, pressure and water flow in the physical model, and the electrical variables of the circuit, voltages, and currents. That is, both models are governed by the same equations referred to a volume (or cell) element and by the same discrete equations for the boundary conditions. Consequently, the errors from the simulation are only attributable to the geometric mesh size, so that for an acceptable number of cells (of the order of 50 in 1D scenarios, 4000 for 3D problems) errors are much less than 1% in linear problems [10]. The powerful computation codes for the numerical resolution of circuits provide the exact solution of these thanks to an optimal selection of the calculation time steps, imposed by the Ngspice code or by the user.
In terms of the first spatial derivative, Equation (9) is written in the form where z and z denote the cell output and input locations in the z-direction, and the same for x and x in the x-direction and for y and y in the y-direction. Dividing the medium into Nx × Ny × Nz volume elements of size Δx × Δy × Δz, and using the nomenclature of Figure 3, the above equation can be written in the form Each of the terms in the above equation is a flow: Each of the terms in the above equation is a flow: so that Equation (12) can be written as a flow balance: In the electric analogy, by establishing associations between variables so that the excess pore pressure (u) is equivalent to the electrical voltage (V) while the temporal (first derivative) and spatial changes (second derivative) of u are two electric currents (J), we have: which is equivalent to Ohm's law, which relates the current (I R ) in a resistance with the voltage (V R ) at its ends (terminals): Thus, the implementation of Equation (15) is carried out by the resistance between the central node and the upper edge of the cell (Figure 3). The value of this resistance is By the same reasoning, the term is implemented by a new resistor of the same value as Equation (17) between the central node and the lower edge of the cell. Finally, the term is equivalent to the equation that relates the electric current at a capacitor (I C ) to the voltage or potential difference at its ends (V C ): Thus, the implementation of this term of Equation (12) is carried out by a capacitor C of value unity connected between the central node and an outer node (which will be the common reference node for all cells). In short, the values of the electrical devices that implements the model of the volume element, whose network is depicted in Figure 4, are Notice that, physically, an increase in the consolidation coefficient decreases the value of the resistances and makes the consolidation process faster, as expected.
Each cell of the domain is connected to the adjacent cells by means of ideal electrical contacts. As regards the boundaries, a second-type Neumann (homogeneous) condition is adopted at the impermeable edges, which is simply implemented by a resistance of infinite value, while a first-type Dirichlet condition (constant pressure) is assumed at the draining boundaries (soil surface and drainage area of the PVD), which is implemented by a constant voltage source of zero value.
Finally, the drain discharge capacity (q w ) limitation caused by the unit gradient of pressure at the drain is modeled by fixing the drain permeability (hydraulic conductivity) to the required value. Note that this solution allows us to increase the flow of water through the drain as the water moves towards the exit border of the drain.

The Simulation Program
The objective of this program (SICOMED_2018) is the simulation of the 3D consolidation problem with partially or totally penetrating PVD, in heterogeneous soils formed by one, two, or three layers and taking into account the effects of the smear zone and the discharge capacity of the drain. It must be an easy-to-use software for the user but able to provide powerful numerical calculations. This requirement entails the creation of a practical and simple graphical interface, where data entry and choice of the different options is performed in an orderly, intuitive, and guided way. Also, the program should include a complete guide [32] that allows the user to understand the successive steps in data entry, simulation and representation of results.
With regard to the simulation possibilities, the user can use SICOMED_2018 both to design and to optimize the installation of PVD in the ground, since the program allows us to know, in a simple way, the great variety of results that will allow engineers to make the decisions that most suit their needs: evolution of settlements on the soil surface, excess pore pressure at any point of the ground, and average degree of settlement, among others. Figure 5 shows a simplified scheme of the SICOMED_2018 flowchart.

The Input Data and Network Design
Data entry is carried out through the following three windows that appear successively as they are completed (Figures 6-8). The start-up screen (Figure 6), allows the entry of the problem geometry (values a, b, and c in Figure 2), the definition of the smear zone, and the desired grid size in plant. In the second window (Figure 7), the geometry of the strata (thickness and vertical grid size), parameters to define the consolidation coefficients (both the undisturbed and the smear zone), the depth of the PVD, and its discharge capacity are introduced. Finally, in the third screen (Figure 8), parameters related to the simulation such as the value of the uniform applied load, initial and final calculation times, the maximum time step between each iteration, and a relative tolerance parameter used by Ngspice to achieve the convergence are introduced. As a unit of measure of time, the year (y) has been chosen, a more appropriate magnitude than the second(s) for this type of problem, while for the other magnitudes the units of the international system have been maintained.

The Simulation Program
The objective of this program (SICOMED_2018) is the simulation of the 3D consolidation problem with partially or totally penetrating PVD, in heterogeneous soils formed by one, two, or three layers and taking into account the effects of the smear zone and the discharge capacity of the drain. It must be an easy-to-use software for the user but able to provide powerful numerical calculations. This requirement entails the creation of a practical and simple graphical interface, where data entry and choice of the different options is performed in an orderly, intuitive, and guided way. Also, the program should include a complete guide [32] that allows the user to understand the successive steps in data entry, simulation and representation of results.
With regard to the simulation possibilities, the user can use SICOMED_2018 both to design and to optimize the installation of PVD in the ground, since the program allows us to know, in a simple way, the great variety of results that will allow engineers to make the decisions that most suit their needs: evolution of settlements on the soil surface, excess pore pressure at any point of the ground, and average degree of settlement, among others.

The Input Data and Network Design
Data entry is carried out through the following three windows that appear successively as they are completed (Figures 6-8). The start-up screen (Figure 6), allows the entry of the problem geometry (values a, b, and c in Figure 2), the definition of the smear zone, and the desired grid size in plant. In the second window (Figure 7), the geometry of the strata (thickness and vertical grid size), parameters to define the consolidation coefficients (both the undisturbed and the smear zone), the depth of the PVD, and its discharge capacity are introduced. Finally, in the third screen (Figure 8), parameters related to the simulation such as the value of the uniform applied load, initial and final calculation times, the maximum time step between each iteration, and a relative tolerance parameter used by Ngspice to achieve the convergence are introduced. As a unit of measure of time, the year (y) has been chosen, a more appropriate magnitude than the second(s) for this type of problem, while for the other magnitudes the units of the international system have been maintained. Figure 5 shows a simplified scheme of the SICOMED_2018 flowchart.

Simulation and Output Data
To increase the flexibility and effectiveness of the program, the user can avoid entering data by working with another file previously saved. For this, SICOMED_2018 includes the options of saving and loading data. Once all the data have been entered or loaded the simulation can start. To do this, SICOMED_2018 creates the network model file in a specific source code that is run in Ngspice. At the end of the simulation the user can access the results given in a graphic way. SICOMED_2018 offers up to six possibilities of results representation (Figure 9), conveniently arranged so that the geotechnical engineer can have all the necessary information in a simple, precise form:

1.
Excess pore pressure in a given column of the soil 2.
Excess pore pressure in a given point of the soil 3.
Average degree of settlement 4.
Local settlements in a given column of the soil 5.
Total settlement in a given point of the surface 6.
Surface settlements animation.
Water 2018, 10, x FOR PEER REVIEW 10 of 20

Simulation and Output Data
To increase the flexibility and effectiveness of the program, the user can avoid entering data by working with another file previously saved. For this, SICOMED_2018 includes the options of saving and loading data. Once all the data have been entered or loaded the simulation can start. To do this, SICOMED_2018 creates the network model file in a specific source code that is run in Ngspice. At the end of the simulation the user can access the results given in a graphic way. SICOMED_2018 offers up to six possibilities of results representation (Figure 9), conveniently arranged so that the geotechnical engineer can have all the necessary information in a simple, precise form: 1. Excess pore pressure in a given column of the soil 2. Excess pore pressure in a given point of the soil 3. Average degree of settlement 4. Local settlements in a given column of the soil 5. Total settlement in a given point of the surface 6. Surface settlements animation. Some of these representations are shown in the illustrative applications of the following section, which have been selected with the aim of presenting all the potentialities and utilities of SICOMED_2018, such as: (i) analysis of classic consolidation problems; (ii) real scenarios of PVD layout optimization; (iii) inclusion of the smearing and drain discharge capacity effects; and (iv) verification of results.

First Scenario: Consolidation of a One-Layer Soil with PVD
This application refers to a fixed one-layered consolidation scenario in which the PVD layout is already imposed. The geometry of the problem and the soil properties are shown in Table 1 (note that the PVD thickness, td, is not given, as this parameter does not have an influence when the smear and drain discharge capacity effects are not considered, as in this case). The final time of the simulation has been set at four years, while the chosen time step is 0.01 years. Finally, a 20 × 16 × 20 grid has been chosen, enough to guarantee precision in the results, with relative errors below 1%.  Some of these representations are shown in the illustrative applications of the following section, which have been selected with the aim of presenting all the potentialities and utilities of SICOMED_2018, such as: (i) analysis of classic consolidation problems; (ii) real scenarios of PVD layout optimization; (iii) inclusion of the smearing and drain discharge capacity effects; and (iv) verification of results.

First Scenario: Consolidation of a One-Layer Soil with PVD
This application refers to a fixed one-layered consolidation scenario in which the PVD layout is already imposed. The geometry of the problem and the soil properties are shown in Table 1 (note that the PVD thickness, t d , is not given, as this parameter does not have an influence when the smear and drain discharge capacity effects are not considered, as in this case). The final time of the simulation has been set at four years, while the chosen time step is 0.01 years. Finally, a 20 × 16 × 20 grid has been chosen, enough to guarantee precision in the results, with relative errors below 1%. At this point it is important to note that SICOMED_2018 has the option to directly enter the values of the consolidation coefficients (c v,z and c v,h ), allowing only the analysis of the problem in terms of the duration of the process (options 1, 2, and 3 in Figure 9), or provide the values of the soil parameters (a v , e o , k v and k h ) from which these coefficients are deduced-Equation (22)-allowing this second option to perform an analysis of the problem in terms of settlement (options 4, 5, and 6 in Figure 9): Once the simulation is run, with a computation time of 10 min on an Intel Core i7-6700 CPU 4.00 GHz computer, the user can access the results (Figure 9). The excess pore pressure can be represented for all the cells of one or two soil columns, given to the program by their (X, Y) coordinates ( Figure 10). As can be seen, the dissipation of the excess pore pressure ( Figure 11) occurs faster in those points of the soil closer to the draining edges, soil surface, and PVD. The program also allows (option 2) the representation of the evolution of this unknown in one, two, or three locations, indicated to the SICOMED_2018 by its (X, Y, Z) coordinates. At this point it is important to note that SICOMED_2018 has the option to directly enter the values of the consolidation coefficients (cv,z and cv,h), allowing only the analysis of the problem in terms of the duration of the process (options 1, 2, and 3 in Figure 9), or provide the values of the soil parameters (av, eo, kv and kh) from which these coefficients are deduced-Equation (22)-allowing this second option to perform an analysis of the problem in terms of settlement (options 4, 5, and 6 in Figure 9): Once the simulation is run, with a computation time of 10 min on an Intel Core i7-6700 CPU 4.00 GHz computer, the user can access the results (Figure 9). The excess pore pressure can be represented for all the cells of one or two soil columns, given to the program by their (X, Y) coordinates ( Figure  10). As can be seen, the dissipation of the excess pore pressure ( Figure 11) occurs faster in those points of the soil closer to the draining edges, soil surface, and PVD. The program also allows (option 2) the representation of the evolution of this unknown in one, two, or three locations, indicated to the SICOMED_2018 by its (X, Y, Z) coordinates.    At this point it is important to note that SICOMED_2018 has the option to directly enter the values of the consolidation coefficients (cv,z and cv,h), allowing only the analysis of the problem in terms of the duration of the process (options 1, 2, and 3 in Figure 9), or provide the values of the soil parameters (av, eo, kv and kh) from which these coefficients are deduced-Equation (22)-allowing this second option to perform an analysis of the problem in terms of settlement (options 4, 5, and 6 in Figure 9): Once the simulation is run, with a computation time of 10 min on an Intel Core i7-6700 CPU 4.00 GHz computer, the user can access the results (Figure 9). The excess pore pressure can be represented for all the cells of one or two soil columns, given to the program by their (X, Y) coordinates ( Figure  10). As can be seen, the dissipation of the excess pore pressure (Figure 11) occurs faster in those points of the soil closer to the draining edges, soil surface, and PVD. The program also allows (option 2) the representation of the evolution of this unknown in one, two, or three locations, indicated to the SICOMED_2018 by its (X, Y, Z) coordinates.   The evolution of the average degree of settlement (the degree of consolidation averaged over all soil columns) is shown in Figure 12. Thus, it can be observed that 90% of the consolidation is reached after 1.30 years.
Water 2018, 10, x FOR PEER REVIEW 12 of 20 The evolution of the average degree of settlement (the degree of consolidation averaged over all soil columns) is shown in Figure 12. Thus, it can be observed that 90% of the consolidation is reached after 1.30 years. Options 4-6 of Figure 9 are also available (Figures 13-15) since parameters av, eo, kv, and kh are given. Note that the soil columns closer to the PVD will reach the final settlement before (Figure 14), although the value of this unknown is the same for all the columns: 0.0927 m.   Options 4-6 of Figure 9 are also available (Figures 13-15) since parameters a v , e o , k v , and k h are given. Note that the soil columns closer to the PVD will reach the final settlement before (Figure 14), although the value of this unknown is the same for all the columns: 0.0927 m. The evolution of the average degree of settlement (the degree of consolidation averaged over all soil columns) is shown in Figure 12. Thus, it can be observed that 90% of the consolidation is reached after 1.30 years. Options 4-6 of Figure 9 are also available (Figures 13-15) since parameters av, eo, kv, and kh are given. Note that the soil columns closer to the PVD will reach the final settlement before (Figure 14), although the value of this unknown is the same for all the columns: 0.0927 m.   The evolution of the average degree of settlement (the degree of consolidation averaged over all soil columns) is shown in Figure 12. Thus, it can be observed that 90% of the consolidation is reached after 1.30 years. Options 4-6 of Figure 9 are also available (Figures 13-15) since parameters av, eo, kv, and kh are given. Note that the soil columns closer to the PVD will reach the final settlement before (Figure 14), although the value of this unknown is the same for all the columns: 0.0927 m.   The last representation option ( Figure 15) is a video animation of the settlement evolution of the whole soil surface. As representation options, the user can choose the time interval to display and the number of frames per second of the animation (for faster or slower display).

Second Scenario: Optimization of PVD Layout
In this application we analyze and accelerate the process of consolidation of a three-layered soil (whose properties are given in Table 2) on which an embankment is going to be built, transmitting to the ground a uniform load of 60 kN/m 2 . Using SICOMED_2018, different layouts (depth of penetration and separation of PVDs) will be simulated (with different final and step times, depending on the case), looking for the optimal solution from the point of view of consolidation time. A 10 × 10 × 30 grid has been chosen (0.2 m high in each cell). To address this problem it is necessary to set a definition for the consolidation time (or characteristic time, to [33]), as well as to fix a target value of this time to be achieved with the use of PVDs. So, we define the characteristic time as the time that takes the soil to reach an average degree of settlement (U ) of 90%. In addition, PVDs must reduce this time to a value below two years. In this way, we proceed to perform the simulation of the problem in the most unfavorable scenario, that is, when there are no PVDs installed in the ground. For this case (Figure 16), 90% of the average degree of settlement is achieved after more than 30 years, a value too high for work in civil engineering. The last representation option ( Figure 15) is a video animation of the settlement evolution of the whole soil surface. As representation options, the user can choose the time interval to display and the number of frames per second of the animation (for faster or slower display).

Second Scenario: Optimization of PVD Layout
In this application we analyze and accelerate the process of consolidation of a three-layered soil (whose properties are given in Table 2) on which an embankment is going to be built, transmitting to the ground a uniform load of 60 kN/m 2 . Using SICOMED_2018, different layouts (depth of penetration and separation of PVDs) will be simulated (with different final and step times, depending on the case), looking for the optimal solution from the point of view of consolidation time. A 10 × 10 × 30 grid has been chosen (0.2 m high in each cell). To address this problem it is necessary to set a definition for the consolidation time (or characteristic time, t o [33]), as well as to fix a target value of this time to be achieved with the use of PVDs. So, we define the characteristic time as the time that takes the soil to reach an average degree of settlement (U) of 90%. In addition, PVDs must reduce this time to a value below two years. In this way, we proceed to perform the simulation of the problem in the most unfavorable scenario, that is, when there are no PVDs installed in the ground. For this case (Figure 16), 90% of the average degree of settlement is achieved after more than 30 years, a value too high for work in civil engineering.  Figure 17 shows the local settlement evolution of an entire soil column, which is different for each stratum. The simulation also provides a final surface settlement value of 0.53 m. In order to reduce the consolidation time, PVDs separated every 2 m (both at x and y spatial directions) are installed. The width of the PVD is 0.1 m, a typical commercial value for this type of drain of geosynthetic origin. As depth of penetration we begin with d = 1 m (partially penetrating drain to the first layer). As the reduction of to is not enough, we increase the penetration up to 4 m (partially penetrating drain to the intermediate stratum) and 6 m (totally penetrating drain). The results for the consolidation time (to) are shown in Table 3. Table 3. Values of to as a function of PVD penetration (a = 1 m; b = 1 m; c = 0.05 m).    In order to reduce the consolidation time, PVDs separated every 2 m (both at x and y spatial directions) are installed. The width of the PVD is 0.1 m, a typical commercial value for this type of drain of geosynthetic origin. As depth of penetration we begin with d = 1 m (partially penetrating drain to the first layer). As the reduction of to is not enough, we increase the penetration up to 4 m (partially penetrating drain to the intermediate stratum) and 6 m (totally penetrating drain). The results for the consolidation time (to) are shown in Table 3. In order to reduce the consolidation time, PVDs separated every 2 m (both at x and y spatial directions) are installed. The width of the PVD is 0.1 m, a typical commercial value for this type of drain of geosynthetic origin. As depth of penetration we begin with d = 1 m (partially penetrating drain to the first layer). As the reduction of t o is not enough, we increase the penetration up to 4 m (partially penetrating drain to the intermediate stratum) and 6 m (totally penetrating drain). The results for the consolidation time (t o ) are shown in Table 3.  Since none of the studied options fits the two-year objective, it is necessary to bring PVDs closer. Thus, for the cases d = 6 m and d = 4 m, a set of simulations in which both the distance between PVDs and the separation between rows of PVDs are reduced will be simulated with SICOMED_2018. Characteristic times for each layout are summarized in Table 4. Since none of the studied options fits the two-year objective, it is necessary to bring PVDs closer. Thus, for the cases d = 6 m and d = 4 m, a set of simulations in which both the distance between PVDs and the separation between rows of PVDs are reduced will be simulated with SICOMED_2018. Characteristic times for each layout are summarized in Table 4. Table 4. Values of to as a function of PVD penetration and surface layout (c = 0.05 m).

Third Scenario: Influence of the Smear Zone and the Discharge Capacity Limitation of the PVD
The third application addresses the use of limited discharge capacity drains, as well as the influence of the smear zone. The geometric characteristics and ground and drain properties are shown in Table 5. The geometric values of the problem (a, b, c, and td) have been selected so that the solutions of SICOMED_2018 can be compared with the results obtained by Hansbo et al. [16] and Barron [34], in the same problem but using a two-dimensional (2D) radial geometry, with an equivalent diameter of the PVD and a diameter of influence of 0.062 m and 0.95 m, respectively. The physical soil and drain properties are identical. Regarding av and e, any pair of values that provides the same consolidation coefficients is valid, having adopted the values of av = 1.224 × 10 −5 m 2 /N and e = 1. Finally, ∆Q = 10.000 N/m 2 . Four consolidation scenarios have been simulated, comparing the average degree of consolidation achieved in three typical times: 0.5, one, and two years (Table 6). In the first case, in which the smear effects are not considered and the drain is supposed to have an infinite discharge capacity, the results reproduce those obtained analytically by Barron [34] under the free strain hypothesis, with relative deviations below 3% for the three times. The other three cases assume the

Third Scenario: Influence of the Smear Zone and the Discharge Capacity Limitation of the PVD
The third application addresses the use of limited discharge capacity drains, as well as the influence of the smear zone. The geometric characteristics and ground and drain properties are shown in Table 5. The geometric values of the problem (a, b, c, and t d ) have been selected so that the solutions of SICOMED_2018 can be compared with the results obtained by Hansbo et al. [16] and Barron [34], in the same problem but using a two-dimensional (2D) radial geometry, with an equivalent diameter of the PVD and a diameter of influence of 0.062 m and 0.95 m, respectively. The physical soil and drain properties are identical. Regarding a v and e, any pair of values that provides the same consolidation coefficients is valid, having adopted the values of a v = 1.224 × 10 −5 m 2 /N and e = 1. Finally, ∆Q = 10.000 N/m 2 . Four consolidation scenarios have been simulated, comparing the average degree of consolidation achieved in three typical times: 0.5, one, and two years (Table 6). In the first case, in which the smear effects are not considered and the drain is supposed to have an infinite discharge capacity, the results reproduce those obtained analytically by Barron [34] under the free strain hypothesis, with relative deviations below 3% for the three times. The other three cases assume the same smear length but consider different PVD discharge capacities. For them, the comparison is established with the results obtained by Hansbo et al. [16], based on their empirical expressions for the average degree of consolidation. Deviations, particularly for relatively high times, can be considered acceptable, taking into account that these are not numerical or analytical results. These deviations grow slightly to 11.6% for small times and low discharge capacities, perhaps because the models degrees of freedom (3D for SICOMED_2018 and 2D radial for Hansbo) are not the same. On the other hand, the results are coherent with respect to the correlation between the average degree of consolidation and the influence of both the smear effects and the discharge capacity limitation of the drain.

Final Comments and Conclusions
SICOMED_2018 has been created as free software for solving 3D consolidation problems with partial or totally penetrating vertical drains, in heterogeneous soils and considering the influence of the smear effects and the discharge capacity limitation of the drain. All the unknowns of interest for the civil engineer, such as the excess pore pressure evolution, the average degree of consolidation, local, and total settlements and even a video animation of the soil surface deformation, can be obtained with SICOMED_2018 in a relatively simple way. Therefore, it can be considered a program with professional, educational, and research applications.
Following network simulation method rules, SICOMED_2018 designs network models (electrical circuits) whose governing equations (with spatial discretization and retaining time as a continuous variable) are formally equivalent to the discretized equations of the consolidation problem. These network models are solved by an electrical circuit resolution code such as Ngspice, directly providing the solutions to the consolidation problem.
The program has been developed in a pleasant communication environment, both for the data entry and the graphical presentation of results, under a classic windows environment programmed with MATLAB, a code whose license is not free but allows the distribution of the programs created with it. On the other hand, for the numerical calculation the powerful computational algorithms integrated in Ngspice, an electrical circuits resolution free code that provides practically exact solutions to the models designed by SICOMED_2018, have been used. The combination of these codes has allowed us to obtain completely free software.
To illustrate the operation of the program, three applications, corresponding to real cases, are presented. The first one deals with the simulation and resolution of a consolidation scenario with the PVD distribution already imposed (a classic problem), showing the evolution of the main variables of interest (excess pore pressure, average degree of consolidation, and local and total settlements). The second one concerns an optimization problem in a soil formed by three strata with different properties, in which it is sought to find the optimum PVD distribution in order to reduce the characteristic time of consolidation to a certain target value. Finally, the third application serves as a verification of the program against the analytical and empirical solutions provided by other authors, studying the influence that the smear zone and the discharge capacity limitation of the drain have on the consolidation problem. Deviations are reduced in all cases, showing the versatility and reliability of the program.
Finally, from the work presented here, a series of challenges that could be addressed in future research are suggested. On the one hand, in the field of soil consolidation, this research could be completed (or improved) by implementing additional options in the program, such as analysis based on real tabulated data, or even the possibility of addressing inverse problems. On the other hand, the wide experience of the 'network simulation' research group of the UPCT opens up the possibility of using a similar methodology in the creation of free software in disciplines as diverse as seepage, atmospheric dispersion of pollutants, and hydrologic engineering, among others. U average degree of settlement (dimensionless) u i excess pore pressure in node i (N/m 2 ) u o initial excess pore pressure (N/m 2 ) V C voltage between the terminals of a capacitor (V) V R voltage between the terminals of a resistance (V) x long horizontal spatial coordinate (m) X value of the spatial coordinate x (m) y wide horizontal spatial coordinate (m) Y value of the spatial coordinate y (m) z vertical spatial coordinate (m) Z value of the spatial coordinate z (m) γ w specific weight of water (N/m 3 ) σ f final effective pressure (N/m 2 ) σ o initial effective pressure (N/m 2 )