Xue Zhang
email: xue.zhang@uon.edu.au

Daichao Sheng

Scott W Sloan

Jeremy Bleyer

Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity

Keywords: Sensitive clays, Progressive failure, Elastoviscoplasticity, Strain softening, PFEM, Mathematical programming, SOCP

This paper presents a Lagrangian formulation of elastoviscoplasticity, based on the Particle Finite Element Method, for progressive failure analysis of sensitive clays. The sensitive clay is represented by an elastoviscoplastic model which is a mixture of the Bingham model, for describing rheological behaviour, and the Tresca model with strain softening for capturing the progressive failure behaviour. The finite element formulation for the incremental elastoviscoplastic analysis is reformulated, through the application of the Hellinger-Reissner variational theorem, as an equivalent optimization program that can be solved efficiently using modern algorithms such as the interior-point method. The recast formulation is then incorporated into the framework of the Particle Finite Element Method for investigating progressive failure problems related to sensitive clays, such as the collapse of a sensitive clay column and the retrogressive failure of a slope in sensitive clays, where extremely large deformation occurs.

Introduction

Sensitive clay is characterized by a decrease in its shear strength when experiencing plastic deformation. A highly sensitive clay may possess sensitivity, defined as a ratio of the undisturbed shear strength and the remoulded shear strength, of the order of magnitude of a hundred. For example, the reported values of the sensitivity of the clay involved in the 1893 Verdalen landslide and the 2012 Byneset landslide are 300 and 120, respectively [START_REF] Thakur | Quickness of sensitive clays[END_REF]. Due to the strong strain-softening behaviour, geostructures built on a layer of sensitive clay often fail in a progressive manner. Moreover, unexpectedly catastrophic failure of the geostructure might also be induced by a small perturbation. Typical examples are the multiple retrogressive slides and spreads in sensitive clays observed in Canada and Scandinavia [START_REF] Locat | Progressive failures in eastern Canadian and Scandinavian sensitive clays[END_REF], which occurred suddenly, covered large areas (more than 1 hectare) and were caused by an initially small slope failure.

Reliable prediction of the progressive failure behaviour of sensitive clays is of critical importance. It can provide guidelines for relevant engineering practice, for example construction on sensitive clays, and also assist in minimizing the degree of destruction caused by potential geohazards (such as the fore-mentioned large landslides). Although numerical simulation is a powerful tool for analyzing complex geotechnical problems, robust modelling of the large deformations induced by progressive failure in sensitive clays is still a formidable task. Indeed, a major challenge is the complex behaviour that is typically exhibited by sensitive clays. An undisturbed sensitive clay usually behaves like a solid body, but may change to be a semi-liquid material after being remolded [START_REF] Thakur | Quickness of sensitive clays[END_REF]. The transformation between these two states is caused by strain softening. Suitable constitutive models must be capable of describing the rheological behaviour of a sensitive clay, since this is crucial for estimating the run-out distance of landslides [START_REF] Jeong | Determining the viscosity and yield surface of marine sediments using modified Bingham models[END_REF][START_REF] De Blasio | Flow models of natural debris flows originating from overconsolidated clay materials[END_REF][START_REF] Locat | Viscosity, yield stress, remolded strength, and liquidity index relationships for sensitive clays[END_REF] as well as capturing the strain-softening behaviour that contributes to the phenomenon of progressive failure [START_REF] Locat | Progressive failures in eastern Canadian and Scandinavian sensitive clays[END_REF][START_REF] Quinn | Development of progressive failure in sensitive clay slopes[END_REF][START_REF] Bernander | Downhill progressive landslides in long natural slopes: triggering agents and landslide phases modeled with a finite difference method[END_REF]. Sensitive clays typically undergo extremely large deformation along localized shear zones due to strain-softening. This feature can cause severe mesh distortion when the traditional finite element method is adopted and result in computational difficulties. Additionally, the free-surface evolution induced by extreme deformation also challenges the use of the traditional FEM because of its use of a fixed mesh topology. Recently, some alternative numerical approaches have been proposed for modelling the progressive failure of sensitive clays involving large deformation.

Wang et. al [START_REF] Wang | Investigation of retrogressive and progressive slope failure mechanisms using the material point method[END_REF] studied retrogressive and progressive slope failure in sensitive clays using the material point method. Dey et. al [START_REF] Dey | Large deformation finite-element modelling of progressive failure leading to spread in sensitive clay slopes[END_REF][START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF][START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF] analyzed the spread in sensitive clay slopes due to progressive failure by implementing a strain-softening model into the ABAQUS Coupled Eulerian Lagrangian approach. Although these procedures reproduced the pronounced progressive failure behaviour of sensitive clays, it is notable that classical rate-independent models were utilised. However, ignoring the rheology of sensitive clays may lead to the inaccurate predictions. Analytical approaches, such as shear band propagation approaches [START_REF] Germanovich | Dynamic growth of slip surfaces in catastrophic landslides[END_REF][START_REF] Puzrin | Simple criteria for ploughing and runout in post-failure evolution of submarine landslides[END_REF][START_REF] Viesca | Nucleation of slip-weakening rupture instability in landslides by localized increase of pore pressure[END_REF][START_REF] Palmer | The Growth of Slip Surfaces in the Progressive Failure of Over-Consolidated Clay[END_REF][START_REF] Puzrin | The growth of shear bands in the catastrophic failure of soils[END_REF], have also been used to study the progressive failure process in catastrophic landslides in nature. Recent developments in the shear band propagation approach for analyzing catastrophic and progressive failure are summarized in [START_REF] Puzrin | Shear band propagation analysis of submarine slope stability[END_REF].

This paper provides an alternative Lagrangian computational approach for the analysis of progressive failure of sensitive clays involving extremely large deformation. An advanced elastoviscoplastic constitutive relationship, which is a combination of the Bingham model and the Tresca model with strain softening, is adopted for describing their complex behaviour.

To solve the resulting elastoviscoplastic problem with strain softening, a generalized incremental Hellinger-Reissner variational theorem [START_REF] Reissner | On a variational theorem in elasticity[END_REF] is proposed which recasts the associated governing equations into an equivalent min-max program. After finite element discretisation, the resulting problem can be converted into a standard second-order cone programming problem which may be solved efficiently using modern optimization algorithms (for example, the primal-dual interior point method [START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF]). Typical advantages of such a solution strategy include the possibility of analyzing the existence, uniqueness, sensitivity and stability of the solution [START_REF] Sivaselvan | Complementarity framework for non-linear dynamic analysis of skeletal structures with softening plastic hinges[END_REF], the natural treatment of the singularities in the Mohr-Coulomb and Drucker-Prager yield criteria [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF][START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF], the straightforward extension from single-surface plasticity to multi-surface plasticity [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF], and the straightforward implementation of contact between deformable and rigid bodies [START_REF] Zhang | Quasi-static collapse of twodimensional granular columns: insight from continuum modelling[END_REF][START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF]. The proposed solution algorithm is incorporated into the framework of the Particle Finite Element Method (PFEM) [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF][START_REF] Oñate | Possibilities of the particle finite element method for fluid-soil-structure interaction problems[END_REF][START_REF] Cremonesi | A Lagrangian finite element approach for the analysis of fluid-structure interaction problems[END_REF] for handling large deformation. The PFEM is a novel continuum approach suitable for simulating problems involving both solid-like and fluid-like behaviour [START_REF] Oñate | Possibilities of the particle finite element method for fluid-soil-structure interaction problems[END_REF][START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF][START_REF] Salazar | Numerical modelling of landslide-generated waves with the particle finite element method (PFEM) and a non-Newtonian flow model[END_REF]. It makes use of particles to represent the material, as in meshfree particle methods, but solves the governing equations via a standard finite element procedure. Consequently, the PFEM inherits both the solid mathematical foundation of the traditional FEM as well as the capability of meshfree particle methods for handling extremely large deformation and freesurface evolution.

The paper is organized as follows. Section 2 presents the governing equations for dynamic analysis of an elastoviscoplastic problem. An incremental mixed variational principle is then proposed in Section 3 for recasting the governing equations into a min-max problem. Finite element discretisation is performed in Section 4, and the Particle Finite Element Method is described briefly in Section 5. Numerical examples are given in Section 6, before conclusions are drawn in Section 7.

Governing equations for elastoviscoplasticity

    σ b u  (1) T in  εu  (2) u on  uu (3)
T t on  N σ t [START_REF] De Blasio | Flow models of natural debris flows originating from overconsolidated clay materials[END_REF] where σ and ε are the Cauchy stress and the strain, b is the body force, u is the

in a plane-strain case. A superposed dot represents differentiation with respect to time.

Assuming the material is elastoviscoplastic, the total strain rate ε can then be split into an elastic strain rate e ε and a viscoplastic strain rate vp ε e vp  ε ε ε [START_REF] Quinn | Development of progressive failure in sensitive clay slopes[END_REF] The elastic strain rate is determined through Hook's law as e  εσ [START_REF] Bernander | Downhill progressive landslides in long natural slopes: triggering agents and landslide phases modeled with a finite difference method[END_REF] where is the elastic compliance matrix. The material is elastic if the stress state is inside the yield domain, namely vp ( ) 0

F    σ ε 0 (8)
where F is the yield function. In contrast, stress states satisfying ( ) 0

F  σ
lead to a viscoplastic strain rate. The classical Bingham model is utilized in this paper for describing the rheological properties of the sensitive clay. Despite its simple form, it performs well for approximating the plastic flow behaviour of these soils, especially Canadian clays [START_REF] Locat | Viscosity, yield stress, remolded strength, and liquidity index relationships for sensitive clays[END_REF]. The total stress thus is rewritten as

vp   σ τ ε ( 9
)
where  is the viscosity coefficient, τ is the stress lying on the boundary of F so that ( ) 0 F  τ

, and the quantity  στ is called the overstress. The viscoplastic strain rate is also normal to the yield surface at τ : vp () F   τ ετ [START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF] where  is the rate of the non-negative plastic multiplier and  τ is the gradient operator. It is clear that the above elastoviscoplastic model reduces to the classical elastoplastic model in the limiting case of 0   .

Laboratory tests show that the undrained shear strength of a sensitive clay decreases with increasing plastic shear strain. For materials exhibiting softening/hardening behaviour, the yield criterion function is expressed by ( , ) 0

F   τ
, where  is a set of hardening/softening variables which relate to the viscoplastic strain in the form of () vp H   ε [START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF] Specifically, for the Tresca yield criterion, we have [START_REF] Germanovich | Dynamic growth of slip surfaces in catastrophic landslides[END_REF] where cohesion softening is adopted to capture the basic post-failure behaviour. Following [START_REF] Potts | Finite element analysis of progressive failure of Carsington embankment[END_REF][START_REF] Troncone | Numerical analysis of a landslide in soils with strain-softening behaviour[END_REF], strain-softening is accounted for by reducing the cohesion u c using a bilinear function (Figure 1 [START_REF] Puzrin | Simple criteria for ploughing and runout in post-failure evolution of submarine landslides[END_REF] in which ij  is the Kronecker delta.

Variational principle

Hellinger-Reissner variational principle

The Hellinger-Reissner (HR) variational principle is of a mixed kind. Unlike the principle of minimum potential energy, in which displacements are considered as the only master field, the Hellinger-Reissner variational principle treats both the displacements and the stresses as the master fields [START_REF] Reissner | On a variational theorem in elasticity[END_REF].

For an elastostatic boundary-value problem, the Hellinger-Reissner functional [START_REF] Reissner | On a variational theorem in elasticity[END_REF] may be expressed as

t T T T T T 1 ( , ) ( )d d d 2                σ u σ σ σ S u b u t u (14)
The stationary value for the Hellinger-Reissner functional cannot be shown to be an extremum. Instead, the point obtained by ( , ) 0

 σu is a saddle point and, consequently, the problem becomes one of a min-max optimisation:

1 min max ( )d d d 2              u σ σ σ σ S u b u t u (15)
where the internal work is maximised with respect to the stresses and the total potential energy is minimised with respect to the displacements.

Generalised Hellinger-Reissner variational principle

A generalised Hellinger-Reissner variational principle is proposed for incremental analysis of elastoviscoplasticity. The governing equations summarised in section 2 are first discretized in time using the standard  -method. Details of the time discretisation, as well as the resulting incremental equations, are documented in Appendix A.

We present here the generalized incremental Hellinger-Reissner variational principle for incremental elasto-viscoplastic analysis. As for elastostatics, the principle is expressed in the form of a min-max program:

t tt t F                                                                   u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u τ n+1 )0  ( 16
)
where r is a set of variables that can be interpreted as dynamic forces. To illustrate the equivalence between the program ( 16) and the incremental form of the governing equations presented in Appendix A, the Karush-Kuhn-Tucker (KKT) optimality conditions associated with ( 16) are now derived. Following [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF][START_REF] Krabbenhoft | An interior-point algorithm for elastoplasticity[END_REF], the inequality constraint is first converted into an equality by adding a positively-restricted variable n+1 s . Then, the inequality on n+1 s is represented by introducing a penalty term in the objective function:

n+1 t T T T T T 1 n+1 n ( , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T T T n n n+1 11 min max d ( )d ( )d 2 1 dd 2 1 ( ) ( )d ( )d 2 ( ) d d d ln d t tt t s                                                                        u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u n+1 n+1 subject to ( ) 0 Fs  τ ( 17
)
where  is a sufficiently small positive constant. The penalty term n+1 ln  s in the objective function imposes the non-negativity requirement on n+1 s naturally, and is known as a logarithmic barrier function. The Lagrangian associated with program ( 17)

is t n+1 n+1 n+1 n+1 T T T T T 1 n+1 n 1 2 T T T 3 n+1 n+1 n+1 T T T T n n n n ( , , , , , ) 1 1 d ( )d ( )d 2 11
d d ( ) ( )d 22
( )d ( ) d d d ln s t t tt s                                                                          u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u n+1 n+1 n+1 d ( ( ) )d Fs         τ (18)
The KKT optimality conditions are found by differentiating the above Lagrangian with respect to the optimisation variables, namely:

TT 1 n+1 n n+1 1 T 1 n+1 n t 1 1 in 1 ( ) on                           σ σ r b 0 u N σ σ t (19) T 3 nn n+1 ( ) ( ) ( ) in t t                   u σ σ τ σ τ 0 σ ( 20
) 3 n n G n+1 n+1 ( ) ( ) ( ) in t t F                  σ τ σ τ τ 0 τ (21) 2 n+1 n+1 in t         r u 0 r (22) n+1 n+1 ( ) 0 in Fs        τ (23) 1 n+1 n+1 n+1 0 in ss s                (24)
It is apparent that the KKT conditions ( 19)-( 22 . The essential boundary condition ( 40) is assumed to hold a priori, and thus is not reflected in the KKT conditions. From condition [START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF] we can also see that the newly introduced variables r are dynamic forces.

Material hardening/softening

The variational principle [START_REF] Puzrin | The growth of shear bands in the catastrophic failure of soils[END_REF] can also be extended to handle more complex models involving hardening/softening yield surfaces following [START_REF] Krabbenhoft | Computational Cam clay plasticity using second-order cone programming[END_REF]. More specifically, the min-max program considering material hardening/softening is expressed as:

n+1 T T T T T 1 n+1 n ( , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T 1 2 T T nn 1 1 min max d ( )d ( )d 2 1 dd 2 1 ( ) ( )d ( )d 2 1 ( ) d d d d 2 t t t t t                                                                          u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u t n+1 n+1 subject to ( , ) 0 F    τ (25)
The underlined term is the newly introduced one with t being a new constitutive modulus associated with hardening/softening. The according KKT condition related to the variable 

is 1 n+1 n+1 n+1 n+1 n+1 ( , ) 0 ( , ) tt FF                        ττ (26)
which is the hardening/softening law, i.e. the evolution law, for the variable 

          εε εε εε (27)
Since Eq. ( 26) cannot be brought to be equal to the actual hardening/softening law (27) using a constant modulus, we therefore use the following tangent modulus as in [START_REF] Krabbenhoft | Computational Cam clay plasticity using second-order cone programming[END_REF]

vp nn n vp nn ( , ) () ( , ) t F dH dF       τ τ ε ετ (28)
which is updated at the beginning of each time step. Such a treatment of material hardening/softening behaviour in mathematical programming has been used successfully for approximating the hardening/softening behaviour in the Cam clay model [START_REF] Krabbenhoft | Computational Cam clay plasticity using second-order cone programming[END_REF].

Finite element formulation

The min-max program [START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF] can now be discretized using finite elements. For the sake of convenience, an intermediate variable e  σ σ τ (overstress) is introduced, which enables the optimization problem [START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF]

                                                                 u σ τ σ r σ σ σ u σ u r r r u σ σ σ σ b u t u σ σ τ τ n+1 , ) 0   (29)
Using standard finite element notations, we have

( , ) 0, 1, 2, , jG F j N                              u σ τ σ r κ σ C σ u B σ u B σ r Dr u A r σ M σ σ f κ H κ u f σ σ τ τκ (31)
where

ee t T T T 2 T T T TT 3 e T T c T e n d , d , d , d , 1 d , d , d d , d
t t t t                                       σ σ u σ r r u r σσ uu σ C N N B B N D N N A N N M N N H N N f N b N t f N σ (32)
The yield conditions are enforced at all Gauss integration points, with

2 2 2 1 ˆ2 1 ŝubject to ˆˆ( , ) 0, 1, 2, , jG F j N                         σ τ σ r κ σ C σ r Dr σ M σ κ H κ σ f B σ B σ A r f 0 σ σ τ τκ (33)
The natural boundary condition [START_REF] Krabbenhoft | Granular contact dynamics using mathematical programming methods[END_REF] r , resulting in [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF] This is obviously the discretised form of the displacement boundary conditions [START_REF] Oñate | The Particle Finite Element Method -An Overview[END_REF].

jG F j N                          σ τ σ
d n+1 u n+1       EU E u 0 r
Interaction between a deformable body and a rigid surface can be achieved in a straightforward manner in the above program according to [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF]. The classical Coulomb model for frictional contact is adopted in this study, which is

0, 0, = 0, | | 0 NN g p pg qp    (36)
As shown in Figure 3, N g is the gap between the material and the rigid surface, p is the contact pressure which is positive corresponding to compression, q is the tangential stress, and  is the friction coefficient between the material and the surface. After enforcing the conditions in [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF] on finite element nodes, the principle reads eu n+1

T T eT e n+1 n+1 ˆˆ( , , , , , ) n is the number of potential contacts. The above program is the final optimization problem to be solved. While it may be solved in a number of ways using either general or specialized methods, it is transformed here into a second-order cone program (SOCP) and then resolved using the high performance optimization solver MOSEK [START_REF] De Blasio | Flow models of natural debris flows originating from overconsolidated clay materials[END_REF]. The transformation of programs of the same type as (37) into a SOCP is straightforward, and has been documented in detail in [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF][START_REF] Krabbenhoft | Computational Cam clay plasticity using second-order cone programming[END_REF]. The main operation is to recast the quadratic terms in the objective function to linear ones, subject to a quadratic constraint, and to reform the yield function as a cone. Due to the attractive advantages presented in the introduction, a variety of mechanics problems have been formulated and solved in such a manner, including computational limit analysis of solids and plates [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF][START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF], static/dynamic analysis of elastoplastic frames and solids [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF][START_REF] Yonekura | Second-order cone programming with warm start for elastoplastic analysis with von Mises yield criterion[END_REF][START_REF] Zhang | Numerical investigation of the cylinder movement in granular matter[END_REF], analysis of steady-state non-Newtonian fluid flows [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF], consolidation analysis [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF], and the analysis of granular contact dynamics [START_REF] Huang | Three-dimensional granular contact dynamics with rolling resistance[END_REF][START_REF] Krabbenhoft | Granular contact dynamics with particle elasticity[END_REF][START_REF] Krabbenhoft | Granular contact dynamics using mathematical programming methods[END_REF].

T eT c d T u n+1 0 =1 T T u c e T 1 n+1 n+1 n+1 n 1 e T T 1 1 1 ˆˆˆmin 2 2 2 1 ˆˆ() 2
1 ˆŝubject to ˆˆ= , = 1, , = n c jj j k k c k gp p k n q                              σ τ σ r κ r σ C σ r Dr σ M σ κ H κ σ f EU r B σ A r Er E ρ f B σ σ σ τ n ρ n ρ * n+1 n+1 , | | 0 ( , ) 0, 1, 2, ,

Particle Finite Element Method

The Particle Finite Element Method (PFEM) is a Lagrangian approach capable of handling general large deformation problems without any real limitation on the magnitude of the deformation [START_REF] Oñate | Possibilities of the particle finite element method for fluid-soil-structure interaction problems[END_REF][START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF][START_REF] Oñate | The Particle Finite Element Method -An Overview[END_REF][START_REF] Cremonesi | A Lagrangian finite element approach for the simulation of water-waves induced by landslides[END_REF]. Its major characteristic is to treat mesh nodes as 'particles' that can move freely, and even separate from, the computational domain to which they originally belong. The basic steps of the utilized PFEM are summarized (see also Figure 4) in the following, with more details given in [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF]:

(1) Suppose that we have a cloud of particles, n C , at time n t ;

(2) Identify the computational domain using the  -shape method [START_REF] Edelsbrunner | #252, and cke, Three-dimensional alpha shapes[END_REF] on the basis of n C ;

(3) Create a finite element mesh, n M , through a triangulation of the recognized domain and discretize governing equations on n M ;

(4) Map the state variables such as stresses, strains, velocities, etc. from the old mesh, n-1 M , to the new mesh, n M ;

(5) Solve the discrete governing equations on the new mesh, n M , through a standard finite element procedure;

(6) Update the position of mesh nodes to arrive at n+1

C and repeat.

To date, a number of challenging problems involving large deformation and free-surface evolution have been tackled by the PFEM. These include the modelling of granular flows [START_REF] Zhang | Quasi-static collapse of twodimensional granular columns: insight from continuum modelling[END_REF][START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF][START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF][START_REF] Dávalos | On the numerical modeling of granular material flows via the Particle Finite Element Method (PFEM)[END_REF][START_REF] Cante | PFEM-based modeling of industrial granular flows[END_REF], landslides [START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF][START_REF] Cremonesi | A basal slip model for Lagrangian finite element simulations of 3D landslides[END_REF], landslide-generated waves [START_REF] Salazar | Numerical modelling of landslide-generated waves with the particle finite element method (PFEM) and a non-Newtonian flow model[END_REF][START_REF] Cremonesi | A Lagrangian finite element approach for the simulation of water-waves induced by landslides[END_REF], multi-fluid flows [START_REF] Becker | A unified monolithic approach for multi-fluid flows and fluid-structure interaction using the Particle Finite Element Method with fixed mesh[END_REF][START_REF] Idelsohn | The Particle Finite Element Method for Multi-Fluid Flows[END_REF][START_REF] Idelsohn | Multi-fluid flows with the Particle Finite Element Method[END_REF], fluid-structure interaction [START_REF] Oñate | Possibilities of the particle finite element method for fluid-soil-structure interaction problems[END_REF][START_REF] Zhu | Improved fractional step method for simulating fluid-structure interaction using the PFEM[END_REF][START_REF] Cremonesi | A Lagrangian finite element approach for the analysis of fluid-structure interaction problems[END_REF], soil-structure interaction [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF][START_REF] Zhang | Numerical investigation of the cylinder movement in granular matter[END_REF], bubble dynamics [START_REF] Mier-Torrecilla | Advances in the simulation of multi-fluid flows with the particle finite element method. Application to bubble dynamics[END_REF], the melting and spreading of polymers [START_REF] Oñate | Melting and spread of polymers in fire with the particle finite element method[END_REF], industrial forming processes, and the flow of fresh cement [START_REF] Cremonesi | Simulation of the flow of fresh cement suspensions by a Lagrangian finite element approach[END_REF]. In this paper, the solution algorithm for elastoviscoplastic analysis with strain softening is incorporated into the PFEM for progressive failure analysis of sensitive clays. It is notable that the governing equations proposed are on the basis of the infinitesimal strain theory which may lead to several errors for large deformation analysis.

The most serious one is the generation of strains as a result of rigid body motion. However, it has been shown in [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF][START_REF] Zhang | Particle finite element method in geomechanics[END_REF]] that this and related errors are relatively minor when the time steps used are small. As such, the price to pay for the convenience of being able to operate with usual infinitesimal strain theory appears to be very small. Indeed, such a strategy has been verified against analytical solutions for penetration problems [START_REF] Hu | A practical numerical approach for large deformation problems in soil[END_REF] and validated qualitatively as well as quantitatively against both quasi-static and dynamic collapse of a granular column [START_REF] Zhang | Quasi-static collapse of twodimensional granular columns: insight from continuum modelling[END_REF][START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF] and the penetration of shallowly embedded pipelines [START_REF] Tian | A simple implementation of RITSS and its application in large deformation analysis[END_REF].

Furthermore, it succeeds in reproducing a real-world flow-like landslide [START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF].

Numerical Examples

This section discusses numerical results for progressive failure analysis of sensitive clays using the proposed approach. Note that finite element analysis of strain-softening materials encounters issues of mesh sensitivity when using rate-independent models because the field equations that describe the motion of the body may lose hyperbolicity. Indeed, the corresponding boundary-value problem becomes ill-posed, with pathologically meshdependent solutions in which the width of the shear bands depends on the mesh size. The application of rate-dependent models is an effective way to circumvent this problem. It has been shown that viscous terms introduce a length scale effect into the initial boundary-value problem, even the rate-dependent model does not explicitly contain a parameter with the dimension of length [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Prevost | Dynamic strain localization in elasto-(visco-)plastic solids, part 2. plane strain examples[END_REF]. Consequently, viscoplastic models result in solutions where the shear bands have a finite width when strain localization occurs. It should be noted, however, that the main objective of this work is to capture the entire failure process in sensitive clays involving large deformation, rather than to predict the thickness of localized shear bands. As noted by Moore [START_REF] Moore | Numerical models for evaluating progressive failure in earth structures-A review[END_REF], the typical thickness of a shear band in clay at failure is between 0.01 and 2 cm, and thus it is impractical to predict both the microscopic and macroscopic soil response using a purely continuum model where a large earth structure is considered. One possible way of accounting for the responses on both the macro and micro levels is through the multiscale computational modelling technique [START_REF] Andrade | Multiscale modeling and characterization of granular matter: From grain kinematics to continuum mechanics[END_REF][START_REF] Li | A generalized Hill's lemma and micromechanically based macroscopic constitutive model for heterogeneous granular materials[END_REF][START_REF] Guo | A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[END_REF][START_REF] Li | Micro-macro homogenization of gradient-enhanced Cosserat media[END_REF].

One-dimensional elasto-viscoplastic problem

To verify the proposed variational principle, we consider an axial bar subject to a prescribed load (Figure 5(a)). The material is represented by a one-dimensional elastoviscoplastic model (Figure 5(b)). If the mass of the bar is sufficiently small, so that any induced inertial forces are negligible, the load produces a uniform stress and strain along the bar and an analytical solution is available. In the following, the material parameters of the bar are assumed to be:

Young's modulus yield stress strengths can be interpreted as strengths of a material at undisturbed and remoulded states. The simulation is conducted using a total of 20, 30, and 40 time increments, respectively, and again the agreement between the numerical and analytical solutions is satisfactory (Figure 9).

Collapse of a sensitive clay column

As the second example, we consider the collapse of a sensitive clay column (Figure 10) in a container which is 50 cm wide and 100 cm high. The container is lifted up quickly leading to the spread of the sensitive clay. Such an experimental test has been widely used for investigating the behaviour of granular matter [START_REF] Lube | Collapses of two-dimensional granular columns[END_REF][START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Balmforth | Granular collapse in two dimensions[END_REF][START_REF] Thompson | Granular column collapses: further experimental results[END_REF][START_REF] Lacaze | Planar collapse of a granular column: Experiments and discrete element simulations[END_REF], but has also been adopted for studying the quickness of sensitive clays [START_REF] Thakur | Quickness of sensitive clays[END_REF]. Here, the problem is considered to deform under planestrain conditions and only half of the geometry is modelled due to the symmetry. The material parameters are as follows: Young's modulus . The frictional coefficient between the clay and the rigid surface is taken as 0.3 and the gravitational acceleration 9.8 g  m/s 2 . The column is discretized using 7,962 6-node triangular elements with [START_REF] Puzrin | The growth of shear bands in the catastrophic failure of soils[END_REF]199 nodes, and the time step utilized is

0.01 t  s.
The collapse procedure of the column obtained from the simulation is illustrated in Figure 11, in which the colour is proportional to the accumulated equivalent plastic strain. The normalized time t refers to clay is also simulated for comparison (Figure 11). As shown, the collapse mechanism for this case is quite different to the previous case of an undisturbed sample. Rather than fail progressively, nearly all the material experiences plastic deformation with the material near the bottom possessing the maximum equivalent plastic strain. Figure 12 shows the curves of the front location and centre height against normalized time for columns of both initially undisturbed and remoulded sensitive clays. The collapse of the column of remoulded clay results in a final deposit with a much smaller height and considerably larger length. Both the sensitive and remoulded clay columns reach their maximum run-out distance at around 2.15 t  . However, the final centre height for the remoulded case is obtained earlier ( ). Notably, the final centre height and length are obtained simultaneously for the undisturbed sample, which differs from that for the remoulded sample. The final mesh topology for both cases, illustrated in Figure 13, verifies that the proposed approach can handle the extreme mesh distortion that accompanies failure for this problem. A video of the collapse of both the remoulded and undisturbed sensitive clay columns is provided in the supplementary materials.

To estimate the mesh sensitivity, the collapse of the initially undisturbed clay was also reanalyzed using three different mesh sizes, where the length of the element edge was set to 1.0 h  cm (1,934 triangles ), 0.75 cm (3,528 triangles) and 0.5 cm (7,962 triangles). The cases are referred to as coarse, medium, and fine meshes. As shown in Figure 14, the curves of locations against time for all three tests agree well with each other. Moreover, all three simulations result in very similar final deposits and shear bands (Figure 15), which proves that a further decrease in the mesh size will not alter the form of the predicted failure mode.

Retrogressive collapse of a slope in sensitive clay

As observed in Scandinavia and eastern Canada, a fast and significant retrogressive collapse of a slope in sensitive clay may be triggered by a small initial slide [START_REF] Locat | Progressive failures in eastern Canadian and Scandinavian sensitive clays[END_REF]. To illustrate the ability of the proposed approach for modelling such a geohazard, we consider the sensitive clay deposit shown in Figure 16. Here, collapse is triggered by removing a rigid triangular block at the toe of the slope (which may be caused by erosion or excavation). The material parameters of the sensitive clay are as follows: Young's modulus . The frictional coefficient between the sensitive clay and the rigid bottom surface is set to 0.1 and the gravitational acceleration 9.8 g  m/s 2 . A total of 18,420 6-noded triangular elements (37,355 mesh nodes) is used to discretize the initial computational domain. The time step is 0.025 t  s and the simulation proceeds until the final deposit is obtained.

The retrogressive failure process from the simulation is illustrated in Figure 17, with the colour being proportional to the accumulated equivalent plastic strain. As illustrated, the erosion leads to the first retrogressive collapse C1 in the slope (Figure 17(b)). Two shear bands initiate from the bottom and propagate towards the top surface and the front inclined surface, respectively, resulting in a graben. During the sliding, one more shear band is generated in the graben dividing it into two elastic parts (Figure 17(c) and(d)). As the disturbed mass due to collapse C1 moves far away from the new slip surface, the second retrogressive collapse C2 occurs (Figure 17(e) and(f)). This mimics the first criterion for the occurrence of retrogressive failure, that the slide debris should be able to flow out of the slide area [START_REF] Thakur | Quickness of sensitive clays[END_REF]. The same as that in C1, two plastic shear bands, also originating from the slope base, are formed in C2 which results in a graben and a horst. The mass in front of the new slope surface continues to move forward leading to the third retrogressive failure of the slope, C3 (Figure 17(g)). After this, a considerable amount of mass is deposited in front of the new slope surface, resisting further collapse (Figure 17(h)). The final configuration of the slope, shown in Figure 17(i), indicates that most of the clay involved in C2 and C3 has been remoulded. This reflects the other criterion for the occurrence of retrogressive failure, which states that the slide debris should be completely remoulded [START_REF] Thakur | Quickness of sensitive clays[END_REF]. Eventually, the retrogressive failure results in a deposit with a run-out distance of 28.71 m and a retrogression distance of 14.76 m (Figure 17(i)). A video of the entire failure process of the slope is available in the supplementary materials.

The velocity of the sliding front and the maximum velocity were also recorded and are depicted in Figure 18. As illustrated, the maximum velocity is not always located at the sliding front. This can be explained by examining the velocity contour (Figure 19). The first retrogressive failure results in the transformation of a part of the gravitational potential energy of the mass into kinetic energy, with the remaining energy being dissipated by plastic shearing (Figure 19(a)). The sliding front thus possesses the maximum velocity due to this transformation. The second retrogressive collapse further releases potential gravitational energy (Figure 19(b) and pushes the materials in front of it, consequently increasing the velocity of the corresponding mass. In contrast, the velocity of the sliding front decreases because of the friction along the basal surface and the effects of plastic dissipation. After a while, the mass at the middle moves faster than the sliding front does as shown in Figure 19(c). The third retrogressive collapse further releases potential gravitational energy (Figure 19(d)). The velocity of the involved mass in collapse C3 is relatively low, however, because a considerable body of clay with low velocity is located in front of the new slope surface. Note that the sliding front already ceased at 11.4 t  s (Figure 19(e)). However, the clay at the middle of the sliding mass continues to be pushed and moves forward, which eventually disturbs the sliding front (Figure 19(f) and (g)).

The value of viscosity, back-calculated from various subaerial and submarine slides by Edgers and Karlsrud [START_REF] Edgers | Soil flows generated by submarine slides -case studies and consequences[END_REF] and Johnson and Rodine [START_REF] Johnson | Debris flow[END_REF], is in the range of 100 to 1499 Pa s  .

We here investigate the effect of the viscosity on the retrogressive failure by analysing the problem with varying viscosity coefficients of Pa s  (Figure 20). This means that an increase of  ( 10   Pa s  ) causes a decrease in the run-out and retrogressive distance, as illustrated in Figure 21, because a higher viscosity results in much more plastic dissipation and, consequently, less potential gravitational energy is converted into kinetic energy. Such a decrease in energy transformation not only leads to a smaller run-out distance, but also causes more clay to be deposited in front of the newly formed slope surface which resists the occurrence of further retrogressive failure.

Conclusions

The progressive failure process of sensitive clays is simulated using the Particle Finite Element Method with an advanced elastoviscoplastic model which is a combination of the Bingham model (for describing rheological behaviour) and the Tresca model with strainsoftening (for capturing progressive failure behaviour). The resulting elastoviscoplastic analysis is reformulated as an optimization problem on the basis of a mixed variational principle and resolved in mathematical programming.

The proposed formulation is verified against the analytical solution of a one-dimensional elastoviscoplastic problem. The capability of the proposed computational approach for modelling progressive failure is illustrated by simulating the collapse of a column of sensitive clay. Additionally, the retrogressive failure of a slope in sensitive clay is reproduced successfully. The simulation results reflect the essential conditions for the occurrence of retrogressive collapse which are that the slide debris should be fully remoulded and flow away from the slide area. Furthermore, the effect of the viscosity of a sensitive clay on the nature of retrogressive collapse is also studied. Numerical results show that an increase of viscosity leads to a decrease in both the run-out distance and the retrogression distance due to the dissipation of a large amount of energy.

Although the problem in this study is simulated under undrained conditions using total stresses, an effective stress analysis can also be performed. This can be achieved by merging the SOCP formulation for consolidation analysis of saturated porous media introduced in [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF],

where rate-independent models are utilised, with the mixed variational principle presented in this paper. The resulting formulation then can be incorporated into the particle finite element method for investigating the hydro-mechanical mechanism in the progressive failure of sensitive clays.

FF       ττ (50)
In summary, the governing equations for incremental analysis of elastoviscoplasticity consist of Eqs. ( 40), ( 41), ( 44)-( 50).

  ) are equivalent to the corresponding incremental equations presented in Appendix A. The last two conditions recover the yield condition and the complementarity condition shown in (

  e σ , τ , r , û , and κ are vectors containing the values of the corresponding field variables at interpolation points, N is a matrix consisting of shape functions, and T uu  BN . The mixed finite element shown in Figure 2 is adopted in this study, where the distribution of the interpolation points for the different variables is depicted. Substituting the above equations into the program (29) leads to

  points. After solving the minimization part over u  of program[START_REF] Potts | Finite element analysis of progressive failure of Carsington embankment[END_REF], and transforming the maximum into a minimum with an opposite sign, we obtain

  the essential boundary condition (44) are enforced, E is an index matrix consisting of entries equal to 0 and 1, d U is a vector consisting of the prescribed displacements at mesh nodes, and the newly introduced u n+1 r represents the nodal reaction force. As shown, both the objective function and the constraints of program (34) are altered (the underlined terms) due to the imposition of the essential boundary conditions. The validity of the above can be checked by differentiating the Lagrangian associated with the program (34) with respect to u n+1

  the rigid boundary, c E is an index matrix of zeros and ones, and c

Y

  the proposed formulation to capture the strain-rate dependence of the stress response and stress relaxation behaviour is examined first. To this end, we set the prescribed strain increase at a constant rate  until time * 0.4 t  s and then hold the strain constant, leading to stress relaxation. The analytical solution of this problem is available[START_REF] Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF] and we consider three different load rates, namely  =0.2, 0.4, and 0.6, respectively, to produce a rate-dependent response (Figure6(a)). The yield stress is set to be constant ( 0 YY   ) in this case (Figure 6(b)) and the time increment is 3 2 10 t    s in all simulations. Figure7illustrates the simulated stress response for different load rates as well as the corresponding analytical solutions. For all cases, the resulting stresses increase in a stable manner until their maximum values are reached. A higher load rate results in a larger maximum stress reflecting the effect of viscosity. At the time * tt  , the stresses for all three cases drop sharply representing stress relaxation behaviour. Eventually, the residual stresses for all cases asymptote towards the initial stress strength of the material, 0 . All the simulated results agree with the analytical solutions, which verifies the proposed variational formulation and finite element implementation.We now consider the details of strain-softening behaviour. The prescribed strain in this case increases with a constant rate  =0.2 (Figure 8(a)); however the yield stress strength Y strain reaches 7% (Figure 8(b)). Such a phenomenon of reduction has been widely observed for materials undergoing plastic deformation. The initial and residual

  being the initial height of the column. For initially undisturbed sensitive clays, lifting the container results in two shear bands dividing the column into three parts (Figure11(a)). The upper part moves downward while the middle part, which is in the shape of a triangle, is pushed out horizontally. After a considerable movement of the middle part, a shear band is formed in the lower part (Figure11(b)) and then one more shear band appears in the upper part (Figure11(c)). The second shear band in the upper part deforms another layer of sensitive clay and the lower part, which was intact, is disturbed significantly because of the shear band formed (Figure11(d)). Further collapse of the column leads to two more layers being squeezed out (Figure11(e) and (f)). Localized shear bands can be observed clearly in the final deposit, with some parts of the column remaining undisturbed throughout the failure process. The collapse of remoulded sensitive

1   1   Pa s  and 10 1  

 11101 Pa s  , 10 Pa s  , 100 Pa s  , and 1000 Pa s  . Other material parameters for these simulations are the same as in the previous case. As shown in Figure 20, retrogressive failure occurs four times for both Pa s  . Thus, the retrogression distances for these two cases, 20.25 m and 20.13 m, are very close (Figure 21), although the run-out distance for Pa s  (38.18 m) is slightly larger than that for 10   Pa s  (36.73 m). When  is increased further, fewer retrogressive failures are induced: for example three times for 100   Pa s  and only twice for 1000  

Figure 2 .Figure 3 .

 23 Figure 2. The mixed triangular element used in the simulation

Figure 4 .Figure 5 .

 45 Figure 4. Steps for the Particle Finite Element Method (after [35])

Figure 6 . 34 Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .Figure 15 .Figure 16 .

 63478910111213141516 Figure 6. Curves of (a) the applied strain rate and (b) the variation of the yield stress for the one-dimensional stress relaxation problem

Figure 17 . 40 Figure 18 .Figure 19 .

 17401819 Figure 17. Retrogressive failure procedures of the slope. Colours are proportional to accumulated equivalent plastic strain

  has been included through the terms e

	problem of program (33). More specifically, the program turns out to be
	eu ˆˆ( n+1 , , , , , )	1 2	T				1 2	T n+1 ˆˆˆmin eT n+1 1 2	e
		T 1 ˆˆ() eT c d T u n+1 2
	to	T	n+1		T	u n+1 ˆŝubject e n+1 1	1	T	n
									1
		e ˆˆ(
			n+1	,	n+1	) 0,	1, 2, ,
	, whereas the imposition of essential boundary conditions for the displacements requires the introduction n+1 f of a new variable u r since the displacement increment u is a field variable for the dual

Acknowledgements

The authors wish to acknowledge the support of the Australian Research Council Centre of Excellence for Geotechnical Science and Engineering and the Australian Research Council Discovery Project funding scheme (Project Number DP150104257).

Appendix A. Time discretization

The momentum conservation equation ( 1) is first discretized in time using the standard  - method as:

2 n

(1 [START_REF] Yonekura | Second-order cone programming with warm start for elastoplastic analysis with von Mises yield criterion[END_REF] where v are velocities, 1

 and 2  are parameters taking values in [0, 1], the subscripts n and n+1 refer to the known and new, unknown states, and

Rearranging the above equations leads to

with the displacement increments

The natural boundary condition is approximated in an analogous manner leading to

and the discretised essential boundary condition is


Ana G Lecuona

COMPLEMENTARY LEGS AND RATIONAL BALLS

In this note we study the Seifert rational homology spheres with two complementary legs, i.e. with a pair of invariants whose fractions add up to one. We give a complete classification of the Seifert manifolds with 3 exceptional fibers and two complementary legs which bound rational homology balls. The result translates in a statement on the sliceness of some Montesinos knots.

introduction

Seifert manifolds are a particularly well understood set of 3-manifolds. They are oriented, closed 3-manifolds admitting a fixed point free action of S 1 and they are classified by their Seifert invariants [START_REF] Orlik | On 3-manifolds with local SO(2) action[END_REF]. They include all lens spaces. They were first introduced in [START_REF] Seifert | Topologie Dreidimensionaler Gefaserter Räume[END_REF] and constitute one of the building blocks in the JSJ decomposition. They have been studied from the most various perspectives and are partially or completely classified from many different points of view: Stein fillability, symplectic fillability, tightness etc.

A relevant and difficult question one can ask about a 3-manifold is whether or not it bounds a rational homology ball [START_REF] Kirby | Problems in low dimensional manifold theory, Algebraic and geometric topology[END_REF]. In the case of Seifert manifolds there are some partial answers spread out in the literature [START_REF] Casson | Some homology lens spaces which bound rational homology balls[END_REF][START_REF] Greene | The slice-ribbon conjecture for 3-stranded pretzel knots[END_REF][START_REF]On the slice-ribbon conjecture for pretzel knots[END_REF][START_REF] Aceto | Dehn surgeries and rational homology balls[END_REF] and there is a complete classification for the lens space case [START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF].

Seifert manifolds can be realized as the double cover of S 3 branched over Montesinos links [START_REF] José | Seifert manifolds that are ramified two-sheeted cyclic coverings[END_REF]. If a Seifert manifold does not bound a rational homology ball, then it is a well known fact [START_REF] Eisenbud | Three-dimensional link theory and invariants of plane curve singularities[END_REF]Lemma 17.2] that the corresponding Montesinos link does not bound an embedded surface in B 4 with Euler characteristic equal to one. In the case of the link being a knot, this tells us that the knot is not slice. The famous slice-ribbon conjecture states that every slice knot bounds an immersed disk in S 3 with only ribbon singularities. A surface singularity is ribbon if it is the identification of two segments, one contained in the interior of the surface and the other joining two different boundary points.

Among the Seifert fibered rational homology spheres Y with Seifert invariants Y = Y (b; (α 1 , β 1 ), . . . , (α n , β n )), b ∈ Z, α i > β i > 0, we will be interested in this note in those having a pair of complementary legs, that is, a pair of invariants verifying βi αi + βj αj = 1. The case n ≤ 3 is of particular interest since for n ≤ 3 there is a one to one correspondence between the set of Seifert manifolds and the set of the corresponding Montesinos links and our results will have an interpretation also in terms of those. Our main result is the following. The list L in the theorem is completely explicit and can be found in Section 3.1. Given any set of Seifert spaces with at most 3 exceptional fibers we will simply say that a Montesinos link belongs to this set if the double cover of S 3 branched over this link belongs to it. In terms of Montesinos links we have the following result.

Theorem 1.2. Each Montesinos link in the set L is the boundary of a ribbon surface Σ such that χ(Σ) = 1.

An immediate corollary of this theorem is that the slice-ribbon conjecture is true inside the family of Montesinos knots whose double covers are Seifert fibered spaces with 3 exceptional fibers and two complementary legs.

The proof of Theorem 1.1 relies very heavily on the analysis previously done in [START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF][START_REF] Lecuona | On the slice-ribbon conjecture for Montesinos knots[END_REF]. In this note we first show that a Seifert space Y with complementary legs is rational homology cobordant to a lens space L. This allows us to study a lattice embedding problem in a very efficient way, which yields a relationship between the Seifert invariants of Y and of L. The classification of lens spaces bounding rational homology balls in [START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF] allows us to make the list L explicit.

1.1. Organization of the paper. In Section 2 we introduce the strict minimum regarding Seifert manifolds, the language of lattice embeddings and Montesinos links. In Section 3 we carry out the proofs of the theorems in the introduction and establish the list L . Finally in Seciton 4 we collect some noteworthy technical remarks about the lattice embedding problem.

Notation and Conventions

2.1. Graphs and Seifert spaces. It is well known [START_REF] Walter | Seifert manifolds, plumbing, µ-invariant and orientation reversing maps, Algebraic and geometric topology[END_REF] that, at least with one of the two orientations, a Seifert fibered rational homology sphere defined by the invariants Y = Y (b; (α 1 , β 1 ), (α 2 , β 2 ), (α 3 , β 3 )) is the boundary of a 4-manifold obtained by plumbing according to the graph:

-a 1,2 -a n2,2 • . . . • -a n1,1 -a 1,1 -a 0 • . . . • • -a 1,3 -a n3,3 • . . . • where α i β i = [a 1,i , ..., a ni,i ] := a 1,i - 1 a 2,i - 1 . . . a ni-1,i - 1 a ni,i , a j,i ≥ 2 and a 0 = -b ≥ 1.
This graph is unique and will be called the canonical graph associated to Y . We will use the expression three legged canonical graph to denote any such graph.

If we are given a weighed graph P , like for example the one above, we will denote by M P the 4-manifold obtained by plumbing according to P and by Y P its boundary. The incidence matrix of the graph, which represents the intersection pairing on H 2 (M P ; Z) with respect to the natural basis, will be denoted by Q P .

The number of vertices in P coincides with b 2 (M P ) and we will call (Z b2(MP ) , Q P ) the intersection lattice associated to P . The vertices of P , which from now on will be identified with their images in Z n and will also be called vectors, are indexed by elements of the set J := {(s, α)| s ∈ {0, 1, ..., n α }, α ∈ {1, 2, 3}}. Here, α labels the legs of the graph and L α := {v i,α ∈ P | i = 1, 2, ..., n α } is the set of vertices of the αleg. The string associated to the leg L α is the n α -tuple of integers (a 1,α , ..., a nα,α ), where a i,α := -Q P (v i,α , v i,α ) ≥ 0. The three legs are connected to a common central vertex, which we denote indistinctly by v 0 = v 0,1 = v 0,2 = v 0,3 (notice that, with our notation, v 0 does not belong to any leg).

Let (Z n , -Id) be the standard negative diagonal lattice with the n elements of a fixed basis E labeled as {e j } n j=1 . As an abbreviation in notation we will write e j • e i to denote -Id(e j , e i ). If the intersection lattice (Z b2(MP ) , Q P ) admits an embedding ι into (Z n , -Id), we will then write P ⊆ Z n and will omit the ι in the notation, i.e. instead of writing ι(v) = j x j e j we will directly write v = j , x j e j .

We will adopt the following notation: for each i ∈ {1, ..., n}, (s, α) ∈ J and S ⊆ P we define

V S := {j ∈ {1, ..., n} | e j • v s,α = 0 for some v s,α ∈ S}.
Given e ∈ Z n with e • e = -1, we denote by π e : Z n → Z n the orthogonal projector onto the subspace orthogonal to e, i.e.

π e (v) := v + (v • e)e ∈ Z n , ∀v ∈ Z n .

Complementary legs.

The key concept in this paper is that of complementary legs. As explained in the introduction, they are pairs of legs in the graph P associated to Seifert invariants with the property βi αi + βj αj = 1. Suppose that in the graph displayed at the beginning of Section 2.1 the legs L 2 and L 3 were complementary. It is well known then that the strings of integers (a 1,1 , ..., a n1,1 ) and (a 1,2 , ..., a n2,2 ) are related to one another by Riemenschneider's point rule [START_REF] Riemenschneider | Deformationen von Quotientensingularitäten (nach zyklischen Gruppen)[END_REF] and that there is an embedding of L 2 ∪ L 3 ⊆ Z n2+n3 [START_REF]Sums of lens spaces bounding rational balls[END_REF]. There essential features of this embedding that we will need in this note are:

• There exists a sequence of contractions (see [START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF] for the definition and general theory of embeddings of linear sets) to the legs L2 and L3 with associated strings (2) and (2). The embedding L2 ∪ L3 ⊆ Z 2 is unique and given by e 1 + e 2 and e 1 -e 2 . • The basis vector e 1 appears only twice in the embedding

L 2 ∪ L 3 ⊆ Z n2+n3 .
• Given a linear weighed graph L with string (b 1 , . . . , b k ) and a embedding L ⊆ Z k , it is possible to combine L and two complementary legs to yield a 3 legged graph P L with a embedding in Z k+n2+n3 . Indeed, the three legs will be the complementary legs L 2 and L 3 , and the leg L 1 with associated string (b 1 , . . . , b k-1 ); the central vertex will have weight b k + 1 and the embedding is the obvious one, adding -e 1 to the embedding of the vertex in L with weight b k . • The embedding in Z k+n2+n3 of the graph P L constructed in the preceding point can be contracted to an embedding of P c L ⊆ Z k+2 which is a graph with L2 and L3 as complementary legs.

A final notion that will appear in the text and is related to complementary legs is that of the dual of a string of integers (a 1 , . . . , a n ), a i ≥ 2. It refers to the unique

string (b 1 , . . . , b m ), b i ≥ 2 such that 1 [a 1 , . . . , a n ] + 1 [b 1 , . . . , b m ] = 1.
2.3. Montesinos links. Seifert spaces admit an involution which presets them as double covers of S 3 with branch set a link. This link can be easily recovered from any plumbing graph P which provides a surgery presentation for the Seifert manifold Y P . As shown in Figure 1 the Kirby diagram obtained from P is a strongly invertible link. The involution on Y P , which is a restriction of an involution on M P , yields as quotient S 3 , boundary of B 4 , and the branch set consists of a collection of twisted bands plumbed together. The twists correspond to the weights in the graph and the plumbing to the edges. The Montesinos link is by definition the boundary of the surface obtained by plumbing bands [START_REF] José | Seifert manifolds that are ramified two-sheeted cyclic coverings[END_REF]. An example is shown in Figure 1. Given a three legged graph P , there is a unique Montesinos link ML P which can be obtained by the above procedure. If we start with a linear graph P ′ , then the same procedure yields a 2-bridge link, which we will denote by K(p, q) where Y P ′ is the lens space L(p, q).

The list L and consequences

The key feature of plumbing graphs with complementary legs is the following proposition. A similar proof in a more general setting can be found in [1, Proposition 4.6] (cf. [9, Remark 6.5]). Proposition 3.1. Let P be a 3 legged canonical graph with two complementary legs L 2 and L 3 and such that the string associated to L 1 ∪ {v 0 } is of the form (a n1,1 , ..., a 1,1 , a 0 ). Then, there is a rational homology cobordism between Y P and the 3-manifold Y P ′ associated to the linear graph with string (a n1,1 , ..., a 1,1 , a 0 -1).

Proof. Since the graph P has two complementary legs, one can attach a 4 dimensional 2-handle h to M P obtaining a 4-manifold whose boundary is (S 1 ×S 2 )#Y P ′ , where P ′ is a linear graph with associated string (a n1,1 , a n1-1,1 , ..., a 1,1 , a 0 -1) [8, Lemma 3.1]. Notice that the dual 2-handle h, attached to (S 1 × S 2 )#Y P ′ , kills the free part of H 1 ((S 1 × S 2 )#Y P ′ ; Z), since Y P = ∂M P is a rational homology sphere. Moreover, attaching a 3-handle k to (S 1 × S 2 )#Y P ′ along one of the essentials spheres of the summand S 1 × S 2 , the free part of H 2 ((S 1 × S 2 )#Y P ′ ; Z) is killed and one obtains Y P ′ . Dually, attaching a 1-handle k to Y P ′ yields (S 1 × S 2 )#Y P ′ .

It follows that, if Y P ′ bounds a rational homology ball W , then Y P is the boundary of the 4-manifold X := W ∪ k ∪ h, where the attaching sphere of h has non-trivial algebraic intersection with the belt sphere of k, since h kills the free part of H 1 ((S 1 × S 2 )#Y S ; Z). We claim that X is a rational homology 4-ball. In fact, since W is a rational homology ball, we have

H 4 (X; Q) = H 3 (X; Q) = 0.
The rest of the homology groups can be easily computed using cellular homology with rational coefficients: calling X m the m-skeleton of the CW-complex X, the cellular chain complex is

• • • / / H 3 (X 3 , X 2 ; Q) d3 / / H 2 (X 2 , X 1 ; Q) d2 / / H 1 (X 1 , X 0 ; Q) d1 / / Q .
Since the attaching sphere of h ∈ H 2 (X 2 , X 1 ; Q) has non-trivial algebraic intersection with the belt sphere of k ∈ H 1 (X 1 , X 0 ; Q), by definition of the boundary operator d 2 , we have h ∈ Ker d 2 , and since W is a rational homology ball it follows

∼ = (a) (b) 2 2 -2 -2 3 3 -1 -1 -1 -1 u Figure 1. Part (a)
shows a 3 legged plumbing graph and its associated Kirby diagram as a strongly invertible link. The bottom picture shows the branch surface in B 4 which is the image of Fix(u). This surface is homeomorphic to the result of plumbing bands according to the initial graph: the gray lines retract onto the sides of the rectangle. that H 2 (X; Q) = 0. Finally, we assume, without loss of generality, that X is connected and hence d 1 ≡ 0. Therefore, we have

H 1 (X; Q) = H 1 (X 1 , X 0 ; Q)/im d 2 . The group H 1 (X 1 , X 0 ; Q), which is a direct sum of finitely many Q, has one more Q summand, coming from the 1-handle k, than H 1 (W 1 , W 0 ; Q). However, k ∈ im d 2 and therefore H 1 (X; Q) = H 1 (W ; Q) = 0.
In fact, an appropriate decomposition of X and W gives the following commutative diagram.

H 2 (W 2 , W 1 ; Q) d2 / / i * H 1 (W 1 , W 0 ; Q) d1 / / i * 0 H 2 (X 2 , X 1 ; Q) d2 / / H 2 (X 1 , X 0 ; Q)
Note that the first row of the diagram is exact because W is connected and H 1 (W ; Q) = 0. Since the attaching sphere of h has non-trivial algebraic intersection with the belt sphere of k, there exists some α ∈ H 1 (W 1 , W 0 ; Q) such that d 2 ( h) = q k + i * (α), where q ∈ Q, q = 0. The exactness of the first row implies that

α = d 2 (β) for some β ∈ H 2 (W 2 , W 1 ; Q). Therefore d 2 ( h -i * (β)) = d 2 ( h) -d 2 i * (β) = q k + i * (α) -i * d 2 (β) = q k,
and, as claimed, X is a rational homology ball.

Vice versa, if Y P bounds a rational homology ball Z, then Y P ′ is the boundary of a 4-manifold V := Z ∪ h ∪ k, where this time the attaching sphere of k has non-trivial algebraic intersection with the belt sphere of h. It follows, arguing as before, that V is a rational homology ball. We conclude that there exists a rational homology cobordism between Y P and Y P ′ .

If P is a canonical 3 legged graph and Y P is the boundary of a rational homology ball W , then there is an embedding of the lattice associated to H 2 (M P , Z) into the standard negative lattice of the same rank. Indeed, since P is canonical the intersection form Q P of M P is negative definite [13, Theorem 5.2] and hence, X P := M P ∪ YP (-W ) is a closed smooth negative 4-manifold. By Donaldson's Theorem [START_REF] Donaldson | The orientation of Yang-Mills moduli spaces and 4-manifold topology[END_REF] the intersection lattice of X P is isomorphic to (Z n , -Id), where n = b 2 (X P ). Moreover, since W is a rational homology ball n = b 2 (M P ) and therefore, via X, there is an embedding of (H 2 (M P , Z), Q P ) into (Z n , -Id), which we simply write as P ⊆ Z P . In the next lemma we study the consequences of this obstruction when the graph P has two complementary legs. Lemma 3.2. Let P be a canonical 3 legged graph with two complementary legs L 2 and L 3 and suppose that Y P is the boundary of a rational homology ball. Then, the linear set L 1 ∪ {v 0 } has associated string of length n 1 + 1 of the form (a n1,1 , ..., a r,1 , 2 [r] ), where 0 ≤ r ≤ n 1 + 1 and the string (a n1,1 , ..., a r,1 -1) corresponds to a linear graph admitting an embedding in Z n1-r+1 and defining a lens space L P ′ .

Proof. Since Y P bounds a rational homology ball, the space Y P ′ from Proposition 3.1 does too and therefore it holds P ′ ⊆ Z P ′ . From Section 2.2 we know that it is always possible to extend the embedding of P ′ ⊆ Z P ′ to an embedding P ⊆ Z P in such a way that V L1 ∩V L2∪L3 = ∅. It follows that if Y P bounds a rational homology ball then there is an embedding P ⊆ Z P with the property V L1 ∩ V L2∪L3 = ∅ which allows us to contract the complementary legs to L 2 = {v 1,2 = e j + e k } and L 3 = {v 1,3 = e j -e k } for some j, k ∈ {1, ..., n}.

Since

v 0 • v 1,2 = v 0 • v 1,3 = 1, then j ∈ V v0 , |e j • v 0 | = 1 and k ∈ V v0 . If v ∈ L 1 was such that V v ∩ {e j , e k } = ∅, then we could not have v • v 1,2 = 0 and v • v 1,3 = 0, therefore V L1 ∩ V L2∪L3 = ∅. Since v 1,1 • v 0 = 1, we have V v0 ∩ V L1 = ∅ and hence v 0 = ṽ0 -e j with V ṽ0 ⊆ V L1 . Moreover, it follows that a 0 = -v 0 • v 0 ≥ 2.
If a 0 ≥ 3 then ṽ0 •ṽ 0 ≤ -2. By definiton of complementary legs |V L2∪L3 | = n 2 +n 3 and therefore P ′ n1+1 := L 1 ∪ {ṽ 0 } ⊆ Z n1+1 . Notice that in this case the r in the statement equals 0.

If a 0 = 2 let us call 1 ≤ r ≤ n 1 + 1 the bigest integer such that a 0,1 = ... = a r-1,1 = 2. It is not difficult to check that the r vectors v 0,1 , ..., v r-1,1 belong to the span of r + 1 basis vectors. If r = n 1 + 1, then |V L1∪{v0} | = n 1 + 2, by definition of complementary legs |V L2∪L3 | = n 2 + n 3 and V L1∪{v0} ∩ V L2∪L3 = {j}. Therefore every 3 legged canonical graph with two complementary legs L 2 and L 3 satisfying a j,1 = 2 ∀j ∈ {0, ..., n 1 } admits an embedding into Z n , where n = n 1 + n 2 + n 3 + 1.

Finally, if r < n 1 + 1 there exists i ∈ {1, ..., n} such that {i} = V vr,1 ∩ V vr-1,1 and since |V L1 ∩ V v0 | = 1 we have that

P ′ n1-r+1 := (L 1 \ r+1 j=2 v j-1,1 ) ∪ π ei (v r,1 ) ⊆ Z n1-r+1 .
Since v r-1,1 • v r,1 = 1, we have π ei (v r,1 ) • π ei (v r,1 ) = -a r,1 + 1 and the statement follows.

Combining Proposition 3.1 and Lemma 3.2 and using the same notation and conventions we immediately obtain: Corollary 3.3. The 3-manifold Y P bounds a rational homology ball if and only if either r = n 1 + 1 or the lens space L P ′ bounds a rational homology ball.

Proof. From the proof of Lemma 3.2 we know that if r = n 1 + 1 then the string associated to P ′ is (2, 2..., 2, 1) and therefore Y P ′ = S 3 which obviously bounds a rational homology ball. On the other hand, if r < n 1 + 1 then the string is of the form (a n1,1 , ..., a r,1 , 2, ..., 2, 1) and therefore Y P ′ = L P ′ where L P ′ is a lens space with associated string (a n1,1 , ..., a r,1 -1). We conclude using Proposition 3.1.

3.1. The List. The above analysis yields a complete list L of Seifert spaces with 3 exceptional fibers and two complementary legs bounding rational homology balls and proves Theorem 1.1. Every element in L is, with one of the two orientations, the boundary of a 4-manifold obtained by plumbing according to a graph -a1,2 -an 2 ,2 Proof. Since P has two complementary legs, we know that there is a ribbon move that changes the link ML P into the disjoint union K(p, q) ⊔ U , where U stands for the unknot and K(p, q) is the 2-bridge link associated to the lens space L P ′ in the statement of Lemma 3.2 [8, Lemma 3.3 and discussion before it]. If Y P is the boundary of a rational homology ball, then, by Proposition 3.1, L P ′ is also the boundary of a rational homology ball. [10, Theorem 1.2] tells us then that K(p, q) ⊂ S 3 is the boundary of a ribbon surface Σ ′ with χ(Σ ′ ) = 1. The statement follows immediately.

Conversely, if there exist a surface Σ and a ribbon immersion Σ S 3 , with ∂Σ = ML P and χ(Σ) = 1, then it is well known that Y P is the boundary of a rational homology ball [START_REF] Eisenbud | Three-dimensional link theory and invariants of plane curve singularities[END_REF].

In this last proposition we have shown that any 3 legged slice Montesinos knot with two complementary legs is actually the boundary of a ribbon disk. We have thus:

Corollary 3.5. The slice-ribbon conjecture holds true for all the Montesinos knots ML P with P a 3 legged canonical graph with 2 complementary legs.

noteworthy remarks on lattice embeddings

• Similar families of 3 legged graphs have been studied in [START_REF] Lecuona | On the slice-ribbon conjecture for Montesinos knots[END_REF]. In that article all graphs Γ such that Γ ⊆ Z Γ define Seifert spaces which bound rational homology balls and the corresponding Montesinos links are the boundary of ribbon surfaces with Euler characteristic 1. In other words, Donaldson obstruction is sufficient to yield a complete classification. The situation is completely analogous in the case of linear graphs (lens spaces) [START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF]. In contrast, not all the sets studied in this article, that is 3 legged graphs P with two complementary legs such that P ⊆ Z P , define Seifert spaces which bound rational homology balls. For example, the Seifert space associated to • In all the cases studied in [START_REF] Lecuona | On the slice-ribbon conjecture for Montesinos knots[END_REF][START_REF] Lisca | Lens spaces, rational balls and the ribbon conjecture[END_REF][START_REF] Greene | The slice-ribbon conjecture for 3-stranded pretzel knots[END_REF][START_REF]On the slice-ribbon conjecture for pretzel knots[END_REF] if a graph Γ yields a 3-manifold Y Γ which bounds a rational homology ball, then all the coefficients appearing in the embedding Γ ⊆ Z Γ can be chosen to be equal to ±1. This is no longer true in the family studied in this note. For example,

Theorem 1 . 1 .

 11 The Seifert fibered manifold Y = Y (b; (α 1 , β 1 ), (α 2 , β 2 ), (α 3 , β 3 )) with β2 α2 + β3 α3 = 1 is the boundary of a rational homology ball if and only if Y belongs to the list L .

-e 5 + 3e 6 -e 1 + e 3 -e 3 + e 4 + e 5 e 5 -e 4 e 1 + e 2 e 1

 63354121 e 2 is not the boundary of a rational homology ball, since the lens space L(36, 19) with associated string (2, 10, 2) does not belong to [10, Theorem 1.2].

2e 4 -e 3 e 3 -e 1 e 1 + e 2 e 1 -e 2 satisfies

 311212 |v 1,1 • e 4 | = 2 and it is the boundary of a rational homology ball, since the lens space L(4, 1) has this property.

  Montesinos links and ribbon surfaces. Given that there is a one-to-one correspondence between the set of Seifert spaces with at most 3 exceptional fibers and the Montesinos links which arise as branch sets, we can translate the analysis on the rational homology balls developed so far into the language of Montesinos links and ribbon surfaces. Theorem 1.2 in the introduction is an immediate consequence of the following proposition. Proposition 3.4. Let P be a canonical graph with two complementary legs. The Montesinos link ML P ⊂ S 3 is the boundary of a ribbon surface Σ with χ(Σ) = 1 if and only if the Seifert space Y P is the boundary of a rational homology ball.

						•	. . .	•
	-an 1 ,1		-ar,1	-2		-2
	•	. . .	•	•	. . .	•
						-a1,3	-an 3 ,3
						•	. . .
							reversal and
	duality, can be found in [8, Remark 3.2] which summarizes [10, Lemmas 7.1, 7.2
	and 7.3].
	3.2.

•

where the strings (a 1,2 , ..., a n2,2 ) and (a 1,3 , ..., a n3,3 ) are related to one another by Riemenschneider's point rule, the number of -2 in the left leg is arbitrary (including the possibility zero) and the string (a r,1 -1, ..., a n1,1 ), up to order