Numerical modeling of settlements for shallow foundations on layered soils

: A numerical and an empirical method for calculating nonlinear load-settlement curves for shallow foundations on sand are examined and used in a new methodology. Both methods have merits and drawbacks. The drawbacks are overcome by the methodology proposed and verified in the paper. This methodology combines the merits of each method. For this purpose, a modification of the empirical method is proposed that accounts for the finite initial soil stiffness at very small strains. Computer programs with sophisticated nonlinear stress-strain relationships, such as Hardening Soil Small in Plaxis 2D, which are versatile in solving complex foundation problems, can cover strains from very small to large. When the foundation soil is layered, it is proposed to fit such a numerical load-settlement curve against the modified empirical relationship for each sand layer separately. This requires cone penetration tests, measurements of the shear wave velocity, and basic laboratory tests. The methodology is described and applied at two locations where load tests on footings were carried out. At one location there was a top layer of clay, which was also taken into account. The results of the application of the proposed methodology are very good.


INTRODUCTION
Load-settlement curves of shallow foundations, spanning the range from very small loads to those that induce soil failure, are nonlinear. Existing methods for the calculation of this nonlinear relationship, both numerical and empirical, are reviewed and commented on. One numerical (Benz 2007) and one empirical (Mayne et al. 2012) method are combined in the paper in order to match the two load-settlement curves by adjusting the preloading parameter in the numerical analysis for each sand layer separately. Then, a complex numerical soil-structure interaction analysis with layered foundation soil including clay layers can be performed. The Mayne et al. (2012) method is modified in the paper in order to include the granular soil behavior at very small strains. A high soil stiffness at very small strains has long been considered as essential for obtaining reliable predictions of soil deformations in geotechnical engineering (e.g. Burland 1989). For this purpose, it is required to know the soil shear modulus at very small strains, which can easily be determined from measurements of the shear wave velocity in the foundation soil. The only other required parameter is the cone resistance from a cone penetration test (CPT).
In the present study, the commercial finite element computer program Plaxis 2D (Brinkgreve et al. 2011) was used for numerical modeling with the constitutive (stress-strain) relationship Hardening Soil Small (HSSmall) developed by Benz (2007). The HSSmall model allows for modeling the proper reduction of stiffness when straining the soil from very small strains up to large strains at failure. The required soil parameter is the reference shear modulus at very small strains (  10 -6 ), which equals the maximum shear modulus at atmospheric confining pressure. Other sand parameters can be determined from correlations with the sand relative density (Brinkgreve et al. 2010), and the relative density in turn is obtained from the correlation with the CPT cone resistance (Jamiolkowski et al. 1985). Parameters for clay can be obtained from correlations with the CPT cone resistance (Mayne and Rix 1993;Sanglerat 1972).
In the second step, the preloading parameter is adjusted for each sand layer by fitting the numerical load-settlement curve for the model consisting only of the material in that layer against the modified Mayne et al. (2012) curve. A footing of size close to that from the database used by Mayne et al. (2012), say, 2 m, is used.
The proposed methodology was verified at two locations with footing load tests on layered foundation soils. The first location was Green Cove Springs, Florida (Anderson et al. 2007), and the second location was Alvin (East), Texas (Tand et al. 1994). At the second location, there was a top clay layer. In the absence of shear wave velocity measurements at either location, the sand stiffness at very small strains was determined. The application of the proposed methodology gave very good results when comparing measured settlements with results of numerical analyses for layered foundation soils.

Numerical methods
Generally, nonlinear load-settlement curves may be predicted in several steps. These are as follows: obtaining a sufficient number of high-quality undisturbed samples representative of the foundation soil, stress-strain testing of these samples in a laboratory, interpreting the test results to obtain parameters for an appropriate constitutive soil model, and using these parameters to run a finite element computer program that incorporates this constitutive Lukić Kristić, I., Prskalo, M., Szavits-Nossan, V.

Numerical modeling of settlements for shallow foundations on layered soils
model. This procedure can be used to solve underlying nonlinear boundary value problems for successively increasing foundation loads. However, it heavily relies on representative undisturbed soil samples that are difficult to obtain and mount into a laboratory apparatus when the foundation soil is a particulate coarse grained material such as sand.
An alternative approach is to correlate parameters of a constitutive model with in situ tests. In that respect, Byrne et al. (1987) used the hyperbolic constitutive model of Kondner (1963) and Duncan and Chang (1970). Byrne et al. (1987) correlated parameters of this model with the relative density of sand, and field penetration test results [standard penetration test (SPT) and CPT] obtained in this sand. Correlations were based on an examination of laboratory and field measurements available in the literature. The examined laboratory triaxial tests were conducted on disturbed sand samples at various relative densities. Then, the researchers carried out numerical finite-element back analyses of field observations of foundation settlements on sand in order to check whether the laboratory test results were representative of the field conditions. They concluded that the stiffness parameters obtained from triaxial tests are appropriate for settlement prediction by their hyperbolic law at low sand relative density, whereas the sand stiffness is considerably underpredicted at high relative densities. Therefore, this procedure proved not to be suitable for predicting foundation settlements at higher relative densities.
A similar but much more sophisticated approach was reported by Benz (2007) 20 years later. He developed an elastic plastic hardening constitutive model called HSSmall that accounts for the remarkable reduction of the high initial stiffness of sand by about 70 % when it is strained from very small strains (in the range of 0.0001 %) up to small strains of about 0.1 %. This high stiffness and its large reduction with straining became apparent in laboratory tests only when appropriate techniques and strain gauges were developed that explained previous systematic overpredictions of field deformations based on the standard laboratory testing of stiff and medium-stiff soils (Burland 1989). The HSSmall model can also account for soil hardening induced by preloading (loading and unloading), thus simulating an overconsolidated soil.
In order to validate his constitutive model, which was incorporated into the Plaxis 2D finite element computer program, among other structures, Benz (2007) analyzed the loadsettlement curve of one of the experimental square footings on sand tested at the location of Texas A&M University (Briaud and Gibbens 1994;. Five load tests were performed on square footings ranging in size from 1 m to 3 m, embedded 0.76 m into an 11-m-thick layer of medium dense, fairly uniform silty fine silica sand, underlain by stiff clay. Load tests were carried out until the soil settlement reached 15 cm for each footing. Extensive tests at the experimental site, including measurements of the shear wave velocity and penetration resistance by cone penetration tests, as well as extensive laboratory tests on samples recovered from the site, including triaxial and resonant column tests, were performed (Briaud and Gibbens 1997). The foundation settlements were measured throughout the load tests. Benz (2007) determined numerous soil parameters for the HSSmall model from the results of triaxial and resonant column tests on samples reconstructed to the field sand density. The reference shear modulus at very small strains was determined both by a resonant column test (laboratory parameters) and as a considerably larger value for in situ conditions (in situ parameters). Using the Plaxis 2D program to model the square 3 m × 3 m footing North as an equivalent circular footing, Benz calculated load-settlement curves for both values of the shear modulus at very small strains. He assumed that the sand at the experimental site was normally consolidated (NC), i.e., that it never underwent a load-unload cycle in its geological history. Fig. 1 compares both newly calculated curves with Benz's parameters and the measured load-settlement data. Although the calculated curve using in situ conditions is closer to the measured data, it still considerably overpredicts the measured settlements. Soil is assumed to be normally consolidated (NC). Measured settlements from Briaud and Gibbens (1997) In order to match the measured data, Benz (2007) assumed that the sand at the site had been temporarily preloaded, i.e., loaded and then unloaded by the same uniform surface pressure, thus rendering the sand overconsolidated (OC). Preloading the sand was obviously only an assumption, since there was no means to prove it independently. Benz was ambiguous by not stating the value of preloading. Therefore, for the purpose of the present research, load-settlement curves were recalculated with the Plaxis 2D program using Benz's parameters. For the newly calculated curves, a preloading of 160 kPa was used in order to make the foundation soil overconsolidated. It is stressed, however, that the value of preloading the soil in the recalculated numerical analyses was reached by trial and error until a good match between the newly calculated and Benz's load-settlement curves was obtained.
The newly calculated load-settlement curves by Plaxis 2D, using the same HSSmall model parameters as Benz's (in situ shear modulus at very small strains), and preloading by 160 kPa, are shown in Fig. 2 of the next subsection for all five experimental A&M footings.
It is seen in this subsection that examples of Byrne et al. (1987) and of Benz (2007) demonstrate that the predictions of load-settlement curves of footings on sand based on laboratory tests of sand samples are unreliable.

Empirical methods
A completely different approach to predict load-settlement curves of foundations on sand was proposed by Akbas and Kulhawy (2009a;2009b). They used measured load-settlement curves of 167 prototype foundations at 37 locations with homogeneous sand as the foundation soil. The prototype foundations were square or rectangular with a width-to-length ratio close to 1. Akbas and Kulhawy (2009a) derived a hyperbolic type of load-settlement relationship for footings on sand having one parameter that is related to the bearing capacity of the foundation sand given by Vesić (1975). The bearing capacity is related to the friction angle of sand, which in turn has to be determined from a suitable correlation with results of some in situ tests. The advantage of this method is in avoiding the need for the laboratory testing of sand samples.
Following the line of avoiding the need for laboratory testing of sand, Mayne et al. (2012) developed a simple expression relating the nonlinear load-settlement curve of a square footing on sand and the average cone resistance of the standard CPT performed at the footing location. Mayne et al. (2012) analyzed full-scale load tests on 31 square or nearly square footings with widths of 0.5 m ≤ B ≤ 6 m situated on 13 different quartz-silica sand sites, and derived a simple one-parameter power law function for load-settlement curves. They correlated the parameter of this function directly to the average CPT cone resistance at all footing sites, and obtained the following relationship between the average applied footing pressure p, the footing settlement s, square footing size B, and average cone resistance q c : with the coefficient of determination r 2 = 0.933. They recommended that the value of q c used in equation (1) should be determined as the average cone resistance of the soil below the footing and down a distance of 2B. Fig. 2 compares load-settlement curves calculated by Eq.
(1) with measured data for the five experimental footings of the A&M test site, which has an average cone resistance of q c = 7.5 MPa, as taken by Mayne et al. (2012). The comparison indicates the level of accuracy that can be expected when using equation (1) in predicting load-settlement curves for footings on sand, and it is very good. Fig. 2 also shows Plaxis calculations after the foundation soil was preloaded at 160 kPa (OC), with the same parameters as those that Benz (2007) used for the 3 m × 3 m footing North (in situ value of the maximum shear modulus). Despite the obvious advantage of using a direct correlation between the load-settlement curve and results of a common field test in avoiding laboratory tests, this approach is not suitable for cases with layered soils or flexible spread foundations that require the numerical modeling of soil-structure interaction. Furthermore, equation (1) implies an infinite soil stiffness starting at a zero load. This is a drawback if it is used in numerical back analyses to match the analytical curve because sand exhibits finite stiffness at very small strains, along with large reductions in stiffness in the range of small strains, as in the case of the HSSmall model discussed in the preceding subsection. To overcome this drawback, a modified loadsettlement function that represents an extension of Mayne et al. (2012) is proposed in the next section.  (1) with average cone resistance of q c = 7.5 MPa, and newly calculated by Plaxis 2D with in situ value of initial sand stiffness and preloading of 160 kPa (overconsolidated)

Development of load-settlement function
Sand at very small strains behaves as an elastic material. Therefore, it is proposed herein to extend the Mayne et al. (2012) equation (1) to account for high foundation soil stiffness at very small strains. The elastic settlement s e of a rigid square footing with size B loaded with a uniform pressure p, embedded by a depth D f into a homogeneous elastic layer of thickness h underneath the footing base, the maximum Young's modulus E 0 constant with depth, and Poisson's ratio  can be calculated from the following expression suggested by Mayne and Poulos (1999): where: I E and I G are correction factors for the foundation embedment and the layer thickness, respectively, and are given by   With an increasing load, plastic strains begin to develop. For additional settlements owing to plastic strains, a power law similar to that of equation (1) is added to equation (2): where p L is the value of pressure p at s/B = 0.1, i.e., it is the average contact pressure referred to as the bearing capacity (e.g., Briaud and Gibbens 1997).
Parameter a can be determined by introducing p = p L into equation (6), giving Parameters b and p L are determined by a linear regression analysis of pressurenormalized settlement pairs measured during each of 16 load tests performed at four locations of sand, as reported in the literature. Then, parameter b is determined as the average value for the 16 footings, and a correlation is set between p L for all footings and the CPT cone resistance at the four locations, as shown in the next subsection.

Database of 16 load tests at four locations, and statistical analysis
The correlation between p L and q c is the same as that of Mayne et al. (2012). Thus, it is expected that the proposed load-settlement curve of equation (8) gives very similar results to those obtained by equation (1). This was checked for the four locations used herein, and was proven correct. However, differences between the two methods emerge in the range of very small and small strains. where G 0 is the maximum shear modulus,  is the soil density, and v s is the average shear wave velocity. Then, E 0 = 2G 0 (1+).
In Fig. 4, there is a correct slope of the tangent to the curve at zero load calculated by the proposed method. This slope corresponds to the sand stiffness at very small strains, E 0 . On the other hand, the tangent to the curve calculated by Mayne et al. (2012) is horizontal, indicating an infinite soil stiffness at very small strains. The addition of the elastic response of the foundation sand for very small loads, as given by Eq. (8), is preferable when the curve is used as a benchmark in back analyses to be performed by numerical modeling that incorporates the behavior of sand at very small and small strains, as described in the next section.

PROPOSED METHODOLOGY FOR CALCULATING LOAD-SETTLEMENT CURVES
For complex shallow foundations on layered soil, none of the previously published methods of calculation for the load-settlement curve can be reliably applied, as explained in previous sections. The Mayne et al. (2012) method works very well for rigid square footings of sizes in the range from 0.5 m to 6 m on a relatively homogeneous sand layer. The only parameter required for this method is the average CPT cone resistance of the foundation soil. If a loadsettlement curve has to be predicted for a similar footing or similar ground conditions, neither the proposed modified load-settlement curve nor sophisticated numerical modeling would enhance the predictions. However, if the ground contains sand layers of markedly different stiffness or even contains layers of fine-grained soil, or the depth below the foundation to a very stiff stratum is smaller than double the foundation width, or the spread foundation is e-Zbornik 19/2020.

Numerical modeling of settlements for shallow foundations on layered soils
flexible (such as foundation slabs and plates that necessitate a numerical soil-structure interaction analysis), direct empirical methods for predicting load-settlement curves such as those described in the present paper are of little use. Similarly, finite-element or likewise numerical methods with complex constitutive models for soils, despite being powerful modeling tools, have drawbacks in the case of coarse soils owing to problems in recovering undisturbed soil samples and their handling in the laboratory. Furthermore, it was shown that an a priori unknown preloading of the foundation soil is necessary in order to match the measured data. This paper advocates a third approach that combines the merits of numerical modeling with the modified Mayne et al. (2012) method given by equation (8). No stress-strain laboratory tests for coarse-grained soil layers are required. In the case of fine-grained soil layers, undisturbed samples should be retrieved for laboratory stress-strain testing. In the absence of laboratory tests, correlations with the CPT cone resistance can be used. The only two in situ tests that should be carried out are the CPT and geophysical measurements of the shear wave velocity. The proposed approach for a reliable numerical soil-structure interaction analysis requires the following steps: 1. Perform field CPT and shear-wave-velocity measurements at the future foundation site. Identify layers in the soil profile. 2. Use the values of the shear wave velocity profile to determine the shear modulus at very small strains for the numerical modeling of each soil layer, and for use in equation (8) for sand. This equation is suggested if the constitutive model used in the numerical analysis has the capability of capturing the nonlinear stress-strain behavior in the range of very small to large strains (such as HSSmall of the Plaxis program). Determine other required parameters for numerical modeling by using correlations with the CPT cone resistance. 3. For each coarse-grained layer in the soil profile, numerically calculate the loadsettlement curve for an imaginary rigid square, or equivalent circular footing having a similar size to those described herein, that rests on the entire foundation soil being a homogeneous layer of the same soil. Adjust the preloading input parameter in order to match the load-settlement curve calculated from equation (8) for each coarsegrained soil layer. 4. Use the model parameters determined in steps 2 and 3 for each soil layer, and perform a numerical analysis on the complex soil-structure interaction problem or settlement analysis for the layered soil profile.

VERIFICATION
Two locations with simple footing configurations and layered foundation soils, where load tests were carried out, were chosen to verify the proposed methodology. The first location is Green For a coarse-grained soil with c = 0, G 0 and G 0 ref are related through (9) where σ 3  is the minor principal stress and p ref is the reference pressure, usually taken as the atmospheric pressure (100 kPa). A similar relationship holds for all reference values above. When shear wave velocities v s are measured throughout the foundation soil profile, it is possible to determine the G 0 profile, and then G 0 ref from equation (9) for each distinct layer. This is the starting point for numerical back analyses. In the absence of v s measurements at the two analyzed locations, correlations between G 0 and q c were used first. Schnaid et al. (2004) (after Schnaid 2009) gave the following correlation for sand: where σ v0  is the in situ vertical effective stress relevant to the CPT cone resistance in a layer, and p a is atmospheric pressure equal to 100 kPa. The coefficient  varies from 110 to 800, the value of 110 being the lower bound for uncemented sand, and the value of 800 the upper bound for cemented sand. In the absence of any relevant indication which value of  to use, in this study a value of  = 400 was adopted for all sand layers at the Green Cove Springs location, and  = 250 for all sand layers at the Alvin (East) site.
Values of q c provided for each sand layer by Anderson et al. (2007) and Tand et al. (1994) were taken as q c1 . In situ vertical effective stresses were calculated at the midpoint of each layer in the layered foundation soil. It was then assumed that for a thick homogeneous foundation soil consisting of the sand of one layer, the CPT cone resistance varies with depth. From equations (9) and (10), it can be deduced that where σ v02  is the in situ vertical effective stress in the thick homogeneous layer at depth (D f + B), where D f is the footing embedment depth. Values of q c2 and σ v02  were used to calculate the relevant G 0 from equation (10) and then E 0 for use in equation (8) (12) where G 0 and q c are in kPa. This correlation was used for the clay layer at the Alvin site.
In order to determine unknown HSSmall parameters (except for G 0 ref when the v s profile is known), the relative density D r (%) for each sand layer was determined first from the correlation given by Jamiolkowski et al. (1985): For sand layers γ uns = γ sat . For the clay layer, assuming a constant shear wave velocity with depth, m = 0 so that G 0 ref = G 0 , and G 0 was determined from equation (12). The constrained modulus E oed ref = E oed was obtained from the following correlation (Sanglerat 1992 After having determined parameters for numerical analyses of all soil types, each sand layer is modeled as if the entire foundation soil consisted of the material of that layer. The two footings were modeled with dimensions equal to those used in the load tests. The sizes of the equivalent base area square footings are within the range of values analyzed by Mayne et al. (2012). The footing material is reinforced concrete.
It was stated in the subsection on numerical modeling that Benz (2007) had to use preloading of the foundation sand in order to make it overconsolidated and match the numerical and measured data. Setting the pre-overburden pressure (POP) among the initial condition input data to the value of preloading used in the analysis has the same effect as loading and unloading the soil by that value prior to simulating the load test. POP has a constant value for a given soil, and it is given by (Brinkgreve et al. 2011) where σ p  is the preconsolidation pressure.
If POP is set to zero, the soil is normally consolidated, and plastic deformations start from the beginning of calculations. When POP is set to a nonzero value, according to the HSSmall model, there are elastic deformations at very small strains. Now, the analytical load-settlement curve according to equation (8) is calculated, and POP in Plaxis is adjusted so that the two curves fit the best. However, this step can be used only when values of G 0 ref for each layer are determined from shear wave velocities. Since this study uses no v s measurements, additional steps must be undertaken.
Whenever the value of POP is changed for a given homogeneous foundation sand, the horizontal effective stresses in the initial state also change. Therefore, for each new value of POP, the horizontal effective stress at depth (D f + B) is read off from the Plaxis presentation of the initial state. Then, G 0 ref for the current calculation is obtained from equation (9), where the value of G 0 is taken as the one previously determined at depth (D f + B) for the use of E 0 in the analytical equation (8).
When the fitting procedure is completed for all sand layers, and parameters are defined for the clay layer at the Alvin site, all parameters for the Plaxis analysis of the layered foundation soil, including POP for each sand layer, are set.

Green Cove Springs
At 1.8 m in diameter, 0.6 m in thickness, and 0.6 m embedded in the ground, reinforced concrete footing was load tested to 222 kPa in Green Cove Springs, Florida (Anderson et al. 2007). The groundwater table was at 1.7 m from the ground surface. The 1.8-m-diameter circular footing was taken as an equal area square footing with the size 1.6 m for the analytical expression, and it was numerically modeled as an axisymmetrical footing.
A description of the four foundation soil layers and their relevant parameters are given by Anderson et al. 2007, along with relative densities calculated from equation (13). For the fourth layer, equation (13) gives the relative density of 0. Thus, a value of 0.1 was adopted in order to calculate parameters for Plaxis. Fig. 5 shows the finite element mesh for the layered foundation soil. Strain-controlled calculations were carried out. Fig. 6 shows the measured data (Anderson et al. 2007) and all calculated results. The full black curve is the result of a Plaxis calculation of settlements for the layered foundation soil. It matches the measured data very well. The full gray curves are the results of settlement calculations by equation (8)    The Plaxis fit to equation (8) is very good for all four sand layers. This was achieved by selecting G 0 ref through the procedure described above in the absence of v s measurements. Then, POP was adjusted for the best fit. In this analyzed case, the advantage of using the modified  (1), but this is a simple example where a load test was carried out so that the final result could be compared to measured data for a layered foundation soil. The proposed methodology is convenient for the prediction of load-settlement curves in geotechnical problems where it is not plausible to use the simple Mayne et al. method. Furthermore, it provides direction for more efficient numerical modeling of shallow foundations.

Alvin (East)
Another simple illustrative example of the application of the proposed methodology is one of the four load tests performed in Alvin, Texas (Tand et al. 1994), which essentially covers two footings. Four belled reinforced concrete footings were load tested at two separate locations (East/West). The two footings at each site were installed 2.7 m apart using a crane-mounted drilling rig. The upper 2.7 m at each site is primarily stiff sandy clay (CL) to clay (CH). The clay is underlain by a thick deposit of medium dense to dense sand (SM/SP) to a depth greater than 15 m. The water table is at 1.9 m below the ground surface.
The East location was chosen in this study. The analysis presented herein covers both footings at the East location with a single model. The two footings at this location have an average base diameter of 1.95 m, which corresponds to an equal area square footing of B = 1.7 m. The footings are embedded at 2.2 m, which equals their thickness. Therefore, there is a clay layer 0.5 m thick below the footing bases. A description of the four foundation soil layers and their relevant parameters are given by Tand et al. 1994. The first layer composed of clay was numerically modeled in the undrained state with the undrained shear strength c u as the input parameter for strength characteristics ( u = 1º for numerical reasons). Fig. 7(a) shows the results of Plaxis fit to equation (8) for each sand layer, again by adjusting POP. In this case, the two curves for separate layers do not fit as well as in the previous case because the settlements were calculated up to 70 mm, whereas they were calculated up to 4 mm for the Green Cove Springs location. Fig. 7(b) shows the measured data (Tand et al. 1994) and final Plaxis calculation for the layered foundation soil (full black curve). This curve matches exceptionally well with the measured data. The Plaxis curves are reproduced for layers 2 and 3, and the Plaxis curve for clay is added for reference in Fig. 7(b). The depth of 2B = 4.4 m is within the third layer.
This example also shows that the approach suggested in the proposed methodology can provide good results. This is especially encouraging for the numerical modeling of shallow foundations because of the usual problems with the determination of numerous input parameters which, according to the proposed procedure, require CPT and shear wave velocity measurements in situ and only standard tests in the laboratory (e.g., direct shear and determination of basic soil characteristics). It is further encouraging because the problem of systematic overprediction of settlements on sand by numerical modeling can be solved by adjusting the POP or likewise the specification of overconsolidation in other computer programs against the analytical calculation given by equation (8) (8) and Plaxis best fit to equation (8), (b) measured data (Tand et al. 1994), Plaxis calculation for layered soil, and Plaxis calculation for first three layers separately

CONCLUSIONS
The nonlinear load-settlement curve for a simple footing on thick, fairly homogeneous sand stratum may best be predicted by empirical methods that correlate it directly to the results of some in situ tests. The method by Mayne et al. (2012) uses the cone resistance of the standard cone penetration test as a single soil parameter and appears to be the most simple and convenient solution. For more elaborate spread foundations or soil strata that cannot be treated as deep and homogeneous, this method is not applicable. Finite element or similar numerical methods, along with appropriate constitutive models for soils, which have the potential to model these geometrically more complex problems, have drawbacks. The constitutive model parameters used in the numerical analysis, if determined from unavoidably disturbed sand samples to be tested in the laboratory, are not reliable enough (Benz 2007). In addition, Benz introduced a hypothetical preload of the foundation sand to match the numerically calculated and measured settlements of an experimental footing.
To mitigate the described deficiencies, a methodology is proposed in the paper that combines the advantages of both the empirical and numerical approaches in predicting loadsettlement curves for shallow foundations. The proposed methodology may be used for layered foundation soil that also contains fine-grained soil layers, or when a reliable nonlinear soil-structure interaction analysis is desired. In the proposed methodology, it is suggested to determine the maximum shear modulus from measurements of the shear wave velocity in the field, and perform a CPT. Then, soil parameters for the numerical soil model for each layer from the soil profile separately, and comprising the entire foundation soil, can be determined from correlations with the CPT cone resistance. For each sand model, a loadsettlement curve is calculated using the proposed modified Mayne et al. (2012) method, which accounts for the finite soil stiffness at very small strains. The numerical load-settlement curve is fitted to match this calculated curve by adjusting the overconsolidation parameter.
The modified load-settlement curve has two parameters that were determined through a statistical analysis of 16  The proposed methodology was verified at two locations where load tests of simple footings resting on nonhomogeneous foundation soil were carried out: Green Cove Springs (Anderson et al. 2007) and Alvin (East) (Tand et al. 1994). The two examples were used to illustrate the potential of the methodology, even though neither location had measurements of the shear wave velocity. The shear modulus at very small strains was estimated from correlations with the CPT cone resistance. The final results at both locations showed very good matching between the measured settlements and those calculated numerically for a footing on layered soil. This is very promising because of problems encountered when soil parameters for a complex constitutive model have to be determined, and because of the systematic overprediction of settlements on sand by standard numerical modeling.