Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media

In this paper, we study the robust linearization of nonlinear poromechanics of unsaturated materials. The model of interest couples the Richards equation with linear elasticity equations, employing the equivalent pore pressure. In practice a monolithic solver is not always available, defining the requirement for a linearization scheme to allow the use of separate simulators, which is not met by the classical Newton method. We propose three different linearization schemes incorporating the fixed-stress splitting scheme, coupled with an L-scheme, Modified Picard and Newton linearization of the flow. All schemes allow the efficient and robust decoupling of mechanics and flow equations. In particular, the simplest scheme, the Fixed-Stress-L-scheme, employs solely constant diagonal stabilization, has low cost per iteration, and is very robust. Under mild, physical assumptions, it is theoretically shown to be a contraction. Due to possible break-down or slow convergence of all considered splitting schemes, Anderson acceleration is applied as post-processing. Based on a special case, we justify theoretically the general ability of the Anderson acceleration to effectively accelerate convergence and stabilize the underlying scheme, allowing even non-contractive fixed-point iterations to converge. To our knowledge, this is the first theoretical indication of this kind. Theoretical findings are confirmed by numerical results. In particular, Anderson acceleration has been demonstrated to be very effective for the considered Picard-type methods. Finally, the Fixed-Stress-Newton scheme combined with Anderson acceleration provides a robust linearization scheme, meeting the above criteria.


Introduction
The coupling of fluid flow and mechanical deformation in unsaturated porous media is relevant for many applications ranging from modeling rainfall-induced land subsidence or levee failure to understanding the swelling and drying-shrinkage of wooden or cement-based materials. Assuming linear elastic behavior, the process can be modeled by coupling the Richards equations with quasi-static linear elasticity equations, generalizing the classical Biot equations [1]. In this work, we consider utilize the equivalent pore pressure [2], which allows a thermodynamically stable formulation [3].
For the numerical simulation of large scale applications, the solution of linear problems is typically the computationally most expensive component. For nonlinear problems, the dominating cost is determined by both the chosen linearization scheme and solver technology. In particular, the choice of a linearization scheme defines the requirements for the solver technology. Commonly, Newton's method is the first choice linearization scheme. However, for the nonlinear Biot equations, as monolithic solver, Newton's method requires the solver technology to solve saddle point problems coupling mechanics and flow equations. Additionally, in practice, constitutive laws employed in the model might be not Lipschitz continuous [4]. Thus, the arising systems are ill-conditioned and require an advanced, monolithic simulator. As the latter might be not available, the goal of this work is to develop a linearization scheme, which is robust and allows the use of decoupled simulators for mechanics and flow equations. For this purpose, we adopt closely related concepts for the linear Biot equations and the Richards equations.
For the numerical solution of the linear Biot equations, splitting schemes are widely used; either as iterative solvers [5] or as preconditioners [6]. In particular, the fixed-stress splitting scheme has aroused much interest, being unconditionally stable in the sense of a von Neumann analysis [7] and a global contraction [8][9][10]. As iterative solver, the scheme has been extended in various ways; e.g., it can be rewritten to a parallel-in-time solver [11], a multi-scale version allowing separate grids for the mechanics and flow problem has been developed [12], and the concept has been extended to nonlinear multi-phase flow coupled with linear elasticity [3]. In the context of monolithic solvers, the scheme has been applied as preconditioner for Krylov subspace methods [13][14][15][16] and as smoother for multigrid methods [17]. All in all, the scheme defines a promising strategy to decouple mechanics and flow equations.
For the linearization of the Richards equation, the standard Newton method has to be used with care, since the Richards equation is a degenerate elliptic-parabolic equation, modeling saturated/unsaturated flow, and additionally material laws might be Hölder continuous. Various problem-specific alternatives have been developed in the literature. We want to point out two particular, simple linearization schemes; the L-scheme and the Modified Picard method. The L-scheme [18], employs diagonal stabilization for monotone, Lipschitz continuous nonlinearities. Global convergence has been rigorously proven for several porous media applications [19][20][21]; in particular also for the Richards equation [22]. The L-scheme can be also applied for Hölder continuous problems [21,23]. Furthermore, for the Richards equations it can be used to define a robust, linear domain decomposition method [24]. The L-scheme linearization has been coupled with the fixed-stress splitting scheme for nonlinear Biot's equations with linear coupling [25]. Less robust, but in some cases more efficient is the Modified Picard method [26], which employs the first order Taylor approximation for the saturation, still allowing for a Hölder continuous permeability.
In this paper, we combine the fixed-stress splitting scheme with the L-scheme, the Modified Picard method and Newton's method. The resulting schemes decouple and linearize simultaneously the mechanics and flow equations, utilizing only a single loop and allowing for separate simulators. We show theoretically linear convergence of the Fixed-Stress-L-scheme, assuming non-vanishing residual saturation, permeability and porosity. However, the theoretical convergence rate might deteriorate in unfavorable situations, leading to either slow convergence or even stagnation in practice. As remedy we apply Anderson acceleration.
Anderson acceleration, originally introduced by [27], in order to accelerate fixed point iterations in electronic structure computation, has been successfully applied in various other fields; in particular for the modified Picard iteration [28]. Reusing previous iterations to approximate directional derivatives it can be related to a multi-secant quasi-Newton method [29] and to preconditioned GMRES for linear problems [30]. For nonlinear problems, it naturally generalizes to a preconditioned nonlinear GMRES. Being a post-processing, it can be combined with splitting methods, still allowing separate simulators unlike preconditioned monolithic solvers. So far, theoretical results guarantee only convergence for contractive fixed-point iterations [31]. Furthermore, those results do not guarantee actual acceleration. Based on a special case, we justify theoretically the ability of the Anderson acceleration to effectively accelerate convergence of contractive fixed-point iterations and moreover stabilize non-contractive fixed-point iterations. Applying Anderson acceleration for diverging methods might allow convergence after all. To our knowledge, this is the first theoretical indication of this kind. Other stabilization techniques could be applied as adaptive step size control, adaptive time stepping or the combination of a Picard-type method with a Newton-type method, following ideas by [32]. These concepts have not been considered in the scope of this work.
We present numerical results confirming the theoretical findings of this work. Indeed, the Fixed-Stress-L-scheme is more robust than the modifications employing Newton's method and the Modified Picard method. Moreover, convergence of the Picard-type methods can be significantly accelerated by the Anderson acceleration. When applied to initially diverging methods, convergence can be reliably retained.
The main, new contributions of this work are: • We propose three linearization schemes incorporating the fixed-stress splitting scheme, coupled with an Lscheme, Modified Picard and Newton linearization of the flow. All schemes allow the efficient and robust decoupling of mechanics and flow equations. For the simplest scheme, the Fixed-Stress-L-scheme, we show theoretical convergence, assuming non-vanishing residual saturation, permeability and porosity, cf. Theorem 7.
• Based on a special case, we justify theoretically the general ability of the Anderson acceleration to effectively ac-celerate convergence and stabilize the underlying scheme, allowing even non-contractive fixed-point iterations to converge, cf. Section 7.3.
• The combination of the proposed linearization schemes and Anderson acceleration is demonstrated numerically to be robust and efficient. In particular, Anderson acceleration allows the schemes to converge in challenging situations. The Fixed-Stress-Newton method coupled with Anderson acceleration shows best performance among the splitting schemes, cf. Section 8.
The paper is organized as follows. In Section 2, the mathematical model for the nonlinear Biot equations is explained. In Section 3, a three-field discretization is introduced, employing linear Galerkin finite elements and mixed finite element for the mechanics and flow equations, respectively. In Section 4, we recall the monolithic Newton method and introduce three splitting schemes, simultaneously linearizing and decoupling the mechanics and flow equations. In Section 5, convergence is proved for the Fixed-Stress-L-scheme. In Section 6, Anderson acceleration is recalled, and in Section 7, the ability of the Anderson acceleration to effectively accelerate convergence and increase robustness is discussed theoretically. In Section 8, numerical results are presented, illustrating in particular the increase of robustness via Anderson acceleration. The work is closed with concluding remarks in Section 9.

Mathematical model -Nonlinear Biot's equations coupling Richards equation and linear elasticity
We consider a nonlinear extension of the classical, linear Biot equations modeling flow in deformable porous media under possibly both fully and partially saturated conditions. Further more we assume: (A1) The bulk material is linearly elastic and deforms solely under infinitesimal deformations.
(A2) There exists two fluid phases -one active and one passive phase (standard assumption for the Richards equation).
(A3) The active fluid phase is incompressible and corresponding fluxes are described by Darcy's law.
(A4) Mechanical inertia effects are negligible allowing to consider the quasi-static balance of linear momentum.
We model the medium at initial conditions by a reference configuration Ω ⊂ R d , d ∈ {2, 3}. Due to the limitation to infinitesimal deformations, the domain of the primary fields is approximated by Ω on the entire time interval of interest (0, T ), with final time T > 0. Finally, the governing equations describing coupled fluid flow and mechanical deformation of a porous medium with mechanical displacement u, fluid pressure p w and volumetric flux q w as primary variables are given by where φ denotes variable porosity, s w denotes fluid saturation, k w denotes fluid-dependent mobility, ρ w and ρ b denote fluid and bulk density, respectively, g is the gravitational acceleration, ε(u) and ∇ · u denote the linear strain and the volumetric deformation, respectively, µ and λ denote the Lamé parameters, α is the Biot coefficient and p E denotes the pore pressure. In the following, we comment briefly on the single components of the mathematical model and refer to [2] for a detailed derivation.
Eq. (1): For the fluid flow, an active and a passive fluid phase are assumed. In other words, the passive phase responds instantaneously to the active phase and therefore has a constant pressure. The behavior of the active fluid phase is governed by mass conservation, identical to volume conservation for an incompressible fluid. The volume is given by the product of porosity φ and saturation s w . The porosity changes linearly with volumetric deformation ∇ · u and pore pressure p E by where φ 0 , u 0 and p w,0 are the initial porosity, displacement and pressure, respectively, and α is the Biot coefficient and N is the Biot modulus. Eq. (4) is a byproduct of the thermodynamic derivation of the effective stress by Coussy [2]. Furthermore, the saturation s w = s w (p) is assumed to be described by a material law s w : R → (0, 1], satisfying s w (p) = 1, p ≥ 0, and having a negative inverse p c : (0, 1] → R + , satisfying s w (−p c (s)) = s, s ∈ (0, 1]. In the literature, the function p c is often referred to as capillary pressure. Eq. (2) The volumetric flux q w is assumed to be described by Darcy's law for multiphase flow. Here, the permeability scaled by viscosity is given by a material law k w = k w (s w ). In practice, the material laws can become Hölder continuous.
Eq. (3): The mechanical behavior is governed by balance of linear momentum under quasi-static conditions, combined with an effective stress formulation. Allowing only for small deformations, we employ the St.Venant Kirchhoff model for the effective stress, determining the poroelastic stress as σ por (u, p w ) = 2µε(u) + λ∇ · uI − αp E (p w ).
As pore pressure, we use the equivalent pore pressure [2] which takes into account interfacial effects. By construction it satisfies d p E = s w (p) d p. As body force we assume solely gravity, where for the sake of simplicity the bulk density ρ b is assumed to be constant. All in all, Eq. (3) acts as compatibility condition to be satisfied at each time.
∂Ω of the boundary of Ω and the outer normal n on ∂Ω, we assume boundary conditions and initial conditions All in all, the nonlinear Biot equations (1)-(3) couple nonlinearly the Richards equation and linear elasticity equations. In the fully saturated regime (p w ≥ 0), the model reduces locally to the classical, linear Biot equations for an incompressible fluid and compressible rock. We note, as long as the fluid saturation is not vanishing, the nonlinear Biot equations (1)-(3) are parabolic, unlike the degenerate elliptic-parabolic Richards equation, cf. Remark 3.

Finite element discretization
We discretize the Biot equations (1)-(3) in space and time by the finite element method and the implicit Euler method, respectively. More precisely, given a regular triangulation T h of the domain Ω, we employ linear, constant and lowest order Raviart-Thomas finite elements to approximate displacement, pressure and volumetric flux, respectively. For the sake of simplicity, we assume zero boundary conditions on ∂Ω = Γ m D = Γ f D . The corresponding discrete function spaces are then given by where P 0 and P 1 denote the spaces of scalar piecewise constant and piecewise (bi-)linear functions, respectively, and elements of V h are zero on the boundary. We note, that the chosen discretization is not stable with respect to the full range of material parameters [33], e.g., for very small permeability. However, the specific choice of the finite element spaces is not essential for the further discussion of the linearization. Furthermore, for the temporal discretization we employ a partition {t n } n of the time interval (0, T ) with (constant) time step size τ = t n − t n−1 > 0.
Then given initial data (p, u) 0 h ∈ W h × V h , at each time step n ≥ 1, the discrete problem reads: Given (p, q, u) n−1 2µ where s k w = s w (p k h ) and p k Here, ·, · denotes the standard L 2 (Ω) scalar product.

Monolithic and decoupled linearization schemes
In the following, we consider four linearization schemes. First, we apply the monolithic Newton method, being commonly the first choice when linearizing a nonlinear problem. Second, we propose a linearization scheme, which employs constant diagonal stabilization in order to linearize and decouple simultaneously the mechanics and flow equations. Furthermore, we introduce two modifications of the latter method, utilizing both decoupling and first order Taylor approximations. For direct comparison, we formulate all schemes in incremental form.
Monolithic schemes vs. iterative operator splitting schemes for saddle point problems. Biot's equations yield a saddle point problem, which in general are difficult to solve. Containing more information on coupling terms, a monolithic scheme for this purpose is per se more stable, whereas for iterative operator splitting schemes, stability is always an issue necessary to be checked. However, in contrast to robust splitting schemes, for a monolithic scheme a fully coupled simulator with advanced solver technology is required. For that, one possibility is to apply a splitting scheme as either an iterative solver or a preconditioner for a Krylov subspace method. The latter is more efficient and robust, cf. e.g. [13] for the linear Biot equations. In case only separate simulators are available, the concept of preconditioning the coupled problem cannot be applied in the same sense. But we note that acceleration techniques as Anderson acceleration can be applied as post-processing to the iterative splitting schemes, acting as preconditioned nonlinear GMRES solvers applied to the coupled problem, cf. Section 6.

Notation of residuals
For the incremental formulation of the linearization schemes, we introduce naturally defined residuals of the coupled problem (1)- (3). Given data (p, q, u) n−1 h ∈ W h × Z h × V h for time step n − 1, the residuals at time step n evaluated at some state (p, Shorter, given a sequence of approximations (p, q, u) n,i

Monolithic Newton's method
We apply the standard Newton method linearizing the coupled problem (6)-(8) in a monolithic fashion.
Scheme. The monolithic Newton method reads: For each time step, given the initial guess (p, q, u) n,0 h = (p, q, u) n−1 h , loop over the iterations i ∈ N until convergence is reached. Given data at the previous time step n −1 and iteration i−1, and set After convergence is reached at iteration N, set (p, q, u) n h = (p, q, u) n,N h .
Properties. Newton's method is known to be quadratically, locally convergent, which makes the method commonly the first choice linearization method. However, in general it is not robust and has the following drawbacks: • In order to ensure convergence, the time step size has to be chosen sufficiently small depending on the mesh size. Then the initial guess is sufficiently close to the unknown solution.
• The need for a good initial guess can be relaxed by using step size control, allowing a bigger time step size. Anderson acceleration applied as post-processing can be interpreted as such, for more details, cf. Section 6.
• Bounded derivatives of constitutive laws have to be available. In practice, nonlinearities employed in the model (1)-(3) are not necessarily Lipschitz continuous, e.g., the relative permeability for soils. In particular, in the transition from the partially to the fully saturated regime, the derivative of the relative permeability modeled by the van Genuchten model [4] can be unbounded. Consequently, the Jacobian might become ill-conditioned.
• Due to the coupled nature of the problem, the linear system resulting from the problem (9)-(11) has a saddle point structure and is therefore ill-conditioned. Hence, an advanced solver architecture is required. In the context of Biot's equations, the application of a fixed-stress type solver or preconditioner [13] yields remedy.

Fixed-Stress-L-scheme -A Picard-type simultaneous linearization and splitting
We propose a novel, robust linearization scheme for Eq. (6)- (8). It is essentially a simultaneous application of the L-scheme linearization for the Richards equation, cf. e.g. [22], and the fixed-stress splitting scheme for the linear Biot equations, cf. e.g. [7], both utilizing diagonal stabilization. In the following, we refer to the scheme as Fixed-Stress-L-scheme (FSL). As derived in Section 5, it can be interpreted as L-scheme linearization of the nonlinear Biot equations reduced to a pure pressure formulation or alternatively as nonlinear Gauss-Seidel-type solver, consisting of cheap iterations allowing separate, sophisticated simulators for the mechanical and flow subproblems. In Section 5, we show convergence of the Fixed-Stress-L-scheme under physical assumptions. Before defining the Fixed-Stress-Lscheme, we recall main ideas of both the L-scheme and the fixed-stress splitting scheme.
Main ideas of the L-scheme. The L-scheme is an inexact Newton's method, employing constant linearization for monotone and Lipschitz continuous terms. For remaining contributions Picard-type linearization is applied. Effectively, this approach is identical with applying a Picard iteration with diagonal stabilization. All in all, no explicit derivatives are required at the price that only linear convergence can be expected. Under mild conditions, this concept has been rigorously proven to be globally convergent for various porous media applications, e.g., [20][21][22]25]. Moreover, as pointed out by [22], the resulting linear problem is expected to be significantly better conditioned than the corresponding linear problem obtained by Newton's method.
Regarding the nonlinear Biot equations, the saturation s w = s w (p) is classically non-decreasing and Lipschitz continuous. Hence, assuming that φ i−1 ≥ 0 on Ω, the above criteria apply to the saturation contribution in Eq. (9). The approximation of the saturation at iteration i is then given by where L ∈ R + is a sufficiently large tuning parameter, usually set equal to the Lipschitz constant L s of s w . As coupling terms are not monotone, Picard-type linearization is applied to the remaining contributions of Eq. (6)- (8).
Main ideas of the fixed-stress splitting scheme. Considering the linear Biot equations, their linearization results in a saddle point problem, thus, requiring an advanced solver technology for efficient solution. For this purpose, physically motivated, robust, iterative splitting schemes are widely-used as, e.g., the fixed-stress splitting scheme, originally introduced by [5]. As it decouples mechanics and flow equations, separate simulators can be utilized for both subproblems, reducing the complexity to solving simpler, better conditioned problems. The robust decoupling is accomplished via sufficient diagonal stabilization, introducing a tuning parameter β FS . Its optimization with guaranteed convergence is a research question on its own [8,9,34]. Concepts can be also extended to multiphase flow coupled with linear elasticity [3].
Scheme. We observe that applying the monolithic L-scheme, as just explained to the nonlinear, discrete Biot equations (6)- (8), results in a linear problem equivalent with that for single phase flow in heterogeneous media, for which the fixed-stress splitting scheme is an attractive solver [9]. Both schemes are realized via diagonal stabilization. Anticipating the dynamics to be mainly governed by the flow problem, cf. Assumption (A4), and the mechanics problem to be much simpler, a simultaneous application of the L-scheme and the fixed-stress splitting scheme yields an attractive linearization scheme incorporating the decoupling of flow and mechanics equations. Written as iterative scheme in incremental form, the resulting Fixed-Stress-L-scheme reads: For each time step, given the initial guess (p, q, u) n,0 h = (p, q, u) n−1 h , loop over the iterations i ∈ N until convergence is reached. For each iteration i, perform two steps:

1.
Step: Set L = L s , the Lipschitz constant of s w , and and set

2.
Step: and set After convergence is reached at iteration N, set (p, q, u) n h = (p, q, u) n,N h . 7 Properties. The Fixed-Stress-L-scheme inherits its properties from the underlying methods. It does not require the evaluation of any derivatives, increasing the speed of the assembly process. It is very robust but guarantees only linear convergence, as shown in Section 5, cf. Theorem 7. Furthermore, the Fixed-Stress-L-scheme requires solely the solution of the mechanical and flow problem, allowing separate simulators. In particular, the overall method utilizes a single loop in contrast to the Newton's method combined with a fixed-stress splitting scheme as iterative solver.

Quasi-Newton modifications of the Fixed-Stress-L-scheme
The Fixed-Stress-L-scheme employs constant linearization for the fluid volume φs w with respect to fluid pressure, utilizing an upper bound for the Lipschitz constant. In many practical situations, this approach is quite pessimistic. Recalling the assumption that the flow problem dominates the dynamics of the system, we expect the simultaneous application of the fixed-stress splitting scheme and more sophisticated flow linearizations to be only slightly less robust than the Fixed-Stress-L-scheme. Independent of the flow linearization, diagonal stabilization is added by the splitting scheme increasing the robustness. In the following, based on the derivation of the Fixed-Stress-L-scheme in Section 5, cf. Remark 4, we couple simultaneously a modified Picard method [26] and Newton's method with the fixed-stress splitting scheme yielding the Fixed-Stress-Modified-Picard method and the Fixed-Stress-Newton method, respectively. The modified Picard method, in particular, is a widely-used linearization scheme for the Richards equation and hence rises interest for its use for the linearization of the discrete, nonlinear Biot equations (6)- (8).

Fixed-Stress-Modified-Picard method.
Applied to the Richards equations, the modified Picard method employs a first order Taylor approximation as linearization for the saturation and a Picard-type linearization for the possibly Hölder continuous permeability. By employing a first order approximation of the fluid volume φs w with respect to fluid pressure instead, and by coupling simultaneously with the fixed-stress splitting scheme, we obtain a linearization scheme for Eq. (6)- (8). For later reference, we denote the resulting scheme by Fixed-Stress-Modified-Picard-scheme. It is essentially identical with the Fixed-Stress-L-scheme but with modified first fixed-stress step (1. Step). We exchange Eq. (13)-(14) with Fixed-Stress-Newton method. In case the permeability is Lipschitz continuous, the simultaneous application of the fixed-stress splitting scheme and linearization of the flow equations via Newton's method yields an attractive linearization scheme for Eq. (6)- (8). For later reference, we denote the resulting scheme by Fixed-Stress-Newton method. It is essentially identical with the Fixed-Stress-L-scheme but with modified first fixed-stress step (1. Step). We exchange Eq. (13)- (14) with We note that the Fixed-Stress-Newton method is also closely related to applying a single fixed-stress iteration as inexact solver for the linear problem (9)-(11) arising from Newton's method.

L 2 (Ω)-type stopping criterion
For the numerical examples in Section 8, we employ a combination of an absolute and a relative L 2 (Ω)-type stopping criterion, closely related to the standard algebraic l 2 -type criterion. Given tolerances ε a , ε r ∈ R + , we denote an iteration as converged if it holds

Convergence theory for simultaneous linearization and splitting via the L-scheme
In the following, we show convergence of the Fixed-Stress-L-scheme (13)-(15) under mild, physical assumptions. For this purpose, we formulate the nonlinear discrete problem (6)-(8) as an algebraic problem, reduce the problem to a pure pressure problem by exact inversion and apply the L-scheme as linearization identical to the Fixed-Stress-L-scheme (13)- (15). Convergence follows then from an abstract convergence result. For simplicity, we assume vanishing initial data and a homogeneous and isotropic material.
Algebraic formulation of the nonlinear, discrete Biot equations. Given finite element bases for W h × Z h × V h , the nonlinear, discrete Biot equations (6)- (8) translate to the algebraic equations We omit the detailed definition of the finite element matrices and vectors used in Eq. (20)- (22), as they are assembled in a standard way. We comment solely on their origin and their properties relevant for further discussion.
• Let p ∈ R n p , q ∈ R n q , u ∈ R n u denote the algebraic pressure, volumetric flux and displacement coefficient vectors In the following, let ·, · denote the classical l 2 -vector scalar product • Let S pp : R n p → R n p ×n p denote a diagonal matrix with element-wise saturation s w on the diagonal, i.e., S pp (p) kk = s w (p k ) for p ∈ R n p , k ∈ {1, ..., n p }.
• Let D pu ∈ R n p ×n u and D pq ∈ R n p ×n q denote the matrices corresponding to the divergence operating on displacement and volumetric flux spaces, respectively, mapping into the pressure space. Let local mesh information be incorporated, analog to the mass matrix M pp .
• Let p E : R n p → R n p correspond to the element-wise equivalent pore pressure p E . For given p ∈ R n p , each component of p E is given by p E (p) k = p E (p k ), k ∈ {1, ..., n p }.
• Let φ 0 ∈ R n p denote the initial porosity vector incorporating local mesh information such that φ 0 + αD pu u + 1 N M pp p E (p) corresponds element-wise to the actual porosity of a deformed material.
• Let K −1 qq : R n p → R n q ×n q denote the volumetric flux mass matrix, weighted by the nonlinear permeability contribution k −1 w (s w ) in Darcy's law, and incorporating local mesh information.
• Let A uu ∈ R n u ×n u denote the stiffness matrix, corresponding to the linear elasticity equations, incorporating local mesh information.
• f p ∈ R n p , f q ∈ R n q and f u ∈ R n u incorporate solution independent contributions as volume effects and Neumann boundary conditions and data at the previous time step. Furthermore, let local mesh information be incorporated.
Compact formulation of the algebraic problem. By inverting exactly Eq. (21) and Eq. (22) with respect to q and u, respectively, and insert into Eq. (20), we obtain the equivalent, reduced problem for p By defining the reduced problem (23) can be written in compact form L-scheme linearization. We linearize the abstract problem (26) using the L-scheme, introducing a sequence {p i } i ⊂ R n p approximating the exact solution p ∈ R n p . Given a user-defined parameter L ∈ R + , we set L pp = LM pp . Then given an initial guess p 0 ∈ R n p , the scheme is defined as follows: Loop over the iterations i ∈ N until convergence is reached. At iteration i, given data p i−1 ∈ R n p , find p i ∈ R n p solving the linear problem Remark 2 (Equivalence to the Fixed-Stress-L-scheme). The L-scheme (27)  Lemma 1 (Convergence of the L-scheme). Assume (26) and (27) both have unique solutions p ∈ R n p and p i ∈ R n p , respectively. Furthermore, let the following assumptions be satisfied: is in some sense monotonically increasing and Lipschitz continuous.
(L2) There exist constants k m , k M ∈ R + satisfying k m q 2 pp for all p,p ∈ R n p , i.e. K is in some sense Lipschitz continuous. Here, the subordinate matrix norm · M qq ,∞ is defined by K M qq ,∞ = sup (L3) There exists a constant q ∞ ∈ R + satisfying M −1 qq f q + D ⊤ p ∞ ≤ q ∞ for the solution of problem (26), i.e., boundedness is satisfied.
If the parameter L and the time step size τ are chosen such that 2 The proof of Lemma 1 is given in Appendix A. The proof is essentially the same as for the Richards equation by [22] but formulated in a slightly more general framework. Assumption (L1)-(L2) are generalized versions of assumptions made in [22] due to the possible global dependence of each component of b = b(p) on p.
Consequence for the Fixed-Stress-L-scheme. By Remark 2, the Fixed-Stress-L-scheme (13)- (15) is equivalent with the L-scheme (27). Therefore, we check Assumption (L1)-(L3) of Lemma 1 particularly for Eq. (23) in order to analyze the Fixed-Stress-L-scheme. We make the following physical assumptions: Employing the properties of b, and making use of the specific choice of the equivalent pore pressure (5), the Jacobian of b is given by Hence, D b (p) = D b (p) ⊤ for all p ∈ P φ≥0 with eigenvalues greater or equal than zero. Consequently, the largest value for the Rayleigh quotient (28) is given by the largest eigenvalue of M −1/2 pp D b (p)M −1/2 pp maximized over p ∈ P φ≥0 . The porosity vector φ 0 is essentially scaled by M pp . Additionally, as shown by [16], D pu A −1 uu D ⊤ pu is norm equivalent with M pp . Hence, also D b (p) is norm equivalent with the standard mass matrix with mesh-independent bounds. Together with employing the assumptions, we see there exists a largest eigenvalue with the value given by the smallest eigenvalue l b of M −1/2 pp D b (p)M −1/2 pp minimized over p ∈ P φ≥0 . From above discussion it follows that l b ∈ R + is mesh-independent. All in all, the proposed thesis follows.
Proof. As b is invertible and D b is symmetric, using the Inverse Function theorem, it holds The result follows from Lemma 2.
Analogously, we obtain: Corollary 5. Let Assumption (F1)-(F2) be satisfied. Then for all p,p ∈ P φ≥0 , it holds Corollary 6. Let Assumption (F1)-(F3) be satisfied. Then, Assumption (L2) is satisfied, in the sense, there exists a constant L K ∈ R + , satisfying for all p,p ∈ P φ≥0 , Furthermore, there exist constants k m , k M ∈ R + , satisfying for all p ∈ R n p , q ∈ R n q k m q 2 Proof. As the underlying k w = k w (s w ) is Lipschitz continuous, together with a scaling argument, it follows, there exists a constantL K ∈ R + satisfying Furthermore, as s w = s w (p) is Lipschitz continuous, and S pp is a diagonal matrix, there exists a constant L s ∈ R + satisfying All in all, with Lemma 2 and Corollary 5, Eq. (30) follows with L K =L K L s l −2 b . Eq. (31) follows directly from Assumption (F3) together with a scaling argument.
All in all, under the assumptions of non-vanishing residual saturation, permeability, and porosity, we obtain convergence for the L-scheme (27), which translates directly to the Fixed-Stress-L-scheme (13)-(15). Theorem 7. Let Assumption (F1)-(F4) be satisfied. Let p ∈ R n p and p i ∈ R n p be the solutions of the nonlinear problem (23) and the L-scheme (27), respectively. Assume they are unique. Let the initial guess p 0 ∈ R n p satisfy B p ( p 0 − p M pp ) ⊂ P φ≥0 , where B p (r) ⊂ R n p denotes the sphere with center p and radius r > 0. Let L and τ be chosen

Acceleration and stabilization by Anderson acceleration
The Fixed-Stress-L-scheme is expected to be a linearly convergent fixed-point iteration with the convergence rate depending on a tuning parameter. Its considered modifications employ a less conservative choice for the tuning parameter with the risk of failing convergence. Consequently, we are concerned with two issues -slow convergence and robustness with respect to the tuning parameter.
All presented linearization schemes in Section 4 can be interpreted as fixed-point iterations x i = F P(x i−1 ) = x i−1 + ∆F P(x i−1 ), where x i denotes the algebraic vector associated with (p, q, u) n,i h and ∆F P(x i−1 ) is the actual, computed increment within the linearization scheme. For such in general, Anderson acceleration [27] has been demonstrated on several occasions to be a suitable method to accelerate convergence. Furthermore, due to its relation to preconditioned, nonlinear GMRES [30], we also expect Anderson acceleration to increase robustness with respect to the tuning parameter for the considered linearization schemes. Both properties are justified by theoretical considerations in Section 7 and demonstrated numerically in Section 8.
Scheme. The main idea of the Anderson acceleration applied to a fixed-point iteration is to utilize previous iterates and mix their contributions in order to obtain a new iterate. The method is applied as post-processing not interacting with the underlying fixed-point iteration. In the following, we denote AA(m) the Anderson acceleration reusing m + 1 previous iterations, such that AA(0) is identical to the original fixed point iteration. We can apply AA(m) to post-process the presented linearization schemes. In compact notation, the scheme reads: For the specific implementation, we follow Walker and Peng [30]. In particular, in Step 4, we solve an equivalent unconstrained minimization problem, which is better conditioned, relatively small and cheap. The main price to be paid is the additional storage of the vectors ∆F P(x i−m−1 ), ..., ∆F P(x i−1 ) and F P(x i−m−1 ), ..., F P(x i−1 ) .
Properties. As post-processing Anderson acceleration does not modify the character of the underlying method, i.e., a coupled or decoupled character remains unchanged. In particular, in contrast to classical preconditioning, no monolithic simulator is required. Hence, Anderson acceleration is an attractive method in order to accelerate splitting schemes.
In many practical applications, effective acceleration can be observed. Though, there is no general, theoretical guarantee for the Anderson acceleration to accelerate convergence of an underlying, convergent fixed-point iteration. Theoretically, even divergence is possible [30]. In the literature, so far, theoretical convergence results are solely known for contractive fixed-point iterations [31]. For nonlinear problems, AA(m) is locally r-linearly convergent with theoretical convergence rate not larger than the original contraction constant if the coefficients α remain bounded. Without assumptions on α, AA(1) converges globally, q-linearly in case the contraction constant is sufficiently small. Both results only guarantee the lack of deterioration but not acceleration.
For a special, linear case, in Section 7, we show global convergence and theoretical acceleration for a variant of AA(1), fortifying the potential of Anderson acceleration. In particular, Corollary 10 predicts the ability of the Anderson acceleration to increase robustness, allowing non-contractive fixed-point iterations to converge. This motivates to apply AA(m) also to accelerate possibly diverging Newton-like methods as the monolithic Newton method and the Fixed-Stress-Newton method with the risk of loosing potential, quadratic convergence.

Theoretical contraction and acceleration for the restarted Anderson acceleration
For a special linear case, we prove global convergence of a restarted version of the Anderson acceleration. In particular, convergence for non-contractive fixed-point iterations and effective acceleration for a class of contractive fixed-point iterations is shown.

Restarted Anderson acceleration
The original Anderson acceleration AA(m) constantly utilizes the full set of m previous iterates. By defining the depth m ⋆ i = {i − 1 mod m + 1, m} in the first step of Algorithm 1 and apart from that following the remaining steps, we define a restarted version AA ⋆ (m) of AA(m), closer related to GMRES(m). In words, in each iteration we update the set of considered iterates by the most current iterate. And in case, the number of iterates becomes m + 1, we flush the memory and restart filling it again. In particular, for m = 1, the algorithm reads: Given: x 0 for i=0,2,4,..., until convergence do From [31], it follows directly, that for F P, a linear contraction, AA ⋆ (1) converges globally with convergence rate at most equal the contraction constant of F P. In the following, we extend the result to a special class of non-contractions.

Convergence result
For the convergence results, cf. Lemma 8 and Corollary 9, 10, we make the following assumptions: (C2) A is symmetric, and hence, A is orthogonally diagonalizable and there exists an orthogonal basis of eigenvectors {v j } j and a corresponding set of eigenvalues {λ j } j satisfying Av j = λ j v j .
(C3) There exists a unique x ⋆ such that F P(x ⋆ ) = x ⋆ , i.e., I − A is invertible.
(C4) The initial iterate x 0 is chosen such that the initial error To avoid a trivial case, we assume λ 1 , λ 2 0.
Then we are able to relate the errors between iterations of AA ⋆ (1), allowing to prove further convergence and acceleration results, cf. Corollary 10 and Corollary 9. All in all, the proof employs solely elementary calculations. However, as we are not aware of a general result of same type in the literature, we present the proof.
Proof. First, an iteration-dependent error propagation matrix is derived, and second, an upper bound for its spectral radius is computed. For this purpose, let us ignore for a moment Assumption (C4).

Iterative error propagation.
As we intend to relate e i+4 with e i , we explicitly write out the first four iterates and the corresponding errors. Given x i , by using b = x ⋆ − Ax ⋆ and x i − x i+1 = e i − e i+1 , we obtain By plugging all together, we obtain It suffices to bound the largest eigenvalue of the error propagation matrix A(A + α (i+3) (I − A))A(A + α (i+1) (I − A)). From Assumption (C2) it follows that {v i } i defines an orthogonal basis of eigenvectors for the error propagation matrix with corresponding eigenvalues {λ j } j defined bỹ Explicit definition of α (i+1) and α (i+3) . The minimization problem in Algorithm 2 can be solved explicitly, by solving adequate normal equations. It follows, that .
After employing simple arithmetics and using Eq. (32), we obtain

Discussion
We make the following comments: • The convergence result in Corollary 10 deals only with positive definite matrices. In Fig. 1b, eigenvalue pairs (λ 1 , λ 2 ) ∈ R × R are displayed satisfying r(λ 1 , λ 2 ) < 1 and therefore guaranteeing AA ⋆ (1) to converge. In particular, AA ⋆ (1) converges also for matrices with two eigenvalues larger than 1 with relatively close distance to each other.
• In practice, we do not experience AA ⋆ (1) or AA(1) to fail as long as Assumption (C4) is valid and |λ 1 | < 1 or |λ 2 | < 1. This observation extends also to arbitrarily large decompositions of e 0 as long as at most one eigenvalue of A satisfies |λ j | > 1. Based on similar observations, we state the following claim: If |λ j | > 1 for exactly m eigenvalues {λ j } j , then AA(m) converges for arbitrary e 0 . We note that the worst case approach used to in order to prove Lemma 8 cannot be applied to prove the general claim. It can be verified numerically that in general the eigenvalues of the error propagation matrix (41) can be larger than 1 even if ρ(A) < 1.
• From Fig. 1b, it follows, the closer the eigenvalues to 1, the slower the convergence of AA ⋆ (1). This is consistent with the interpretation of Anderson acceleration as secant method. The Richardson iteration does only damp slowly directions corresponding to eigenvalues close to 1. Hence, a directional derivative in these directions cannot be approximated well purely based on the iterations of fixed-point iterations. Quite contrary to directions corresponding to small or large eigenvalues relative to 1.
• The theoretical convergence result has been obtained from a worst case analysis. Practical convergence rates might be lower than predicted, depending on the weights of the initial error.
• Convergence of AA(m) is not guaranteed to be monotone, when applied for non-contractive fixed-point iterations.

Numerical results -Performance study
In this section, we will show three numerical examples, with increasing complexity, comparing the linearization schemes, presented in Section 4, coupled with Anderson acceleration. In particular, we confirm numerically the parabolic character of the nonlinear Biot equations, cf. Remark 3, the convergence result for the Fixed-Stress-Lscheme, cf. Theorem 7, as well as the acceleration and stabilization properties of the Anderson acceleration, cf. Section 7. All numerical results have been obtained using the software environment DUNE [36][37][38].

Test case I -Injection in a 2D homogeneous medium with Lipschitz continuous constitutive laws
We consider a two-dimensional, homogeneous, unsaturated porous medium (−1, 1) × (0, 1), in which a fluid is injected at the top, cf. Figure 2. Due to the symmetry of the problem, we consider only the right half Ω = (0, 1) ×(0, 1), discretized by 50 × 50 regular quadrilaterals. As initial condition, we choose a constant displacement and pressure field with u(0) = 0 and p w (0) = p 0 , satisfying the stationary version of the continuous problem (1)- (3). In order to avoid inconsistent initial data, we ramp the injection at the top with inflow rate q inflow (t) = q ⋆ × min {t 2 , 1.0} for given q ⋆ ∈ R. Apart from the inflow at the top, we consider no flow at the remaining boundaries, no normal displacement at left, right and bottom boundary and no stress on the top. The boundary conditions are displayed in Figure 2. Physical and numerical parameters. For the constitutive laws, governing saturation and permeability, we use the van Genuchten-Mualem model [4], defining where a vG and n vG are model parameters associated to the inverse of the air suction value and pore size distribution, respectively, k abs is the intrinsic absolute permeability and µ w is the dynamic fluid viscosity. Values chosen for model parameters and numerical parameters are displayed in Table 1. The parameters have been chosen such that the initial saturation is s w,0 = 0.4 in Ω and a region of full saturation (s w = 1) is developed after seven time steps. Furthermore, the constitutive laws for saturation and permeability are Lipschitz continuous (with L s = 0.12). We consider three different values for the Biot coefficient, controlling whether the Richards equation or the nonlinear coupling terms determine the character of the numerical difficulties. The simulation result for strong coupling (α = 1.0) at final time t = 1 is illustrated exemplarily in Figure 3.   Performance of linearization schemes. We consider the four linearization schemes introduced in Section 4, coupled with Anderson acceleration as post-processing. Abbreviations used in this section are introduced in Table 2.
We use the average number of iterations per time step as measure for performance, cf. Table 3. In particular, we disregard the use of CPU time as performance measure due to a not finely-tuned implementation. We just note, that a single iteration of a splitting method is significantly faster than a single monolithic Newton iteration.
First of all, all plain linearization schemes (AA(0)) succeed to converge for all three coupling strengths. This is consistent with Remark 3, demonstrating that the nonlinear Biot equations do not adopt the degeneracy of the Richards equation and remain parabolic in a fully saturated regime. Not surprisingly, the monolithic Newton method requires fewest iterations. So at first impression, it seems to be the preferred method. However, as stressed above, an advanced monolithic solver or an fixed-stress type iterative solver are required for efficient solution independent of the coupling strength, i.e., additional costs are hidden. On the other hand, the remaining linearization schemes allow separate simulators from the beginning. As the Fixed-Stress-L-scheme does not utilize an exact evaluation of derivatives, the Fixed-Stress-Newton method and the Fixed-Stress-Modified-Picard perform better for all three coupling strengths. Solely the performance of the Fixed-Stress-Newton method shows weak dependence on the coupling strength. The remaining methods show improved convergence behavior for increasing coupling strength, due to the decreasing numerical complexity of the problem itself following from Remark 3.   When applying Anderson acceleration, we observe that Anderson acceleration slows down the convergence of the monolithic Newton method, which is consistent with considerations in Section 6. In contrast, Anderson acceleration speeds up significantly the convergence of the Picard-type methods (Fixed-Stress-L-scheme and the variation FSL/2, and Fixed-Stress-Modified-Picard). Largest acceleration effect can be seen for largest considered depth. For the Fixed-Stress-Newton method, the effect of Anderson acceleration depends on the numerical character of the problem. This is due to the fact, that for weak coupling, the method is essentially identical with Newton's method.
Regarding the Fixed-Stress-L-scheme, according to Theorem 7, optimally the diagonal stabilization parameter has to be chosen as small as possible. However, smaller values do not necessarily lead to faster convergence, as can be observed by comparing the plain Fixed-Stress-L-scheme and the plain FSL/2-scheme. Yet when utilizing Anderson acceleration, robustness with respect to the tuning parameter is increased, and eventually the FSL/2-scheme converges faster than the Fixed-Stress-L-scheme. In particular, it performs as good as the Fixed-Stress-Modified-Picard method.
All in all, the theory has been confirmed. The Fixed-Stress-L-scheme converges despite the simple linearization approach and Anderson acceleration is able to accelerate Picard-type schemes. Moreover, the latter has been shown to stabilize the Fixed-Stress-L-scheme, allowing to choose a small tuning parameter leading to improved convergence behavior. Considering the cost per iteration, despite some additional iterations, we finally recommend the use of the Fixed-Stress-Newton method with Anderson acceleration with low depth. It is cheap and allows separate simulators. For strongly coupled problems or in the absence of exact derivatives, the Fixed-Stress-L-scheme with small tuning parameter is an attractive alternative to the Fixed-Stress-Newton method.

Test case II -Injection in 2D homogeneous medium with Hölder continuous permeability
In the following, we reveal the limitations of the considered linearization schemes. Moreover, we demonstrate the stabilization property of Anderson acceleration, allowing non-convergent methods to converge. For this purpose, we repeat test case I with modified physical parameters. In particular, we choose the saturation to be Lipschitz continuous with same Lipschitz constant as in test case I. In contrast, the permeability is chosen to be only Hölder continuous. Hence, the derivative becomes unbounded in the transition between partial and full saturation, causing potential trouble for the Newton-type methods. Again, we choose the initial pressure and the maximal inflow rate such that s w,0 = 0.4 and a region of full saturation (s w = 1) is developed after seven time steps. The simulation result for strong coupling at final time t = 1 is illustrated in Figure 4.
As mentioned in Remark 5, due to lack of regularity for the permeability, each of the considered methods faces difficulties. For Newton-type methods (Newton, Fixed-Stress-Newton), the derivative of the permeability is evaluated, which might be unbounded. Effectively, for the Fixed-Stress-L-scheme, this also means that L K → ∞ or in practice L K becomes very large. Hence, by Theorem 7, the time step size has to be chosen sufficiently small or possibly L has to be chosen larger to guarantee convergence. We note that for chosen initial saturation the permeability is significantly lower than for test case I. Consequently, the theoretical convergence rate for the plain Fixed-Stress-L-scheme (13)- (15) deteriorates. Due to round off errors stagnation is possible.  Performance of linearization schemes. The average number of iterations per time step is presented in Table 4. In contrast to test case I, not all plain linearization schemes (AA(0)) converge. For weak coupling, all Newton-like methods (Newton, Fixed-Stress-Newton) diverge with the Fixed-Stress-Newton method being slightly more robust due to added fixed-stress stabilization. The Fixed-Stress-L-scheme stagnates and shows to be slightly more robustness than the Newton-type methods. The Fixed-Stress-Modified-Picard method is least robust and stagnates already after two time steps. For strong coupling, all methods converge, which is consistent with Remark 3. If convergent, the schemes sorted by required number of iterations are the monolithic Newton method, the Fixed-Stress-Newton method, the Fixed-Stress-Modified-Picard and the Fixed-Stress-L-scheme, meeting our expectations. By utilizing Anderson acceleration, convergence can be observed for all coupling strengths and all linearization schemes besides Newton's method for α = 1. In particular, all previously failing schemes converge. This confirms the possible increase of robustness by Anderson acceleration, postulated in Section 7. Similar observations as before are made for the splitting schemes under Anderson acceleration. All in all, the theory has been confirmed.
As before, for increasing depth, the performance of Newton's method deteriorates. For strong coupling stagnation is observed. For weak coupling, for several time steps practical stagnation is observed with eventual convergence after a very large number of iterations. This is consistent with the fact that Anderson acceleration can also lead to divergence for increasing depth [30]. Hence, Anderson acceleration has to be applied carefully for the monolithic Newton method.
Motivated by test case I, we apply the Fixed-Stress-L-scheme with a decreased tuning parameter. For this test case, the plain FSL/2-scheme fails for all coupling strengths. As the FSL/2-scheme is a priori less robust as the Fixed-Stress-L-scheme, this has been expected. Utilizing Anderson acceleration, the FSL/2-scheme eventually converges. In particular, convergence is always faster than for the corresponding Fixed-Stress-L-scheme. This again demonstrates the ability of the Anderson acceleration to increase robustness and to relax assumptions for practical convergence.
According to the theory for the Fixed-Stress-L-scheme, a larger tuning parameter or a lower time step size could enable convergence, e.g., for α = 0.1, AA(0). However, we do not consider those strategies here, as they lead to worse convergence rates and utilizing Anderson acceleration should be anyhow preferred.
Concerning the best splitting method, we again recommend the use of the Fixed-Stress-Newton method combined with Anderson acceleration with low depth. It is cheap, robust and allows separate simulators.

Test case III -Unsteady seepage flow through a 2D homogeneous levee
We consider unsteady seepage flow through a simple, two-dimensional, homogeneous levee, enforced by a flood. The levee consists of a lower and upper part (lower 5 [m] and upper 10 [m], respectively), cf. Figure 5. Initially, the water table lies at the interface between lower and upper part. The initial fluid pressure is a hydrostatic pressure with p = 0 at the water table. The reference configuration, defined by the domain, is initially already consolidated under the influence of gravity. As u is the deviation of the reference configuration, effectively, no gravity is applied in the mechanics equation, but only in the flow equation. on the left, a hydrostatic pressure boundary condition is applied. On the right side, we apply approximate seepage face boundary conditions, based on the previous time step; i.e., given a fully saturated cell at the previous time step, a pressure boundary condition p = 0 is applied on corresponding boundary for the next time step, otherwise a no-flow boundary condition is applied for the volumetric flux. On the remaining boundary, no-flow boundary conditions are applied for all time. For the mechanics, no displacement in normal direction is assumed on the boundary of the lower part of the levee. On the boundary of the upper part and the interface, zero effective stress is applied. The boundary conditions are visualized in Figure 5.
Physical and numerical parameters. The domain is discretized by a regular, unstructured, simplicial mesh with approximately 67,000 elements and 201,000 nodes. Compared to the previous test cases, we employ more realistic material parameters. Values chosen for model parameters and numerical parameters are displayed in Table 1. We note, the resulting permeability is only Hölder continuous. The saturation history and deformation at four times is displayed in Figure 6. We observe steep saturation gradients during the flooding. Furthermore, both consolidation and swelling can be observed. All in all, the levee is pushed to the right.
Performance of linearization schemes. We consider the same linearization schemes as in the previous test cases, all but FSL, i.e., the Fixed-Stress-L-scheme with L = L s + β FS . Based on previous observations, we expect FSL/2 coupled with Anderson acceleration to be more efficient than FSL. The average number of iterations per time step is presented in Table 5. First of all, we observe that all plain linearization schemes fail in the same phase of the simulation (after around 50 time steps). The reason for that lies mainly in the steep saturation gradients. As before, Anderson acceleration can yield remedy. However, for this test case, the simple combination of Newton's method and Anderson acceleration does not converge for any considered depth. Newton's method combined with Anderson acceleration is still not convergent for AA (1). For increasing depth the robustness decreases again, which is consistent with observations from the previous test cases. For the remaining linearization schemes convergence can be obtained. In particular, the Fixed-Stress-Newton method combined with AA(1) converges with the least amount of iterations. The Picard-type methods are slower, but show again more robustness with respect to increasing depth, whereas the Fixed-Stress-Newton method diverges eventually for m = 10. Here, the Picard type methods require at least depth m = 3 for successful convergence. After all, we conclude that the diagonal stabilization is essential for the success of the linearization schemes. The stabilization is added via both the fixed-stress splitting scheme and the L-scheme. Consequently, we expect also the monolithic Newton method to be convergent when adding sufficient diagonal stabilization.  Table 5: Performance for test case III. Average number of (nonlinear) iterations per time step for Newton's method, the Fixed-Stress-Newton method, the Fixed-Stress-Modified-Picard method and the Fixed-Stress-L-scheme; both plain and coupled with Anderson acceleration for different depths (m = 1, 3, 5, 10). Minimal numbers per linearization type are in bold. Failing linearization due to stagnation at time step n is marked by → [n]. Failing linearization due to divergence at time step n is marked by ր [n].

Concluding remarks
In this paper, we have proposed three different linearization schemes for nonlinear poromechanics of unsaturated materials. All schemes incorporate the fixed-stress splitting scheme and allow the efficient and robust decoupling of mechanics and flow equations. In particular, the simplest scheme, the Fixed-Stress-L-scheme, employs solely constant diagonal stabilization. It has been derived as L-scheme linearization of the Biot equations reduced to a pure pressure formulation. Under mild, physical assumptions, also needed for the mathematical model to be valid, it has been rigorously shown to be a contraction. This also has been verified numerically. Exploiting the derivation of the Fixed-Stress-L-scheme allows modifications including first order Taylor approximations. In this way, we have introduced the Fixed-Stress-Modified-Picard and the Fixed-Stress-Newton method.
The derivation of the Fixed-Stress-L-scheme provides two particular side products. First, it reveals the close relation of the L-scheme and the fixed-stress splitting scheme. Second, the nonlinear Biot equations can be shown to be parabolic in the pressure variable. This holds in particular in the fully saturated regime unlike for the Richards equation.

23
The theoretical convergence rate of the Fixed-Stress-L-scheme might deteriorate for unfavorable situations, leading to slow convergence or even stagnation in practice. Similarly, the Fixed-Stress-Modified-Picard and Fixed-Stress-Newton methods are prone to diverge for Hölder continuous nonlinearities. In order to accelerate or retain convergence, we apply Anderson acceleration, which is a post-processing, maintaining the decoupled character of the underlying splitting methods. The general increase of robustness and acceleration of convergence via the Anderson acceleration has been justified theoretically considering a special linear case. To our knowledge, this is the first theoretical indication of this kind, considering non-contractive fixed-point iterations.
In practice, Anderson acceleration has shown to be very effective for the considered Picard-type methods, confirming the theoretical considerations. After all, we recommend the combination of the Fixed-Stress-Newton method and the Anderson acceleration, being very robust even for Hölder continuous nonlinearities. In case analytical derivatives are not available, we recommend the combination of the Fixed-Stress-L-scheme with a decreased tuning parameter and Anderson acceleration. Without Anderson acceleration, convergence might not be guaranteed. Including it, does not only retain but it also significantly accelerates convergence. This is interesting, as the optimal tuning parameter is not necessarily known a priori and can be more safely approached under the use of Anderson acceleration.
As outlook, with focus on large scale applications, the performance of the linearization schemes should be analyzed under the use of parallel, iterative solvers; in particular, as due to added stabilization, the arising linear systems are expected to be better conditioned than for the monolithic Newton method. Additionally, Anderson acceleration should be further studied in the context of possibly non-contractive fixed point iterations. Examples are (i) the linearization of degenerate problems including Hölder continuities, which are known to be difficult to solve [23], and (ii) numerical schemes employing a tuning parameter. Based on the numerical results in this paper, the approach seems very promising.