Figure 2 shows three different finite element meshes (i.e., coarse, fine, and optimized) that were generated to examine the effect of element size at the interface zone on the load transfer between the soil and facing structure. The coarse mesh (Figure 2A) had the highest element aspect ratio within the real interface zone; the optimized mesh (Figure 2C) had the lowest element aspect ratio in the region where the analysis is focused (and fewer total number of elements), and the element aspect ratio of the fine mesh (Figure 2B) was between that of the coarse and optimized meshes. When zero-thickness elements were used at the interface between the soil and facing, the interface virtual thickness was 18 mm as mentioned earlier. When using continuum elements to simulate the soil-facing interaction, the same 18 mm-value of real zone thickness was modeled (see Figure 2—Details).

Figure 3 shows the normal and shear stresses acting at the interface between the facing and backfill soil with the three different meshes and using the strength/stiffness reduction factor of R _i  = 0.8 where both zero-thickness elements and continuum elements are considered. The use of R _i  = 0.8 results in an interface friction angle of approximately 38°, which is similar to the measured friction angle between smooth concrete and sand (Potyondy 1961; Samtani and Nowatzki 2006). The numerical modeling showed that for the cases examined with zero-thickness elements, the finite element mesh had a minor effect on the normal and shear stresses at the interface between the soil and facing. However, when using continuum elements, both interface normal and shear stresses fluctuated once soil plastic flow occurred for all three meshes. The results also showed that the optimized mesh with the lowest interface continuum element aspect ratio experienced the smallest stress fluctuation amplitudes.

The total horizontal and vertical forces acting at the interface are shown in Supplementary Tables S2, S3, respectively. Both total horizontal and vertical forces using zero-thickness elements are in good agreement with results using continuum elements, and the finite element mesh had a minor effect on the total horizontal and vertical forces for both zero-thickness elements and continuum elements.

Figure 4 shows the normal and shear stresses at the interface between the facing structure and backfill soil for the three different strength/stiffness reduction factors investigated. The modeling results showed that for the continuum elements, increasing the strength/stiffness reduction factor (i.e., increasing the interface stiffness) resulted in greater amplitude of both normal and shear stress fluctuations in the plastic region when other conditions were equal. However, as presented in Supplementary Table S3 in the Supplemental Material, the total vertical loads (i.e., equivalent force from shear stresses) at the interface between the facing and backfill soil from the continuum elements are in good agreement with those from simulations with zero-thickness elements.

Continuum element interfaces were used in CODE_BRIGHT to simulate the soil-facing interaction with 18-mm-thick real zone. As in the previous cases, the structure (facing concrete panel) was modeled as linear elastic with elastic modulus of 32 GPa, Poisson ratio of 0.15, and a unit weight of 25 kN/m³. Soil material was modeled with Drucker–Prager failure criterion (i.e., circular cone in the principal stress space). In the program, the soil friction angle (ϕ) is defined with the critical state slope M-line in order to obtain correct shear strength. Using the Mohr–Coulomb failure criteria that circumscribe the Drucker–Prager failure criteria (see Supplementary Figure S3 in the Supplemental Material), strength is defined as follows: M = 6 ⁡ sin ⁡ ϕ 3 − sin ⁡ ϕ   ( triaxial   compression   state ) (14)

For triaxial compression behavior of material, M _compression, that is, σ ₁ > σ ₂ = σ ₃, where σ ₁, σ ₂, and σ ₃ are the main stress states: σ ₁ (major stress) ≥ σ ₂ ≥ σ ₃ (minor stress)). Assuming Drucker–Prager circle failure criteria that are inscribed within the Mohr–Coulomb failure criteria (see Supplementary Figure S3), the strength is defined as follows: M = 6 ⁡ sin ⁡ ϕ 3 + sin ⁡ ϕ   ( triaxial   extension   state ) (15)

For triaxial extension behavior of material, M _extension, that is, σ ₁ = σ ₂ > σ ₃. Thus, the stress states using the Mohr–Coulomb failure criteria can be defined by the M-parameter corresponding to the compression or extension condition.

As noted earlier, proper material constitutive modeling assuming M-values requires to determine the material state in terms of intermedia main stress scenario (i.e., σ ₂ = σ ₃ in compression, otherwise extension if σ ₂ = σ ₁). In common scenarios and in the ones assumed in the current study, the compression state corresponds to the most suitable case. Thus, soil model properties are the ones presented in Table 1, with proper choice of the compression M-parameter value: M _compression = 1.808.

Table 2 shows the same parameters as in Supplementary Table S1 but including M-values for compression, extension, and intermedia (M _average) case states. While results demonstrated that the interface material M-value falls within compression and extension values, default M _compression values were assumed for modeling the soil-facing interaction in this study.

Four different finite element 2D meshes were generated in CODE_BRIGHT: unstructured or irregular (but optimized) mesh with linear–triangular elements (Figure 7A), linear–triangular structured-mesh (Figure 7B), bilinear–quadrilateral structured-fine-mesh (Figure 7C), and bilinear–quadrilateral structured-coarse-mesh (Figure 7D). With these elements, it was possible to examine the effect of element size at the interface zone on the load transfer between the soil and facing structure.

Among these three meshes, the irregular mesh presented in Figure 7A had the highest element aspect ratio far from the analyzed soil–structure zone, but the optimized shape becomes finer with smaller aspect ratio in the interface zone. The structured fine mesh presented in Figure 7B has fine definition in all regions. However, triangular linear elements with analytical integration may not be the best ones because of the simple linear interpolation definition of this element type. This type of element did not work very well when shear strains occur with limits on volumetric strain development. Figures 7B,C have same number of nodes, even though the triangular mesh scheme has double the number of elements. However, both Figures 7C,D are for quadrilateral structured elements, which implies bilinear interpolation and numerical integration with quadrature of 4 Gauss points, and are probably better to perform soil–structure interactions for regular shear strain scenarios (despite the larger computational efforts required).

Settlement resulting under q = 100-kPa surcharge for the different meshes previously presented is shown in Figure 8 (R _i  = 0.8). Despite the different meshing element types, good agreement was obtained between the four alternatives, for example, a maximum settlement value of 2.33 cm using triangular linear elements (both optimized and fine mesh cases) and 2.32 cm using quadrilateral bilinear elements.

Supplementary Figure S4 in the Supplemental Material presents the total shear strain evolution under q = 50 and 100-kPa surcharges (R _i  = 0.8). Because results are scaled from 100-kPa surcharge in each meshing case, similar color distribution is obtained between meshing type with the exception of the quadrilateral coarse mesh (greater zone affected due to elements size). However, despite the similar responses in soil settlement with different mesh types used (Figure 8), a difference in strains from approximately 13% using the structured triangle mesh up to 19% using the structured quadrilateral fine mesh was obtained.

Results for load transfer from backfill soil to facing panel are plotted in Figure 9 for two cases of interface strength/stiffness reduction factors (R _i  = 0.8 and R _i  = 0.3). Despite variations noted in previous results, the triangular structured mesh improves results by avoiding the fluctuations in both normal and shear stress results when unstructured (but optimized) meshes are considered. Using quadrilateral elements, small differences were obtained between meshing size cases. These differences were even smaller when a softer (and weaker) interface reduction factor was selected.

Figure 10 presents the normal and shear load transfer from backfill soil to facing panel for different critical state slope M-parameter defining interface strength (see explanation in Section 3.2.1). Different results were obtained because of the M-value selected based on the triaxial state. Triaxial compression and extension state trends (i.e., M _compression and M _extension) generate a region of possible results. Differences between both compression and extension states are, however, not dramatic (e.g., similar differences were obtained using FLAC with the interface modeled by springs or continuum material; Figure 5). Results using an average value of M are also plotted (M _average) and these falls between both boundary triaxial states. The modeling results for compression and extension triaxial states using the continuum interface in CODE_BRIGHT and the continuum interface in PLAXIS and FLAC are presented in Figure 11 (R _i  = 0.8 and R _i  = 0.3 cases). Despite the fluctuations in the PLAXIS model response, similar results were obtained for FLAC and PLAXIS models used by Yu et al. (2014,2015). CODE_BRIGHT, PLAXIS, and FLAC results are shown to fall within the region bounded by compression and extension triaxial M-values results.

Complementary to Supplementary Figure S4, Figure 12 shows the normal and shear stresses at the interface between the facing structure and backfill soil for three different strength/stiffness reduction factors investigated in CODE_BRIGHT modeling with quadrilateral coarse mesh elements (mesh case shown in Figure 7D). As before, the modeling results showed that increasing the strength/stiffness reduction factor R _i (i.e., increasing the interface stiffness) resulted in a smaller shear strains and relative displacements in the affected region: from approximately 13% for the rigid interface case to 22% for the R _i  = 0.3 case. Note that, despite that the maximum displacement was similar at the right boundary contour (free displacement condition), the displacement distributions change markedly using the three strength/stiffness interaction factors.

Results for the normal and shear stress transfer from backfill soil to facing panel using different interface strength/stiffness reduction factor and 100- and 50-kPa surcharge loading are found in Figure 13. Despite that the stress magnitude was different between the q-loading cases analyzed, very similar distributions were obtained. Significant variations in these distributions were obtained using the R _i values considered. Different R _i values influenced the resulting soil plastic zone (i.e., peak shear stress value location).

Modeling of an 18-mm-thick (t _i ) interface zone between dissimilar materials in full-height earth-retaining walls using continuum elements can be problematic because of the large difference in shape and size geometry between the different components within the retaining wall. It may be necessary to increase the interface thickness to accommodate other domain geometries. If this is the case, then the properties of the thicker interface must be adjusted to maintain the same (or similar) normal and shear strength/stiffness.

In this section, two interface thicknesses of t _i ^*  = 50 and 100 mm were examined using quadrilateral structured-coarse mesh case (Supplementary Figure S5 in the Supplemental Material, which correspond to complementary cases from Figure 7B—t _i  = 18-mm-thick case). To keep the same interface stiffness, the new shear modulus (G _i ^* ) of the interface was calculated as follows: G i ∗ = ( G i / t i ) t i ∗ (16) where the Poisson ratio is the same for both interface thickness cases (i.e., v _i ^*  = v _i  = 0.45), and the new oedometer modulus (E _oed,i ^* ) and elastic modulus (E _i ^* ) can be calculated using Eqs 4, 5, respectively, with the new shear modulus (G _i ^* ) and Poisson ratio (v _i ^*  = 0.45). Table 3 presents the equivalent interface properties for the case studies assumed and the additional interface thickness cases examined (as before, t _i ^*  = 50 and 100 mm).

The computed settlements and shear strain (deviatoric invariant) for both interface thicknesses assumed (t _i  = 0.05 m and t _i  = 0.10 m) are presented in the Supplemental Material in Supplementary Figure S6 (case with R _i  = 0.8) and Supplementary Figure S7 (R _i  = 0.3). As shown, very similar responses were obtained for the three interface thickness cases (t _i  = 18 mm case results presented in Figure 8D, Supplementary Figure S4 for R _i  = 0.8, and Figure 12 for R _i  = 0.3) under q = 50- and 100-kPa surcharge scenarios. As shown for previous cases, the reduction of the interface strength/stiffness interaction factor leads to an increase in the affected shear strain localization zone.

Results of normal and shear stress transfer from the backfill soil to the facing panel are presented in Figure 14 for the three interface thickness cases under q = 10-, 50-, and 100-kPa surcharge cases, and R _i  = 0.8 and R _i  = 0.3 interface strength/stiffness interaction cases. Results were obtained from a cross-section location at 9-mm distance from facing panel (i.e., at x = 0.509-m distance from left boundary), which corresponds to a line located within the t _i  = 18-mm interface thickness case. The data show very similar responses for the three interface thickness cases.

Supplementary Figure S8 in Supplemental Material presents the same previous results but for the cross-section location in the middle of each interface thickness: at 9-mm distance for ti = 18-mm-thick case, at 25-mm distance for ti = 50-mm-thick case, and at 50-mm distance for ti = 100-mm-thick case. It can be seen that there are practically no differences from the previous fixed cross-section location results.

The numerical results demonstrate that the increased interface thickness cases had a minor effect on the total vertical load at the interface between the facing and backfill soil if the equivalent interface stiffness was kept the same. Thus, a real interface zone between the dissimilar materials using continuum elements with a thickness greater than the virtual interface thickness using zero-thickness elements can be generated to model the soil–structure interactions and give similar numerical outcomes if the soil property values within the real interface zone are properly calculated based on the same interface stiffness.