Assessment of the Bearing Capacity of Foundations on Rock Masses Subjected to Seismic and Seepage Loads

Rubén Galindo 1,* , Ana Alencar 1, Nihat Sinan Isik 2 and Claudio Olalla Marañón 1 1 Departamento de Ingeniería y Morfología del Terreno, Universidad Politécnica de Madrid, 28040 Madrid, Spain; at.santos@alumnos.upm.es (A.A.); claudio.olalla@upm.es (C.O.M.) 2 Department of Civil Engineering, Faculty of Technology, Gazi University, 06560 Ankara, Turkey; nihatsinan@gazi.edu.tr * Correspondence: rubenangel.galindo@upm.es


Introduction
Natural events such as earthquakes and floods generate impulsive loads and seepages in the ground, endangering civil constructions. Infrastructures must be built to resist natural risks in adequate safety conditions and be economically viable. Therefore, an optimized design that reliably considers the negative effects that these natural phenomena introduce to the foundations is highly advisable to guarantee the sustainability of the infrastructures.
Today, research in the area of seismic bearing capacity is very much in demand because of the devastating effects of foundations under earthquake conditions. A high number of failures have occurred where field conditions have indicated that the bearing capacity was reduced during seismic events.
In the estimation of the bearing capacity, when the effect of the earthquake is considered, it is usual to adopt the pseudo-static hypothesis where the seismic force acts as an additional body force within the soil mass. The vertical and horizontal acceleration are applied both on the ground and in the structure. Thus, the limit conditions can be evaluated by introducing pseudo-static equivalent forces, corresponding with the inertial forces in the soil during the seismic excitation. Such an approach is based on the hypothesis of a synchronous motion for the soil underneath the footing, a hypothesis that is acceptable only in the case of small footing widths and large values for the soil stiffness. According to mass with a homogeneous and isotropic behavior. The linearization of the failure criterion implies incorporating approximate methods and requires iterative procedures to ensure an optimized upper or lower bound for the solution; therefore, it is desirable to be able to address the problem using the non-linear criterion directly.
In addition, most of the formulations are limited to flat ground, the need to analyze the bearing capacity of shallow foundations on sloping ground of moderate slopes being very common in dam and bridge foundations. Finally, although the numerical solutions allow solving complex problems with singular considerations when a seismic load or filtration acts, in the face of non-linear criteria, a complete analysis of numerical convergence is necessary, which complicates the practical applicability and the rapid design of foundations in rock masses.
An analytical method for the calculation of shallow foundations that solves the internal equilibrium equations and boundary conditions combined with the failure criterion was proposed by Serrano and Olalla [22] and Serrano et al. [23], applying the Hoek and Brown [19] and the Modified Hoek and Brown failure criterion [20], respectively. It is based on the characteristic lines method [24], with the hypothesis of the weightless rock, strip foundation and associative flow law. The formulation of the bearing capacity proposed by Serrano et al. [23] introduces a bearing capacity factor, which makes the failure pressure proportional to the uniaxial compressive strength of the rock (UCS).
In the present study, the analytical formulation of Serrano et al. [23] and design charts were developed to study the bearing capacity when there is an increase in load induced by forces of seismic origin or filtration in rock masses, where it is necessary to use a non-linear failure criterion and the own weight of the ground is generally negligible compared to the resistant components, considering, as usual, the possibility of the inclination of the ground to the edge of the foundation. Besides, a numerical model was created through a finite difference method, assuming a similar hypothesis for the analytical solution, and it was observed that the results obtained by both methods were quite similar.

Mathematical Model
As is generally known, in rock mechanics, the non-linear Modified Hoek and Brown failure criterion is the most used, and it is applicable for a rock mass with a homogeneous and isotropic behavior, meaning that by the inexistence or by the abundance of discontinuities, it has the same physical properties in all directions.
In this research, the Modified Hoek and Brown failure criterion [21,25] was used, and it was formulated as a function of the major principal stress (σ 3 ) and minor principal stress (σ 1 ) according to the following equation: The uniaxial compressive strength (UCS) is σ c , while the parameters m, s and a can be evaluated with (2)-(4) and depend on the intact rock parameter (m o ), quality index of the rock mass (geological strength index (GSI)) and damage in the rock mass due to human actions (D), which in shallow foundations, is usually equal to zero.
Serrano et al. [23] proposed an analytical formulation for estimating the ultimate bearing capacity of the strip footing for a weightless rock mass, based on the characteristic method, which allows solving the internal equilibrium equations in a continuous medium together with the boundary equations and Sustainability 2020, 12, 10063 4 of 21 those that define the failure criterion. This solution is based on the Modified Hoek and Brown failure criterion [20], taking into account the associated plastic flow rule.
According to this analytical formulation, the ground surface that supports the foundation is composed of two sectors ( Figure 1): Boundary 1 (free), with the inclination i 1 , where the load acting on a surface is known (for example, the self-weight load on the foundation level or the load from installed anchors), and Boundary 2 (foundation), where the bearing capacity of the foundation should be determined (acting with the inclination of i 2 ). Serrano et al. [23] proposed an analytical formulation for estimating the ultimate bearing capacity of the strip footing for a weightless rock mass, based on the characteristic method, which allows solving the internal equilibrium equations in a continuous medium together with the boundary equations and those that define the failure criterion. This solution is based on the Modified Hoek and Brown failure criterion [20], taking into account the associated plastic flow rule.
According to this analytical formulation, the ground surface that supports the foundation is composed of two sectors ( Figure 1): Boundary 1 (free), with the inclination i1, where the load acting on a surface is known (for example, the self-weight load on the foundation level or the load from installed anchors), and Boundary 2 (foundation), where the bearing capacity of the foundation should be determined (acting with the inclination of i2). The solution based on the characteristic lines method requires the equation of the Riemann invariants (Ia) [26] fulfilled along the characteristic line: In this equation, the instantaneous friction angle at the boundary 2 (ρ2) is the only unknown, because the other variables can be defined at Boundary 1: the instantaneous friction angle at The solution based on the characteristic lines method requires the equation of the Riemann invariants (I a ) [26] fulfilled along the characteristic line: In this equation, the instantaneous friction angle at the boundary 2 (ρ 2 ) is the only unknown, because the other variables can be defined at Boundary 1: the instantaneous friction angle at Boundary 1 (ρ 1 ) and the angle (Ψ 1 ) between the major principal stress and the vertical axis in this sector ( Figure 1). Thus, expressing Ψ 2 (the angle between the major principal stress and the vertical axis in Boundary 2, as indicated in Figure 1) as a function of ρ 2 , it is possible to estimate the ultimate bearing capacity.
Through the analytical method [23], the bearing capacity was obtained by (7).
The resistant parameters β a and ζ a were applied to make dimensionless the calculation of the Modified Hoek and Brown failure criterion. β a represents the characteristic strength, which has the same units as the UCS and was used to make the pressures dimensionless, while ζ a (the "tenacity coefficient") is a dimensionless coefficient that, multiplied by β a , corresponds to the tensile strength.
A a , k and the exponent a are constants for the rock mass and depend on the rock type (m), UCS and GSI.
N β is the bearing capacity factor, and it can be calculated, according to the problem statement, as follows.
The angle of internal friction (ρ 1 ) can be obtained by iteration from the load at Boundary 1. From the value of ρ 1 and by the iteration of (5), the value of the internal friction angle at Boundary 2 (ρ 2 ) can be calculated.
Finally, knowing ρ 2 , the bearing capacity factor (N β ) can be calculated, and using, again, parameters β a and ζ a , the ultimate bearing capacity (P h ) was estimated as an expression that depended on the instantaneous friction angle at Boundary 2 (ρ 2 ), the inclination of the load on the foundation (i 2 ) and the exponent of the Modified Hoek and Brown criterion (a; k = (1 − a)/a): cos(i 2 )

Consideration of Pseudo-Static Load: Mathematical Transformation
In the pseudo-static approach, static horizontal and vertical inertial forces, which are intended to represent the destabilizing effects of the earthquake or seepage, are calculated as the product of the seismic/seepage coefficients and the distributed load applied to the boundaries. In the case of the rock mass, the weight collaboration is usually negligible compared to the resistance of the ground, and therefore, the inertial forces are applied both to the foundation and to the free boundary.
The vertical seismic/seepage coefficient k v is supposed to be a fraction of one horizontal k h , and in particular, the vertical acceleration is thus assumed to be in phase with the horizontal acceleration.
The present study is divided into three parts: (a) The first one considered the horizontal (k h ) and vertical (k v ) components of the pseudo-static load on both boundaries, with a free boundary inclined by α at the edge of the foundation, which resembles the hypothesis of a seism. (b) In the second part, only the horizontal component (k h ) on the foundation boundary was adopted (it being possible to consider both the horizontal and vertical components on the free boundary depending on the direction of the seepage), with the free boundary inclined by α at the edge of the foundation. This hypothesis is more similar to the presence of a seepage (both hypotheses are represented schematically in Figure 2, and they are solved and shown in new charts including additional horizontal and vertical loads). (c) The final section is the comparison of the analytical result with that obtained numerically through the finite difference method. The application of the analytical method [23] with an increase in the horizontal and/or vertical loads produced by inertial forces can be carried out by means of a parametric transformation from the incidence angles of the loads in the static configuration to the final configuration including the seismic or seepage loads. Thus, for a general case of a pseudo-static force on the two boundaries, the starting point is the load acting in the static hypothesis (subscript 0), as is represented in Figure 3a. In this case, the inclinations of the load on the foundation (p) and of the load on the free boundary (q) are 02 and 01 , respectively. However, considering the pseudo-static load, the inclinations of the loads in both boundaries are different (i2 for the foundation boundary and i1 for the free boundary). Figure 3b allows the deduction of the mathematical transformations of these angles from the static configuration to the final pseudo-static configuration as a function of the angle (α) of the free boundary and of the horizontal (kh1 or kh2, depending on the boundary) and vertical (kv1 or kv2, depending on the boundary) components of the pseudo-static load. These transformations are expressed in (10) and (11).
Boundary 2 (Foundation) Boundary 1 (Free) The application of the analytical method [23] with an increase in the horizontal and/or vertical loads produced by inertial forces can be carried out by means of a parametric transformation from the incidence angles of the loads in the static configuration to the final configuration including the seismic or seepage loads. Thus, for a general case of a pseudo-static force on the two boundaries, the starting point is the load acting in the static hypothesis (subscript 0), as is represented in Figure 3a. In this case, the inclinations of the load on the foundation (p) and of the load on the free boundary (q) are i 02 and i 01 , respectively. However, considering the pseudo-static load, the inclinations of the loads in both boundaries are different (i 2 for the foundation boundary and i 1 for the free boundary). Figure 3b allows the deduction of the mathematical transformations of these angles from the static configuration to the final pseudo-static configuration as a function of the angle (α) of the free boundary and of the horizontal (k h1 or k h2 , depending on the boundary) and vertical (k v1 or k v2 , depending on the boundary) components of the pseudo-static load. These transformations are expressed in (10) and (11).
Therefore, for the pseudo-static calculation, the analytical formulation of the characteristics method expressed by (5) can be used using the load inclinations on each boundary obtained through the transformations indicated in Equations (10) and (11) as a function of the horizontal and vertical components of the added inertial load and of the angle of inclination of Boundary 1.
In the case of the seismic load, it is considered that k h1 = k h2 and k v1 = k v2 , while in the case of the seepage load, k v1 = 0; thus, for clearer notation, it is denoted that for the seismic load, the horizontal and vertical components of the pseudo-static load are k h and k v , respectively, and for the seepage load, the horizontal component of the pseudo-static load on the foundation boundary is called i a . expressed in (10) and (11).  Therefore, for the pseudo-static calculation, the analytical formulation of the characteristics method expressed by (5) can be used using the load inclinations on each boundary obtained through the transformations indicated in Equations (10) and (11) as a function of the horizontal and vertical components of the added inertial load and of the angle of inclination of Boundary 1.
In the case of the seismic load, it is considered that kh1 = kh2 and kv1 = kv2, while in the case of the seepage load, kv1 = 0; thus, for clearer notation, it is denoted that for the seismic load, the horizontal and vertical components of the pseudo-static load are kh and kv, respectively, and for the seepage load, the horizontal component of the pseudo-static load on the foundation boundary is called ia.

Calculation Cases and Representation of Analytical Results
Once the mathematical transformation of the load angles in the boundaries, for the problem presented in Figure 3a, according to (10) and (11) for the seismic and seepage load has been carried out, the analytical formulation using the method of the characteristic lines can be applied. The results are presented as graphs, which allow the estimation of the bearing capacity considering the presence of a pseudo-static load.
The charts are clustered according to the exponent "a" of the Modified Hoek and Brown criterion, the inclination α of Boundary 1 and kv of the foundation boundary, and they were developed based on io2 and the horizontal component of the pseudo-static load of the foundation boundary (kh). It is noted that for high confining pressures in Boundary 1, it is not always possible to obtain a bearing capacity value; this limit is demarcated by the non-equilibrium line. Figure 3. Scheme of the seismic load estimation: (a) static configuration and (b) pseudo-static inclination on Boundary 1 (external boundary of foundation). Note: In this figure, the subscripts "v" and "h" refer to the vertical and horizontal projections of the load.

Calculation Cases and Representation of Analytical Results
Once the mathematical transformation of the load angles in the boundaries, for the problem presented in Figure 3a, according to (10) and (11) for the seismic and seepage load has been carried out, the analytical formulation using the method of the characteristic lines can be applied. The results are presented as graphs, which allow the estimation of the bearing capacity considering the presence of a pseudo-static load.
The charts are clustered according to the exponent "a" of the Modified Hoek and Brown criterion, the inclination α of Boundary 1 and k v of the foundation boundary, and they were developed based on i o2 and the horizontal component of the pseudo-static load of the foundation boundary (k h ). It is noted that for high confining pressures in Boundary 1, it is not always possible to obtain a bearing capacity value; this limit is demarcated by the non-equilibrium line.
For the development of the new charts, three values of the GSI (geological strength index) were adopted (8,20 and 100), which generated exponents "a" of the Modified Hoek and Brown criterion equal to 0.5, 0.55 and 0.6 from (4). Based on (7), the values of the rock type (m o ) and the uniaxial In the graphs developed for k v > 0 on the foundation boundary, representative of the seismic load, four k h values were adopted and correlated two values of k v , two slope angles for the free boundary (α) and three initial inclination angles for the load on the foundation boundary (i o2 ), representing the inclination angle of the load without considering the pseudo-static load. In Table 1, the values used in the analysis are indicated. Table 1. Geometric parameters adopted in the model (k v > 0). On the other hand, in the charts developed with k v = 0 on the foundation boundary, representative of the seepage loads, three values of horizontal load were used, in those cases called i a (additional inclination); three values of the slope (α) and another three values of i o2 were also used, which are shown in Table 2. Table 2. Geometric parameters adopted in the model (k v = 0).  Tables 1 and 2 are the relation between the horizontal and the vertical load. The final inclination of the load in the foundation boundary, in the cases of seepage loads, can be obtained directly as a simplification of (10):

Charts for Estimation of Bearing Capacity Considering k h and k v > 0 (Seismic Case)
These design charts allow obtaining the bearing capacity factor N β (9), and they are presented in Figures 4-6. Each graph represented corresponds to determined values of the exponent "a" of the Modified Hoek and Brown criterion of the rock mass (a function of the GSI of the quality of the rock mass), the inclination α of the free boundary and the ratio of the vertical and horizontal components of the pseudo-static load (k v /k h = 1 or 0.5). In each graph, different curves corresponding to the angles of the inclination of the static load on Boundary 2 (i o2 ) and horizontal component of the pseudo-static load (k h ) are presented, so that in the abscissa, the known value of the normalized main major stress is presented, normalized on Boundary 1 (σ* 01 ), estimated through (13). It is dimensionless and corresponds to the load acting on Boundary 1, and its value depends on the inclination angle of the load i 1 obtained by (11) (Figure 1).
Among the graphs, it is easy to appreciate that as there is a greater slope for the free boundary (α), a lower k v /k h ratio and a lower exponent "a" from the Modified Hoek and Brown failure criterion, the value of N β follows a declining trend.  9 Non-equilibrium boundary         Table 3 shows the results of the bearing capacity (P h ) of some studied cases. It should be noted that the value of P h is not directly proportional to N β , meaning that a higher value of N β does not necessarily mean that the value of P h will be higher. Besides, the same value of N β is associated with different values of P h depending on the other geotechnical parameters, as shown in (7). On the other hand, it is observed that under equal conditions (only varying the GSI and, consequently, the parameter "a"), the greater the GSI, the lower the value of N β ; however, as expected, the bearing capacity of the rock mass is higher. It should be noted that in each graph, an area is indicated, in the lower right part, in which it was not possible to obtain the mechanical balance due to excessively high load conditions on the free boundary.

Charts for Estimation of Bearing Capacity Considering Only a Horizontal Pseudo-Static Load on Foundation (Seepage Case)
The same representations were realized as in the previous section, where the load capacity factor N β of the analytical Equation (9) can be obtained in Figures 7-9. In this case, each graph represented corresponds to a determined value of the exponent "a", the inclination α of Boundary 1 and considering a value of k v = 0 on the foundation boundary. In the same way, in each graph, the different curves correspond to different inclination angles for the load on Boundary 2 (i o2 ) and the horizontal component of the seepage pseudo-static load (i a ). In this case, of the seepage load, the abscissa of the normalized main major stress normalized on Boundary 1 (σ* 01 ) estimated through (13) corresponds to the transformation of the original inclination of the load from the static configuration to the pseudo-static situation, where the horizontal and vertical components will appear according to the seepage trajectories considered.     In this case, since the seepage is typical of dam foundations, it is more realistic to consider cases of moderate ground slopes, showing the charts for slopes of the free boundary (α) equal to 0 • , 5 • and 10 • .

Numerical Validation
In order to compare the results obtained by applying the chart proposed in the present study with those numerically estimated (through FDM: finite difference method), the same hypotheses (weightless rock, strip foundation and associative flow law) were used in the rock-foundation model.
The 2D models were used to calculate the cases by the finite difference method, employing the commercial FLAC software [27], applying the plane strain condition (strip footing). Two-dimensional numerical models have been used by many researchers to solve problems of foundations under dynamic loads [28] and in rock masses [29]. Figure 10 shows the model used and where the boundaries were located, at a distance that did not interfere with the result; in all the simulations, the associative flow rule, weightless rock mass and smooth interface at the base of the foundation were adopted. The modified Hoek-Brown constitutive model available in FLAC V.7 was used, which corresponds to (1).

Numerical Validation
In order to compare the results obtained by applying the chart proposed in the present study with those numerically estimated (through FDM: finite difference method), the same hypotheses (weightless rock, strip foundation and associative flow law) were used in the rock-foundation model.
The 2D models were used to calculate the cases by the finite difference method, employing the commercial FLAC software [27], applying the plane strain condition (strip footing). Two-dimensional numerical models have been used by many researchers to solve problems of foundations under dynamic loads [28] and in rock masses [29]. Figure 10 shows the model used and where the boundaries were located, at a distance that did not interfere with the result; in all the simulations, the associative flow rule, weightless rock mass and smooth interface at the base of the foundation were adopted. The modified Hoek-Brown constitutive model available in FLAC V.7 was used, which corresponds to (1). It is assumed that the bearing capacity is reached when the continuous medium does not support more load because an internal failure mechanism is formed. In the case under study, due to the presence of a vertical and a horizontal force, the vertical force was considered unknown; therefore, in the calculation, a constant horizontal load was applied, while the vertical load increased until reaching failure. Thus, the inclination of the load applied was also unknown, because it depended on the ratio between the vertical and the horizontal components.
Therefore, to estimate the bearing capacity for a determinate load inclination, it is necessary to carry out a series of calculations to find the corresponding combination of the horizontal (σh) and vertical (σv) components. Figure 11 shows the results for vertical loads obtained for the case studied (kv = 0, io2 = 0°, mo = 15, UCS = 100 MPa, GSI = 65 and foundation width B = 2.25 m) as a function of the equivalent load inclination. It is assumed that the bearing capacity is reached when the continuous medium does not support more load because an internal failure mechanism is formed. In the case under study, due to the presence of a vertical and a horizontal force, the vertical force was considered unknown; therefore, in the calculation, a constant horizontal load was applied, while the vertical load increased until reaching failure. Thus, the inclination of the load applied was also unknown, because it depended on the ratio between the vertical and the horizontal components.
Therefore, to estimate the bearing capacity for a determinate load inclination, it is necessary to carry out a series of calculations to find the corresponding combination of the horizontal (σ h ) and vertical (σ v ) components. Figure 11 shows the results for vertical loads obtained for the case studied The vertical load was applied through velocity increments, and the bearing capacity was determined from the relation between the stresses and displacements of one of the nodes; in this case, the central node of the foundation was considered. In Figure 12a, the displacement of the central node of the foundation (abscissa) with respect to the load applied to the ground from the foundation is represented. In this figure is observed that the maximum load that the ground supports is limited to the asymptotic value of the curve represented.
Additionally, a convergence study was carried out, consisting of the analysis of the values of the bearing capacity obtained under the different increments of the velocity used. With a decrease in the value of the velocity increments, the result converged towards the final value that is the upper limit in the theoretical method (27b). The vertical load was applied through velocity increments, and the bearing capacity was determined from the relation between the stresses and displacements of one of the nodes; in this case, the central node of the foundation was considered. In Figure 12a, the displacement of the central node of the foundation (abscissa) with respect to the load applied to the ground from the foundation is represented. In this figure is observed that the maximum load that the ground supports is limited to the asymptotic value of the curve represented. The vertical load was applied through velocity increments, and the bearing capacity was determined from the relation between the stresses and displacements of one of the nodes; in this case, the central node of the foundation was considered. In Figure 12a, the displacement of the central node of the foundation (abscissa) with respect to the load applied to the ground from the foundation is represented. In this figure is observed that the maximum load that the ground supports is limited to the asymptotic value of the curve represented.
Additionally, a convergence study was carried out, consisting of the analysis of the values of the bearing capacity obtained under the different increments of the velocity used. With a decrease in the value of the velocity increments, the result converged towards the final value that is the upper limit in the theoretical method (27b). In the example studied, a seepage case was studied, where kv = 0, io2 = 0°, mo = 15, UCS = 100 MPa, GSI = 65, α = 0° and B = 2.25 m were adopted (which did not influence the result because of the assumption of a weightless rock mass). In addition, with io2 = 0°, i2 = arctan(ia) (see (12)). Table 4 shows the results obtained numerically (PhFDM) and analytically (PhS&O) (first chart proposed Figure 7) considering different values of ia. According to the margin error ratio observed in Table 4, less than 5%, it can be concluded that the two calculation methods have very similar results. In the numerical calculation, to estimate the bearing capacity, a stress path was formed until was reached, taking into account the whole wedge of the ground below the foundation. Therefore, the graphic output of the vertical component of the total stress tensor at failure was used to understand how the failure mechanism changed depending on ia. Figure 13 shows that the stress turned horizontally with an increase in ia, the rotation being larger when the horizontal force was wider. Additionally, it is noted that the vertical stress was well distributed in the case that ia = 0. Additionally, a convergence study was carried out, consisting of the analysis of the values of the bearing capacity obtained under the different increments of the velocity used. With a decrease in the value of the velocity increments, the result converged towards the final value that is the upper limit in the theoretical method (27b).
In the example studied, a seepage case was studied, where k v = 0, i o2 = 0 • , m o = 15, UCS = 100 MPa, GSI = 65, α = 0 • and B = 2.25 m were adopted (which did not influence the result because of the assumption of a weightless rock mass). In addition, with i o2 = 0 • , i 2 = arctan(i a ) (see (12)). Table 4 shows the results obtained numerically (P hFDM ) and analytically (P hS&O ) (first chart proposed Figure 7) considering different values of i a . According to the margin error ratio observed in Table 4, less than 5%, it can be concluded that the two calculation methods have very similar results. In the numerical calculation, to estimate the bearing capacity, a stress path was formed until was reached, taking into account the whole wedge of the ground below the foundation. Therefore, the graphic output of the vertical component of the total stress tensor at failure was used to understand how the failure mechanism changed depending on i a . Figure 13 shows that the stress turned horizontally with an increase in i a , the rotation being larger when the horizontal force was wider. Additionally, it is noted that the vertical stress was well distributed in the case that i a = 0. Sustainability 2020, 12, x FOR PEER REVIEW 20 of 22

Conclusions
An optimized design that adequately considers the negative effects that natural phenomena such as earthquakes and floods introduce to foundations is highly advisable to guarantee the sustainability of infrastructures.
Applying a pseudo-static approach that considers the seismic force and the seepage as an additional body force, a series of parameterized charts for estimating the bearing capacity of shallow foundations on rock masses were proposed. The charts were calculated through the analytical solution proposed by Serrano et al. [23] by means of a previous mathematical transformation of the angles of incidence in the boundaries of the pseudo-static loads produced by the inertial action.
Each chart was made according to the exponent "a" of the Modified Hoek and Brown criterion, the inclination α of the free boundary and kv of the foundation boundary (kv = 0 in the seepage case), and they were developed based on io2 and the horizontal component of the pseudo-static load of the foundation boundary (called kh in the seismic case and ia in the seepage case). It is noted that for high confining pressures in Boundary 1, it is not always possible to obtain a bearing capacity value; this limit is demarcated by the non-equilibrium line.
As expected, it was observed that the bearing capacity decreased with an increase in the pseudo-static load and the original inclination of the load. It was also observed that the value of the bearing capacity was not directly proportional to the bearing capacity factor Nß; therefore, a higher value of Nß is not associated with a greater value of Ph. In addition, the same value of Nß generates different values of Ph, depending on the other geotechnical parameters.
A validation using the finite difference method was carried out in a particular case. The numerical and analytical results, according to the example studied, show a variation of less than 5%. In the numerical graphic output of the vertical component of the total stress tensor, it is observed that the stress turned horizontally with an increase in ia. In addition, it is noted that the vertical stress was well distributed in the case that ia = 0.

Conclusions
An optimized design that adequately considers the negative effects that natural phenomena such as earthquakes and floods introduce to foundations is highly advisable to guarantee the sustainability of infrastructures.
Applying a pseudo-static approach that considers the seismic force and the seepage as an additional body force, a series of parameterized charts for estimating the bearing capacity of shallow foundations on rock masses were proposed. The charts were calculated through the analytical solution proposed by Serrano et al. [23] by means of a previous mathematical transformation of the angles of incidence in the boundaries of the pseudo-static loads produced by the inertial action.
Each chart was made according to the exponent "a" of the Modified Hoek and Brown criterion, the inclination α of the free boundary and k v of the foundation boundary (k v = 0 in the seepage case), and they were developed based on i o2 and the horizontal component of the pseudo-static load of the foundation boundary (called k h in the seismic case and i a in the seepage case). It is noted that for high confining pressures in Boundary 1, it is not always possible to obtain a bearing capacity value; this limit is demarcated by the non-equilibrium line.
As expected, it was observed that the bearing capacity decreased with an increase in the pseudo-static load and the original inclination of the load. It was also observed that the value of the bearing capacity was not directly proportional to the bearing capacity factor N β ; therefore, a higher value of N β is not associated with a greater value of P h . In addition, the same value of N β generates different values of P h , depending on the other geotechnical parameters.
A validation using the finite difference method was carried out in a particular case. The numerical and analytical results, according to the example studied, show a variation of less than 5%. In the numerical graphic output of the vertical component of the total stress tensor, it is observed that the