## 2 Problem Definition
For further details regarding the theoretical background and the numerical implementation  
Figure 1: Problem geometry and notation. of this method, readers are referred to [49; 50; 51] for the lower bound finite element limit theory, and to [49; 51; 52; 53] for the upper bound one. Notably, a striking advantage of the FELA approach is that preliminary assumptions related to the failure surface are not needed [54]. Meshes consisting of triangular elements have been adopted, independently of the values of the geometrical parameters introduced above. To start the iterative mesh adaption procedure, each model has been assumed to be made of 1000 elements. Through h-adaptivity, elements have been automatically generated or reduced in size where shear stress concentration and dissipation is going to affect the solution, until a limit of 10,000 elements has been attained. Such initial and final numbers of elements have been chosen on the basis of the results presented in [55; 56], the latter to specifically keep the computational burden of each analysis below a critical threshold. The Mohr-Coulomb failure criterion has been used to model the soil behavior, while the strip footing has been modeled as a weightless rigid plate. The ultimate load bearing capacity has then been defined through a multiplier for a uniform load applied to the footing, and automatically increased until the formation of collapse mechanisms. In the analyses, a strip footing with a width \(B=1\) m has been placed on soil with a specific weight of 20 kN/m\({}^{3}\). To avoid any perturbation to the solution induced by the domain boundaries, a model with a width of \(40B\) and a depth of \(20B\) has been considered. The contact of the soil with the bedrock has been modelled by constraining both the horizontal and the vertical displacements along the bottom surface. Along the side boundaries, only the horizontal components of the displacements have instead been restrained. The pseudo-static method has been adopted to assess the seismic performance of the foundation over the cavity. For this purpose, the value of the horizontal acceleration coefficient \(k_{h}\) has been exploited, and once the upper and lower bounds on \(q_{u}\) or \(Q\) at collapse have been computed, the relevant mean value has been assumed as the representative result for each model. An example of the adopted models is depicted in Figure 2, where the typical FELA adaptive mesh of the Optum G2 software, along with the adopted boundary conditions, are sketched. To validate the modeling approach, the results of some studies in the literature have been considered for comparison purposes. First, as seen in [57], the bearing capacity of the foundation has been computed via the lower bound FELA method, in case of horizontal seismic loadings. The relevant results are compared in Figure 3, in terms of the bearing capacity factor \(N_{\gamma}\), to show how they perfectly match the available ones. Second, a comparison with the studies of [57; 58; 59; 60; 61; 62; 63; 64] has been made through the data gathered in Table 2. Such a comparison is provided in terms of the bearing capacity factor, at a varying value of the internal friction angle \(\varphi\). The results again show a good agreement with those  
Figure 2: Typical FELA mesh for the studied problem with a close-up to show the adaptive mesh refinement around the void, where the collapse mechanism is expected to form. published, and correctly bound almost all of them. Third, in Table 3, a further comparison is provided for the bearing capacity factor \(N_{c}\), allowing for the presence of a void at different depths \(D\) beneath the strip footing central line, see [10; 12; 19]. It has been assumed that the opening width is equal to the surface footing width, and is located within an utterly cohesive soil at an undrained state with cohesion varying in the range 60-300 kPa. Results are shown again to agree well with the reference ones.