Three-Dimensional Upper Bound Solution to Estimate Soil Thrust of a Track System on Saturated Clay Slopes under Undrained Conditions

: This study proposes a three-dimensional upper bound solution for estimating the soil thrust of tracked vehicles on saturated clay slopes. The present study considered block, triangular wedge, and trapezoidal wedge failure modes to formulate an upper bound solution for each. The analytical solution for soil thrust was determined as the minimum upper bound solution among those for each failure mode. This analytical solution was validated through numerical simulations that modeled track-ground interactions. Parametric studies, based on the upper bound solution, assessed the impact of track system shape, vehicle weight, undrained shear strength, and ground slope on soil thrust. The analytical solutions and parametric studies provide a rapid method for assessing vehicle operability on clay slopes and offer references for designing tracked vehicles suitable for site conditions.


Introduction
Recently, there has been an increased focus on seabed resources [1], which has spurred the development of related mining technologies, such as the piping lifting system (e.g., [2,3]).In seabed mining, vehicles like backhoes, bulldozers, and excavators are essential for collecting and transporting resources, significantly impacting the efficiency of the process.Typically, these mining vehicles are heavy and equipped with a track system that offers better flotation and traction than wheels on the seabed, which typically consists of soft clay.Consequently, tracked vehicles are predominantly used in deep-sea mining as part of the piping lifting system, and their performance is extensively evaluated (e.g., [4][5][6][7]).
When a tracked vehicle moves, its weight and engine torque are transmitted to the ground via the track system as vertical and horizontal forces.These forces induce shear deformation in the ground, generating a shear force known as soil thrust [8].The vehicle advances as a reaction force overcomes this shear resistance; thus, soil thrust serves as a key performance indicator for the tracked vehicle, and its accurate estimation is important in seabed mining.
Bekker [9] initially introduced a soil thrust solution by assuming a rectangular block failure in the ground caused by the track system.Based on this, various researchers have conducted model experiments to demonstrate that block failure typically occurs in sandy soils and have developed an improved method for evaluating soil thrust (e.g., [10][11][12][13]).Additionally, further studies have both experimentally and numerically analyzed soil thrust considering the properties of the soil and the shape of the track system on sandy ground (e.g., [14,15]).
For a track system on clay ground, numerous experimental and numerical studies have focused on evaluating soil thrust based on the conditions of the tracked vehicle and the ground [16][17][18][19][20][21].Baek et al. [22] experimentally demonstrated that clay under a track system can experience a wedge failure mode as well as the block failure mode.Wang et al. [23] conducted a numerical analysis considering variables such as the strength parameters of clay and specifications of tracked vehicles, confirming the occurrence of the wedge failure mode.These findings suggest that under a track system, the failure mode in clay varies with the specifications of the track vehicle and the undrained shear strength of clay; thus, the estimation of the reasonable failure mode is one of the key factors in the evaluation of soil thrust of a track system over clay ground.
Woo and Baek [24] proposed upper bound solutions for soil thrust of a track system over clay ground, focusing on different failure modes (e.g., block, triangular wedge, and trapezoidal wedge) under the plane strain condition; subsequently, Woo and Beak [25] expanded on this research by developing a three-dimensional theoretical solution for soil thrust of an individual track system over clay ground; however, these studies only focused on tracked vehicles operating on flat terrain.To address this limitation, Na and Woo [26] introduced a theoretical solution for the soil thrust of a track system operating over a clay slope; however, this study only considered the soil thrust in the two-dimensional plane strain condition.
This study introduces a three-dimensional upper bound solution for soil thrust on a clay slope.The present study considered block, triangular wedge, and trapezoidal wedge failure modes, developing theoretical solutions for each.The seabed is typically composed of clayey deposits with undrained shear strengths ranging from 1.9 to 40 kPa [27]; thus, this study assumed the full saturation of clay.As the upper vehicle typically moves faster than the rate at which pore water drains in the clay, this study assumes undrained conditions within clay ground.According to Grečenko [28], the total thrust from tracked vehicles is the sum of the thrusts from individual track systems; therefore, this study concentrates on a theoretical solution for the soil thrust of a single-track system.This study also compared the theoretical solutions with numerical solutions for a track system on a clay slope.Finally, this paper presents parametric studies considering the shape of a track system, using dimensions from the commercial track system applied in Baek et al. [22], alongside the weight of the vehicle, slope angle, and the undrained shear strength of clay.

Assumption
Following Chen [29], in this study, the present upper bound solution assumes the following: (1) the undrained condition with fully saturated clay ground; (2) the kinematically admissible failure mechanism (or velocity field); (3) the occurrence of plastic dissipation limited to the failure surfaces derived from the admissible velocity field; (4) full dissipation of the external power to generate the plastic slip along the failure surface.Based on these assumptions, the present study considered three failure modes: block; triangular wedge; and trapezoidal wedge failure modes.In each failure mode, the soil block moves along the failure surfaces, and plastic dissipation is generated along the bottom and side failure planes.For a saturated clay under undrained conditions coupled with the Tresca plasticity model (the Mohr-Coulomb plasticity model with zero friction angle), the plastic dissipation rate along the failure surface D is given by Woo and Beak [25]: where c u is the undrained shear strength of the clay; V is the velocity along the failure surface, and A is the area of the failure surface.In the calculation of total external power, this study considered the vehicle weight and soil thrust as vertical and horizontal external forces, respectively.The weight of the soil block was ignored as the magnitude was very small compared to the vehicle weight.Based on these assumptions, the following sections derive the upper bound solutions for the soil thrust for each failure mode.

Block Failure Mode
Figure 1 schematically illustrates the block failure mechanism when a vehicle ascends a slope with angle β.In Figure 1, surface B represents the bottom failure plane, while surfaces S 1 and S 2 are the lateral failure planes.The dimensions H, L, and D denote the height, length, and width of a single-track system, respectively.The force W g , aligned with the gravitational force, represents the weight transmitted from the upper vehicle.The force F b , parallel to the slope, indicates the soil thrust.In the block failure mode, the brick-shaped soil block moves at velocity V.
Appl.Sci.2024, 14, 4335 3 of 14 small compared to the vehicle weight.Based on these assumptions, the following sections derive the upper bound solutions for the soil thrust for each failure mode.

Block Failure Mode
Figure 1 schematically illustrates the block failure mechanism when a vehicle ascends a slope with angle β.In Figure 1, surface B represents the bottom failure plane, while surfaces S1 and S2 are the lateral failure planes.The dimensions H, L, and D denote the height, length, and width of a single-track system, respectively.The force Wg, aligned with the gravitational force, represents the weight transmitted from the upper vehicle.The force Fb, parallel to the slope, indicates the soil thrust.In the block failure mode, the brickshaped soil block moves at velocity V.In these conditions, the total external power P is given by which corresponds to the sum of inner products between external forces (Wg and Fb) and velocity V.As there exist one bottom and two side failure planes, the total dissipation rate D is expressed as where Db and Ds represent plastic dissipation rates along the bottom failure surface B and side failure surfaces (S1 and S2), respectively.Along surface B, plastic dissipation rate Db is written as Along two side failure surfaces, plastic dissipation rate Ds is Thus, the total dissipation rate D is expressed as By assuming that all energy is dissipated to generate plastic deformation, the total external power P is identical to the total (internal) plastic dissipation rate D; thus, equating P = D gives the soil thrust Fb as which shows that the soil thrust from the block failure mode is deterministic.In these conditions, the total external power P is given by which corresponds to the sum of inner products between external forces (W g and F b ) and velocity V.As there exist one bottom and two side failure planes, the total dissipation rate D is expressed as where D b and D s represent plastic dissipation rates along the bottom failure surface B and side failure surfaces (S 1 and S 2 ), respectively.Along surface B, plastic dissipation rate D b is written as Along two side failure surfaces, plastic dissipation rate D s is Thus, the total dissipation rate D is expressed as By assuming that all energy is dissipated to generate plastic deformation, the total external power P is identical to the total (internal) plastic dissipation rate D; thus, equating P = D gives the soil thrust F b as which shows that the soil thrust from the block failure mode is deterministic.

Triangular Wedge Failure Mode
Figure 2 schematically describes the failure mechanism of the triangular wedge failure mode.In Figure 2, surface B is an inclined bottom failure plane forming an angle θ with the ground surface, while surfaces S 1 and S 2 are two triangular side failure planes.In the triangular wedge failure mode, the prism-shaped soil block moves at velocity V.

Triangular Wedge Failure Mode
Figure 2 schematically describes the failure mechanism of the triangular wedge failure mode.In Figure 2, surface B is an inclined bottom failure plane forming an angle θ with the ground surface, while surfaces S1 and S2 are two triangular side failure planes.In the triangular wedge failure mode, the prism-shaped soil block moves at velocity V.In this case, total external power P is expressed as where Ft represents the soil thrust under the triangular wedge failure mode.For the triangular wedge failure mode, plastic dissipation rate Db along the bottom failure surface B is defined in whereas, along the two triangular-shaped side failure surfaces S1 and S2, plastic dissipation rate Ds is expressed as Consequently, the total dissipation rate D is formulated as Assuming that total external power is entirely dissipated to produce the plastic deformation along failure surfaces, the soil thrust Ft can be written as which shows that the soil thrust Ft is a function of the geometric variable θ.Therefore, the stationary upper bound solution is identified at the angle θ where ∂Ft/∂θ = 0; thus, which can be reorganized by which is a third-order polynomial equation of sinθ.Although there is an analytical way to solve Equation ( 14) for θ, the present research numerically solves it using the bisection method, which is very efficient when the solution range is definite (in this problem, 0 < θ < tan −1 (L/H)).In this case, total external power P is expressed as where F t represents the soil thrust under the triangular wedge failure mode.For the triangular wedge failure mode, plastic dissipation rate D b along the bottom failure surface B is defined in whereas, along the two triangular-shaped side failure surfaces S 1 and S 2 , plastic dissipation rate D s is expressed as Consequently, the total dissipation rate D is formulated as Assuming that total external power is entirely dissipated to produce the plastic deformation along failure surfaces, the soil thrust F t can be written as which shows that the soil thrust F t is a function of the geometric variable θ.Therefore, the stationary upper bound solution is identified at the angle θ where ∂F t /∂θ = 0; thus, which can be reorganized by Appl.Sci.2024, 14, 4335 which is a third-order polynomial equation of sinθ.Although there is an analytical way to solve Equation ( 14) for θ, the present research numerically solves it using the bisection method, which is very efficient when the solution range is definite (in this problem, 0 < θ < tan −1 (L/H)).

Trapezoidal Wedge Failure Mode
Figure 3 schematically illustrates the failure mechanism in the trapezoidal wedge failure mode.In this mode, a trapezoid-shaped soil block forming an angle θ with the ground surface moves at velocity V.

Trapezoidal Wedge Failure Mode
Figure 3 schematically illustrates the failure mechanism in the trapezoidal wedge failure mode.In this mode, a trapezoid-shaped soil block forming an angle θ with the ground surface moves at velocity V.Under the trapezoidal wedge failure mode, the total P power is expressed as The dissipation rates Db (along the slip surface B) and Ds (along the two side failure surfaces S1 and S2) are given as The total dissipation rate D is then Assuming that the total work is fully dissipated, the following equation is derived.
Thus, the soil thrust Fp is Given that cu, H, L, D, and Wg are known in Equation ( 20), the least upper solution obtained by setting the gradient of Fp with respect to θ should be zero: which yields Equation ( 22) represents a nonlinear equation of sinθ.This study employed the bisection method to solve it, with θ ranging from tan −1 (L/H) to 90°.Once θ is determined, the soil thrust Fp can be calculated using Equation (20).Under the trapezoidal wedge failure mode, the total P power is expressed as The dissipation rates D b (along the slip surface B) and D s (along the two side failure surfaces S 1 and S 2 ) are given as The total dissipation rate D is then Assuming that the total work is fully dissipated, the following equation is derived.
Thus, the soil thrust F p is Given that c u , H, L, D, and W g are known in Equation ( 20), the least upper solution obtained by setting the gradient of F p with respect to θ should be zero: Appl.Sci.2024, 14, 4335 6 of 13 which yields Equation ( 22) represents a nonlinear equation of sinθ.This study employed the bisection method to solve it, with θ ranging from tan −1 (L/H) to 90 • .Once θ is determined, the soil thrust F p can be calculated using Equation (20).

Upper Bound Solution
The upper bound theorem states that the collapse loads derived from any admissible velocity field are always higher than the actual collapse load [29]; therefore, the lowest load among those calculated from the various admissible velocity fields serves as a reasonable approximation of the actual collapse load.In this study, the soil thrust F x is determined as the minimum value among F b from the block failure mode, F t from the triangular wedge failure mode, and F p from the trapezoidal failure mode, as shown in Equation ( 23) and Figure 4.

Upper Bound Solution
The upper bound theorem states that the collapse loads derived from any admissible velocity field are always higher than the actual collapse load [29]; therefore, the lowest load among those calculated from the various admissible velocity fields serves as a reasonable approximation of the actual collapse load.In this study, the soil thrust Fx is determined as the minimum value among Fb from the block failure mode, Ft from the triangular wedge failure mode, and Fp from the trapezoidal failure mode, as shown in Equation ( 23) and Figure 4.

Numericla Modeling
Woo and Beak [25] verified the analytical solution of soil thrust through numerical analysis and model experiments on a flat ground surface in their study; following this, this study performed numerical analyses to validate the proposed theoretical solution for a track system on a clay slope.The interaction between the track and the ground was simulated using ABAQUS Explicit.Table 1 lists the material properties of both the ground and the track system considered in the numerical analyses.The track system was modeled as a linear elastic material employing the material properties of steel.In the numerical analysis, this study described the mechanical responses of the saturated clay using a linear elastic-perfect plastic constitutive model coupled with the Tresca yield criterion (Mohr-Coulomb yield criterion with zero friction angle).The elastic behavior was modeled as a linear elastic model with Young's modulus E = 20,000 kPa and Poisson's ratio ν = 0.49 (to describe the undrained responses), while its plastic behavior was described using the Tresca model with its yield strength corresponding to the undrained shear strength of clay assuming the perfect plastic condition.To examine the influence of the undrained shear strength c u of clay, this study conducted the numerical analyses by varying its values to 5, 7.5, 10, 15, 20, and 25 kPa.Figure 5 depicts the geometric modeling of a track system and subgrade clay ground in the numerical analyses of this study.To improve the efficiency of the numerical analysis, a half section of the entire domain considering x-symmetricity was modeled, as shown in Figure 5.In Figure 5, the ground slope angle β was set at 10 • ; regarding the track system, its height, length, and width were set to 0.035, 0.124, and 0.062 m, respectively.To speed up the numerical analyses, a mass scaling factor of 400 was applied to all simulations.The outer boundaries of the ground (except the symmetric boundary) were modeled using infinite elements to minimize the disturbance caused by reflection.
this study performed numerical analyses to validate the proposed theoretical solution for a track system on a clay slope.The interaction between the track and the ground was simulated using ABAQUS Explicit.Table 1 lists the material properties of both the ground and the track system considered in the numerical analyses.The track system was modeled as a linear elastic material employing the material properties of steel.In the numerical analysis, this study described the mechanical responses of the saturated clay using a linear elastic-perfect plastic constitutive model coupled with the Tresca yield criterion (Mohr-Coulomb yield criterion with zero friction angle).The elastic behavior was modeled as a linear elastic model with Young's modulus E = 20,000 kPa and Poisson's ratio ν = 0.49 (to describe the undrained responses), while its plastic behavior was described using the Tresca model with its yield strength corresponding to the undrained shear strength of clay assuming the perfect plastic condition.To examine the influence of the undrained shear strength cu of clay, this study conducted the numerical analyses by varying its values to 5, 7.5, 10, 15, 20, and 25 kPa.5 depicts the geometric modeling of a track system and subgrade clay ground in the numerical analyses of this study.To improve the efficiency of the numerical analysis, a half section of the entire domain considering x-symmetricity was modeled, as shown in Figure 5.In Figure 5, the ground slope angle β was set at 10°; regarding the track system, its height, length, and width were set to 0.035, 0.124, and 0.062 m, respectively.To speed up the numerical analyses, a mass scaling factor of 400 was applied to all simulations.The outer boundaries of the ground (except the symmetric boundary) were modeled using infinite elements to minimize the disturbance caused by reflection.The numerical analysis consists of two stages: loading; and pushing.In the loading stage, the surface traction (σ c = 11.2 kPa = W g /(LD)) was uniformly applied as traction on the upper surface of the track.The perpendicular and tangential components of traction with respect to the surface were σ c sinβ and σ c cosβ, respectively.The interface between the track and clay was modeled as the frictionless contact model.In the pushing stage, the front panel of the track was pushed at 1 mm/s velocity parallel to the ground surface (that had an angle β with the horizontal surface); this stage was maintained for 7.5 s, and the final displacement of the track was 7.5 mm.
To verify the convergence of the numerical simulations, this study carried out multiple simulations under consistent conditions while adjusting the mesh density.Figure 6 illustrates three different mesh configurations considered in this study.The section around the track system has the smallest and most uniform mesh; the mesh spacing increases gradually toward the outer boundaries.In Figure 6a, which has the lowest mesh density, the smallest node spacing is 3.6 mm; the total numbers of nodes and elements are 43,865 and 39,530, respectively.In Figure 6b, the intermediate mesh density case, the smallest node spacing is 2.4 mm; the total numbers of nodes and elements are 96,131 and 88,684, respectively.In Figure 6c, the highest mesh density case, the smallest node spacing is 1.8 mm; the total numbers of nodes and elements are 162,893 and 152,324, respectively.front panel of the track was pushed at 1 mm/s velocity parallel to the ground surface (that had an angle β with the horizontal surface); this stage was maintained for 7.5 s, and the final displacement of the track was 7.5 mm.
To verify the convergence of the numerical simulations, this study carried out multiple simulations under consistent conditions while adjusting the mesh density.Figure 6 illustrates three different mesh configurations considered in this study.The section around the track system has the smallest and most uniform mesh; the mesh spacing increases gradually toward the outer boundaries.In Figure 6a, which has the lowest mesh density, the smallest node spacing is 3.6 mm; the total numbers of nodes and elements are 43,865 and 39,530, respectively.In Figure 6b, the intermediate mesh density case, the smallest node spacing is 2.4 mm; the total numbers of nodes and elements are 96,131 and 88,684, respectively.In Figure 6c, the highest mesh density case, the smallest node spacing is 1.8 mm; the total numbers of nodes and elements are 162,893 and 152,324, respectively.

Comparison of Upper Bound Solutions with Numerical Analysis Result
Figure 7 shows the soil thrust (both upper bound and numerical solutions) with respect to cu with the identical σc, β, and dimensions H, L, and D, as cited by Beak et al. [22].In the numerical analysis results, the soil thrust is calculated by doubling the sum of the reaction force at all nodes on the front panel of the track, where velocity was constrained.Once the soil thrust converged in the numerical analyses, the final soil thrust was estimated by averaging soil thrust values for 1.5 s.
In Figure 7, solid, dashed, and dotted lines represent Fb (from the block failure mode), Ft (from the triangle wedge failure mode), and Fp (from the trapezoidal failure mode), while the thick gray solid line represents the upper bound solution Fx, which is the minimum between Fb, Ft, and Fp.In Figure 7, when cu is less than about 8 kPa, Fb is less than Ft and Fp; thus, in this case, Fx = Fb.If cu is greater than 8 kPa, Ft is smaller than Fb and Fp; thus, in this case, Fx = Ft.Figure 7 also plots the numerical solutions overlapped with the upper bound solutions.In Figure 7, the green, blue, and red symbols correspond to the numerical solutions with the lowest, intermediate, and highest mesh densities.As mesh density increases, the numerical solution more closely approximates the analytical solution proposed in this study; the general rule in the numerical simulation is that the numerical analysis with the finer mesh approaches the exact solution; thus, the analytical solution in this study is very close to the exact solution.

Comparison of Upper Bound Solutions with Numerical Analysis Result
Figure 7 shows the soil thrust (both upper bound and numerical solutions) with respect to c u with the identical σ c , β, and dimensions H, L, and D, as cited by Beak et al. [22].In the numerical analysis results, the soil thrust is calculated by doubling the sum of the reaction force at all nodes on the front panel of the track, where velocity was constrained.Once the soil thrust converged in the numerical analyses, the final soil thrust was estimated by averaging soil thrust values for 1.5 s. Figure 8 shows the distribution of equivalent plastic strain in the ground from the results of the numerical simulation with the highest mesh density.When cu = 4 kPa (Figure 8a), the numerical solution clearly shows the block failure, which is identical to the failure mode determined by the upper bound solution.When cu =15, 20, and 25 kPa (Figure 8df), the numerical solutions definitely depict the triangular wedge failure mode, which is identical to the failure mode determined by the upper bound solution; in Figure 8d-f, as cu increases, the length of the triangular wedge decreases.When cu is 7.5 and 10 kPa, as shown in Figure 8b,c, as the failure mode transitions from the block to the triangular In Figure 7, solid, dashed, and dotted lines represent F b (from the block failure mode), F t (from the triangle wedge failure mode), and F p (from the trapezoidal failure mode), while the thick gray solid line represents the upper bound solution F x , which is the minimum Appl.Sci.2024, 14, 4335 9 of 13 between F b , F t , and F p .In Figure 7, when c u is less than about 8 kPa, F b is less than F t and F p ; thus, in this case, F x = F b .If c u is greater than 8 kPa, F t is smaller than F b and F p ; thus, in this case, F x = F t .Figure 7 also plots the numerical solutions overlapped with the upper bound solutions.In Figure 7, the green, blue, and red symbols correspond to the numerical solutions with the lowest, intermediate, and highest mesh densities.As mesh density increases, the numerical solution more closely approximates the analytical solution proposed in this study; the general rule in the numerical simulation is that the numerical analysis with the finer mesh approaches the exact solution; thus, the analytical solution in this study is very close to the exact solution.
Figure 8 shows the distribution of equivalent plastic strain in the ground from the results of the numerical simulation with the highest mesh density.When c u = 4 kPa (Figure 8a), the numerical solution clearly shows the block failure, which is identical to the failure mode determined by the upper bound solution.When c u =15, 20, and 25 kPa (Figure 8d-f), the numerical solutions definitely depict the triangular wedge failure mode, which is identical to the failure mode determined by the upper bound solution; in Figure 8d-f, as c u increases, the length of the triangular wedge decreases.When c u is 7.5 and 10 kPa, as shown in Figure 8b,c, as the failure mode transitions from the block to the triangular wedge modes, numerical analyses simulated both failure modes concurrently.This implies that, as c u increases, rather than an abrupt change in the failure mode, both failure modes occur simultaneously.As a result, following the comparison with the numerical solution, the proposed upper bound solution can accurately predict the soil thrust value and the corresponding failure mode.Figure 8 shows the distribution of equivalent plastic strain in the ground from the results of the numerical simulation with the highest mesh density.When cu = 4 kPa (Figure 8a), the numerical solution clearly shows the block failure, which is identical to the failure mode determined by the upper bound solution.When cu =15, 20, and 25 kPa (Figure 8df), the numerical solutions definitely depict the triangular wedge failure mode, which is identical to the failure mode determined by the upper bound solution; in Figure 8d-f, as cu increases, the length of the triangular wedge decreases.When cu is 7.5 and 10 kPa, as shown in Figure 8b,c, as the failure mode transitions from the block to the triangular wedge modes, numerical analyses simulated both failure modes concurrently.This implies that, as cu increases, rather than an abrupt change in the failure mode, both failure modes occur simultaneously.As a result, following the comparison with the numerical solution, the proposed upper bound solution can accurately predict the soil thrust value and the corresponding failure mode.

Parametric Study
This study analyzed the influence of parameters (geometry, vehicle design, and ground strength) on the soil thrust over clay slope using the proposed upper bound theoretical solution.Table 2 lists the parameters considered; this study considered the ground slope as a geometry parameter and undrained shear strength cu as a ground strength parameter.For the vehicle parameters, this study included ground contact pressure σc, track height H, track length L, and track width D. Table 2 also presents the values and ranges of the considered parameters.
According to Costello et al. [30], 95.5% of the seabed has slopes of 6° or less, while only 4.5% exceed this angle; consequently, the slope angle (β) was set within a range of 0° (which represented the flat surface) to 10° (which corresponded to the extreme slope).For the ground contact pressure (σc), this study established a range of 5 to 50 kPa to cover the diverse types of underwater heavy vehicles.For the undrained shear strength cu, following Spagnoli, Freudenthal, Spagnoli et al. [27], this study set a range of cu from 2 to 40 kPa.Following Baek et al. [22], the height, length, and width of a track system were set to 0.035, 0.124, and 0.124 m, respectively.Additionally, to assess the impact of these aspect ratios,

Parametric Study
This study analyzed the influence of parameters (geometry, vehicle design, and ground strength) on the soil thrust over clay slope using the proposed upper bound theoretical solution.Table 2 lists the parameters considered; this study considered the ground slope as a geometry parameter and undrained shear strength c u as a ground strength parameter.For the vehicle parameters, this study included ground contact pressure σ c , track height H, track length L, and track width D.  According to Costello et al. [30], 95.5% of the seabed has slopes of 6 • or less, while only 4.5% exceed this angle; consequently, the slope angle (β) was set within a range of 0 • (which represented the flat surface) to 10 • (which corresponded to the extreme slope).For the ground contact pressure (σ c ), this study established a range of 5 to 50 kPa to cover the diverse types of underwater heavy vehicles.For the undrained shear strength c u , following Spagnoli, Freudenthal, Spagnoli et al. [27], this study set a range of c u from 2 to 40 kPa.Following Baek et al. [22], the height, length, and width of a track system were set to 0.035, 0.124, and 0.124 m, respectively.Additionally, to assess the impact of these aspect ratios, the upper bound solution was calculated and compared while varying the dimensions of height, length, and width.
This study parametrically studied the responses of the stationary upper bound solution F x (=min(F b , F t , F p )) with respect to the considered parameters (β, c u , σ c , H, L, and D).
Figure 9 illustrates how soil thrust responds to the changes in β and σ c , while c u and H, L, and D remain constant (dimensions cited from Beak et al. [22]).Figure 9a plots the soil thrust versus the angle β, where the block failure predominates due to low undrained shear strength (c u = 2 kPa), while Figure 9b illustrates the soil thrust against the angle β, where the triangle wedge failure predominates because of high undrained shear strength (c u = 40 kPa).As the angle β increases for all failure modes, there is a noticeable decrease in soil thrust; this decrease is further amplified as σ c increases.In very soft ground conditions (Figure 9a), soil thrust becomes negative with higher β and σ c , indicating that the tracked vehicle cannot climb the slope and begins to slip downward.On the other hand, in stiff ground conditions (Figure 9b), where the triangular wedge failure is dominant, the soil thrust remains positive, enabling the vehicle to climb the slope.Figure 10 displays how soil thrust changes with cu and σc while β is constant, and dimensions H, L, D from Beak et al. [22] are also fixed.As σc increases, the kink point, where the failure mode transitions from the block failure mode to triangle wedge failure mode, moves to the upper left side.This implies that the tendency of the block failure mode increases with the vehicle load.In Figure 10 (and Figure 9a), when σc increases with low cu, negative soil thrust occurs, causing the vehicle to slip down.Before the kink points (where block failure mode happens), with fixed cu, higher soil thrust is generated as σc is lower (cu = 2 kPa in Figure 10); however, after the kink points (where triangular wedge failure mode happens), with fixed cu, the higher soil thrust is generated as σc is higher (cu = 35 kPa in Figure 10).Figure 10 displays how soil thrust changes with c u and σ c while β is constant, and dimensions H, L, D from Beak et al. [22] are also fixed.As σ c increases, the kink point, where the failure mode transitions from the block failure mode to triangle wedge failure mode, moves to the upper left side.This implies that the tendency of the block failure mode increases with the vehicle load.In Figure 10 (and Figure 9a), when σ c increases with low c u , negative soil thrust occurs, causing the vehicle to slip down.Before the kink points (where block failure mode happens), with fixed c u , higher soil thrust is generated as σ c is lower (c u = 2 kPa in Figure 10); however, after the kink points (where triangular wedge failure mode happens), with fixed c u , the higher soil thrust is generated as σ c is higher (c u = 35 kPa in Figure 10).Figure 10 displays how soil thrust changes with cu and σc while β is constant, and dimensions H, L, D from Beak et al. [22] are also fixed.As σc increases, the kink point, where the failure mode transitions from the block failure mode to triangle wedge failure mode, moves to the upper left side.This implies that the tendency of the block failure mode increases with the vehicle load.In Figure 10 (and Figure 9a), when σc increases with low cu, negative soil thrust occurs, causing the vehicle to slip down.Before the kink points (where block failure mode happens), with fixed cu, higher soil thrust is generated as σc is lower (cu = 2 kPa in Figure 10); however, after the kink points (where triangular wedge failure mode happens), with fixed cu, the higher soil thrust is generated as σc is higher (cu = 35 kPa in Figure 10).Subsequently, this study examined the effect of dimensions (H, L, and D) of the track system on vehicle slippage on a clay slope.As shown in Figure 9, the most critical case is that the heavy vehicle climbs up a slope consisting of very soft clay, where the block failure mode is dominant; thus, the present analysis only considered the block failure mode Subsequently, this study examined the effect of dimensions (H, L, and D) of the track system on vehicle slippage on a clay slope.As shown in Figure 9, the most critical case is that the heavy vehicle climbs up a slope consisting of very soft clay, where the block failure mode is dominant; thus, the present analysis only considered the block failure mode with fixed slope (β = 8 • ) and σ c (=40 kPa) and with varied undrained shear strength c u and dimensions (H, L, and D) of the track.
Figure 11 illustrates soil thrust versus c u with the different dimensions of the track system.The dashed line represents soil thrust for the default dimensions (H = 0.035 m, L = 0.124 m, and D = 0.124 m, cited from Beak et al. [22]), while the solid colored lines indicate soil thrust against c u for modified dimensions ((3H, L, D), (H, 3L, D), and (H, L, 3D)).As illustrated in Figure 11, increasing H (solid blue line) leads to a rise in both the intercept and slope of the soil thrust curve with respect to c u ; however, increasing L or D (solid green or orange line) increases the slope of the soil thrust curve while reducing its intercept, which does not cut negative soil thrust region along the horizontal axis (or c u axis).With the definition of σ c (=W g /(LD)), Equation ( 7) can be rewritten as In Equation ( 24), an increase in H causes the first term c u L(D + 2H) to increase, while the second term σ c LDsinβ remains unchanged; accordingly, the soil thrust F b increases.As L and D increase, however, the first and second terms increase together, and the second term σ c LDsinβ is also proportional to σ c .Since soil thrust is calculated by subtracting these terms, increases in L and D do not significantly boost soil thrust; in fact, they sometimes result in a reduction as the second term σ c LDsinβ increases.Therefore, when operating underwater heavy tracked vehicles on very soft ground where negative soil thrust is a concern, it is necessary to select a more effective strategy adjusted to site conditions: an increase in the track height or a reduction in the vehicle weight.
the second term σcLDsinβ remains unchanged; accordingly, the soil thrust Fb increases.As L and D increase, however, the first and second terms increase together, and the second term σcLDsinβ is also proportional to σc.Since soil thrust is calculated by subtracting these terms, increases in L and D do not significantly boost soil thrust; in fact, they sometimes result in a reduction as the second term σcLDsinβ increases.Therefore, when operating underwater heavy tracked vehicles on very soft ground where negative soil thrust is a concern, it is necessary to select a more effective strategy adjusted to site conditions: an increase in the track height or a reduction in the vehicle weight.

Conclusions
This study proposes a three-dimensional upper bound solution for estimating soil thrust on a clay slope for tracked vehicles.It assumes that the failure modes of clayey soil include block, triangular wedge, and trapezoidal wedge failures.An upper bound solution for each mode was formulated and provided.Using the principle of the least action, the minimum upper bound solution was selected as the analytical soil thrust.This analytical solution was validated through numerical analyses that simulated track-ground interactions.
The parametric analysis of the analytical solution on clay slope shows that only block and triangular wedge failures occur, with trapezoidal wedge failure absent.As the ground slope angle increases, there is a corresponding decrease in soil thrust.For a constant slope angle, the reduction in soil thrust became more noticeable as the vehicle weight increased;

Conclusions
This study proposes a three-dimensional upper bound solution for estimating soil thrust on a clay slope for tracked vehicles.It assumes that the failure modes of clayey soil include block, triangular wedge, and trapezoidal wedge failures.An upper bound solution for each mode was formulated and provided.Using the principle of the least action, the minimum upper bound solution was selected as the analytical soil thrust.This analytical solution was validated through numerical analyses that simulated trackground interactions.
The parametric analysis of the analytical solution on clay slope shows that only block and triangular wedge failures occur, with trapezoidal wedge failure absent.As the ground slope angle increases, there is a corresponding decrease in soil thrust.For a constant slope angle, the reduction in soil thrust became more noticeable as the vehicle weight increased; in extreme cases, this can result in negative soil thrust, causing the vehicle to slip and be unable to move forward.
Although further validation through model and field experiments with consideration of various clay ground and tracked vehicles is desirable to confirm the applicability of the proposed solution, the analytical solution and parametric analysis results presented in this study offer a rapid way to assess vehicle operability on a clay slope in real time.

Figure 1 .
Figure 1.Schematic diagram of the block failure mode for a single-track system.

Figure 1 .
Figure 1.Schematic diagram of the block failure mode for a single-track system.

Figure 2 .
Figure 2. Schematic diagram of the triangular wedge failure mode for a single-track system.

Figure 2 .
Figure 2. Schematic diagram of the triangular wedge failure mode for a single-track system.

Figure 3 .
Figure 3. Schematic diagram of the trapezoidal wedge failure mode for a single-track system.

Figure 3 .
Figure 3. Schematic diagram of the trapezoidal wedge failure mode for a single-track system.

Figure 4 .
Figure 4. Soil thrust Fx versus the undrained shear strength cu of clay with fixed dimensions of a track system (height H, length L, and width D), as referred to in Beak et al. [22]; ground slope β and weight Wg transmitted from the upper vehicle.

Figure 4 .
Figure 4. Soil thrust F x versus the undrained shear strength c u of clay with fixed dimensions of a track system (height H, length L, and width D), as referred to in Beak et al. [22]; ground slope β and weight W g transmitted from the upper vehicle.

Figure 5 .
Figure 5. Geometrical modeling of a track system over clay slope for the numerical simulations in this study.

Figure 5 .
Figure 5. Geometrical modeling of a track system over clay slope for the numerical simulations in this study.

Figure 6 .
Figure 6.Mesh configurations applied to the numerical analysis: (a) lowest mesh density; (b) intermediate mesh density; and (c) highest mesh density.

Figure 6 .
Figure 6.Mesh configurations applied to the numerical analysis: (a) lowest mesh density; (b) intermediate mesh density; and (c) highest mesh density.

14 Figure 7 .
Figure 7. Upper bound solutions (Fb, Ft, Fp) and minimum upper bound solution Fx overlapped with numerical solutions with different mesh densities.

Figure 7 .
Figure 7. Upper bound solutions (F b , F t , F p ) and minimum upper bound solution F x overlapped with numerical solutions with different mesh densities.

Figure 7 .
Figure 7. Upper bound solutions (Fb, Ft, Fp) and minimum upper bound solution Fx overlapped with numerical solutions with different mesh densities.

Figure 9 .
Figure 9. Soil thrust versus ground slope angle with varying ground contact pressure σc with (a) very soft clay ground (cu = 2 kPa) and (b) stiff clay ground (cu = 40 kPa).

Figure 9 .
Figure 9. Soil thrust versus ground slope angle with varying ground contact pressure σ c with (a) very soft clay ground (c u = 2 kPa) and (b) stiff clay ground (c u = 40 kPa).

Figure 9 .
Figure 9. Soil thrust versus ground slope angle with varying ground contact pressure σc with (a) very soft clay ground (cu = 2 kPa) and (b) stiff clay ground (cu = 40 kPa).

Figure 10 .
Figure 10.Soil thrust versus undrained shear strength cu with different ground contact pressure σc.

Figure 10 .
Figure 10.Soil thrust versus undrained shear strength c u with different ground contact pressure σ c .

Figure 11 .
Figure 11.Examination of the influence of dimensions of a track system.

Figure 11 .
Figure 11.Examination of the influence of dimensions of a track system.

Table 1 .
Material parameters in the numerical analyses in this study.

Table 1 .
Material parameters in the numerical analyses in this study.

Table 2
also presents the values and ranges of the considered parameters.

Table 2 .
Parameters and their ranges in the parametric study.