Iterative splitting schemes for a soft material poromechanics model

We address numerical solvers for a poromechanics model particularly adapted for soft materials, as it generally respects thermodynamics principles and energy balance. Considering the multi-physics nature of the problem, which involves solid and fluid species, interacting on the basis of mass balance and momentum conservation, we decide to adopt a solution strategy of the discrete problem based on iterative splitting schemes. As the model is similar (but not equivalent to) the Biot poromechanics problem, we follow the abundant literature for solvers of the latter equations, developing two approaches that resemble the well known undrained and fixed-stress splits for the Biot model. A thorough convergence analysis of the proposed schemes is performed. In particular, the undrained-like split is developed and analyzed in the framework of generalized gradient flows, whereas the fixed-stress-like split is understood as block-diagonal $L^2$-type stabilization and analyzed by means of a relative stability analysis. In addition, the application of Anderson acceleration is suggested, improving the robustness of the split schemes. Finally, we test these methods on different benchmark tests, and we also compare their performance with respect to a monolithic approach. Together with the theoretical analysis, the numerical examples provide guidelines to appropriately choose what split scheme shall be used to address realistic applications of the soft material poromechanics model.


Introduction
Poromechanics addresses the behavior of fluid-saturated permeable porous materials, and in particular the interaction of their mechanical deformation and the fluid flow. Since its origin in the context of civil engineering [1][2][3][4], most commonly known as Biot's theory of poroelasticity, it has been used for countless applications of societal and industrial relevance, e.g., in reservoir geomechanics, hydrology and soil mechanics, and material sciences (see the review [5] and the references therein). More recently, it has also captured the attention of researchers interested in the behavior of highly deformable, soft biological tissues [6][7][8]; a prominent example is the perfusion of the heart [9][10][11][12].
As the classical theory of poroelasticity and resulting models were originally developed for civil applications, they are in general inadequate for biomechanics. This difficulty can be more evidently appreciated when considering soft tissues undergoing large deformations and perfusion with potentially moderate flow rates [13][14][15][16]. Ultimately, this has called for more general formulations obeying the fundamental principles of continuum mechanics and thermodynamics [17], which renders these models applicable to a broader range of scenarios.
Among various advances, we particularly highlight the development of a general, thermodynamically consistent poromechanics model by Chapelle and Moireau [18], which also serves as basis for this work. The model is based on a thermodynamic derivation combined with thermodynamically consistent constitutive laws. It couples the balance of linear momentum for the solid and fluid phases including the viscous dissipation governed by the interaction of both phases due to friction, as well as the conservation of mass. Most importantly, in contrast to the classical quasi-static Biot consolidation model, the aforementioned model satisfies an energy-dissipation identity, predicting the dissipation of the combination of the kinetic and Helmholtz free energy. A further difference between the two approaches is that the former considers the absolute fluid velocity instead of the relative one.
The analysis of the well-posedness, stability and numerical approximation of this class of poromechanics models is still largely open. Among recent advances, we highlight the development and analysis of an implicit time discretization preserving the dissipation-energy identity at the discrete level [18]; an energy-preserving implicit-explicit time discretization incorporating a (non-iterative) operator splitting, decoupling solid and flow computations [19]; an energy-stable space and time discretization for a linearized model with focus on quasi-incompressible solids [20]; and finally, a space and time discretization for the same linearized model, exploiting a generalized saddle point structure and ultimately suggesting the use of Taylor-Hood type finite elements [21].
Motivated by the success of block-partitioned solvers for the related, classical quasi-static Biot equations, the main objective of this work is to develop and analyze for the first time iterative coupling strategies for the general, thermodynamically consistent poromechanics model proposed in [18]. Similar to previous theoretical works in this context, see for example [20,21], a linearized model is considered for the numerical analysis.
In general, solvers decoupling different physics allow the employment of methods tailored to the separate sub-problems, as flow and elasticity. However, a sequential-implicit solution requires iterating until convergence at each time step. In contrast, fully-implicit approaches, solving the fully-coupled problem at once, yield unconditional stability but require advanced and efficient preconditioners. Here it is worth to mention that robust iterative coupling strategies can effectively guide the design of scalable preconditioners for the monolithic solution by Krylov subspace methods.
For robust iterative coupling, in general, a problem-specific strategy is required; yet, we can learn from the well-studied, related Biot equations. For the latter, solvers based upon a sequential-implicit solution of the flow and mechanics sub-problem have been studied since over two decades [22]. The most popular iterative schemes are the undrained split [23] and the fixed-stress split [22], both relying on additional stabilization to one of the sub-problems. Due to suitable choices of stabilization, both have been shown to be unconditionally stable [24][25][26] with theoretical convergence rates depending on stabilization and model parameters, but independent on mesh properties; inf-sup stability of the discretization even allows robust convergence in the fluid-incompressible and quasi-impermeable regime [27]. Moreover, the fixed-stress split has been successfully generalized to several complex extensions of the quasi-static Biot equations. In view of biomedical applications, we emphasize work on large deformations [28]. For optimal performance of the iterative solvers, the choice of the stabilization is well-known to be vital. This choice does depend on several factors [27] as problem parameters, but also boundary conditions and geometry, which are difficult to quantify. To alleviate this, it has been shown in [29] that Anderson acceleration [30] greatly relaxes the requirement of optimal stabilization. Furthermore, utilizing the fact that stabilized split schemes are equivalent to a preconditioned Richardson iteration [25], they provide a basis to design efficient block preconditioners for the fully-implicit approach [31], next to alternative efficient preconditioners [32][33][34][35][36]. In this context, the need for optimal stabilization is similarly relaxed. Hence, after all, stabilization parameters derived in theoretical analyses offer a practical choice.
In this work, we develop and analyze splitting schemes for a linearization of the general poromechanics model [18], previously introduced and analyzed in [20,21]. This (linearized) model resembles Biot's equations, but presents fundamental differences, most importantly, new terms in the momentum equations of the fluid and solid phases due to inertia, and a structurally different saddle-point structure, compared to a double saddle point structure of the Biot equations. Still, iterative coupling concepts can be adapted to the new setting. Ultimately, we present schemes similar to the undrained split and the fixed-stress split. In particular, the undrained-like split is developed and analyzed in the framework of generalized gradient flows and alternating minimization following [37], whereas the fixed-stress-like split is understood as block-diagonal L 2 -type stabilization and analyzed by means of a relative stability analysis. In practice, additional application of Anderson acceleration is suggested, motivated by associated works in the literature and the here presented numerical examples.
This work is structured as follows. In Section 2, we present the general model of interest and its linearized version. In Section 3 and Section 4, we present respectively the alternating minimization split and the diagonally L 2 -stabilized split. The convergence of both schemes is analyzed in Section 5. In Section 6, an extensive numerical study is presented which validates the theoretical results. Finally, we close with concluding remarks in Section 7.

The thermodynamically consistent poromechanics model
The purpose of this work is to develop efficient solution strategies for the linearized and discretized version of the thermodynamically consistent poromechanics model originally developed by Chapelle and Moireau in [18], further described below. Two main steps are essential to reach this objective. One is the discretization of the equations (in this work we consider finite difference schemes in time and finite elements for the space discretization) and the other is the linearization of the model through a Newton-Raphson method. It is natural to operate the linearization on the discrete version of the problem, obtaining a discrete tangent problem to which the solution strategies proposed later on will be applied. This can be named the discretize then linearize strategy.
We remark that in the definition of the tangent problem the shape derivatives are neglected, namely the physical domain Ω t is identified with the reference one Ω 0 . As in this case the tangent problem depends only on the Fréchét derivatives of the mathematical operators that govern the nonlinear problem, the discrete tangent problem obtained by means of the discretize then linearize approach is equivalent to the one that would be derived from the linearize then discretize strategy. The latter strategy corresponds to address the linearization of the continuous problem first, giving rise to a fully continuous tangent problem. Then, we address the numerical discretization of such problem and we develop the numerical solvers for it, based on the splitting into several sub-problems. We choose to follow the latter approach, because it is much simpler as it allows us to work with the strong formulation of the equations.

The general model for finite deformations
The model assumes that two phases, a fluid (f ) and a solid (s), coexist at each point of the domain of interest. Let us denote by φ the volume fraction of the fluid. We use Lagrangian (reference) and Eulerian (physical) coordinate frames, denoting by Ω 0 the domain in the Lagrangian frame and by Ω t the same domain in the deformed configuration. In the same way, we denote with the subindex (0) the operators defined in the Lagrangian frame. For example, given the displacement field in the Lagrangian frame, namely u s (x 0 , t) such that x = x 0 + u s (x 0 , t) for any x ∈ Ω t , x 0 ∈ Ω 0 , the deformation gradient tensor is F = I + ∇ 0 u s and its determinant is J = detF. We also introduce the symbol J s = J(1 − φ). One of the primary variables of the model is the added mass m = ρ f (Jφ − φ 0 ) that is the fluid mass added to the system due to pore deformation. To define the constitutive laws of the model, we introduce Ψ(F, J s ) which is a suitable free energy of the solid.
In view of the linearization of the problem, we formulate the equations on the following abstract form: Find the velocity of the solid phase v s , the velocity of the fluid phase v f and the added mass (per unit volume) m, such that where the operators S(·), F(·), M(·) correspond to the momentum conservation in the solid and fluid phases, and the mass balance, respectively. More precisely, referring to the strong formulation of the model presented in [20], the operators S(·), F(·), M(·) correspond to the following sub-problems: complemented by the following constitutive laws Given v s , m and f in Ω t , find v f in Ω t such that Here, ρ s and ρ f constitute (spatially and temporally) constant densities of the solid and fluid phases, respectively, k f denotes the fluid mobility (absolute permeability divided by the fluid viscosity). Potentially, φ 0 , k f , and κ s are spatially varying, and the source f is varying in space and time.
The problem must be complemented by boundary and initial conditions. For the boundary constraints many options are possible, as discussed for example in [20]. For the sake of simplicity, we present here only one of the possible variants. Let us split the whole boundary ∂Ω t into two distinct non-intersecting parts, Γ D t and Γ N t , where we enforce Dirichlet and Neumann type conditions, respectively. Let v D s , v D f , t, be assigned velocities and traction for boundary conditions, and let v 0 s , v 0 f be the assigned initial values, under the assumption that Ω t = Ω 0 at t = 0, We define the boundary and initial conditions as follows, in Ω 0 × {0}.

Derivation of the tangent problem
Using the previous abstract form of the problem, we formally derive the tangent problem. To this purpose, we denote by D u A the derivative of a generic operator A with respect to the field u. We point out that such derivative should account for the classical Fréchét derivative of the operator, combined with the shape derivatives due to deformations of the domain. The central hypothesis in the definition of the tangent problem is that we neglect the shape derivatives, limiting ourselves to account for the Fréchét ones. In other words, we identify the physical domain, Ω t , with the reference one, Ω 0 (and for simplicity we drop the subindices 0, t, denoting Ω t and Ω 0 both by Ω). In this setting, the nonlinear problem is approximated, at the point v s , v f , m, by the following linear problem, called the tangent problem: given v s , v f , m, such that the boundary and initial conditions of the nonlinear problem are satisfied, calculate δv s , δv f , δm, solution of the following system of linear equations, where D vs S etc. represent Fréchét derivatives at the point v s , v f , m, and the system must be solved using boundary and initial conditions of the same type of the nonlinear problem, but with homogeneous (null) data.
In [20] an approximate yet explicit expression of the tangent problem is provided. More precisely, the nonlinear problem is linearized around the configuration at rest, namely v s , v f , m = 0. As a result we have m = 0 and φ = φ 0 = 0. Concerning the fluid phase, Newtonian and incompressible behavior is assumed, which yields σ f (v f ) = 2µ f ε(v f ), being ε(v) = 1 2 (∇v + ∇v T ) the symmetric deformation gradient. As in [20], we denote by v s , v f , m the increments with respect to such state and use an additive decomposition of the free energy, with a Saint-Venant Kirchhoff component for the mechanics and a quadratic potential for the volumetric deformation of the solid phase J s , which reads where E = 1 2 F T + F + F T F denotes the Green-Lagrangian strain tensor, also µ, λ are the Lamé constants and κ s is the bulk modulus. Under small deformations we have that E ≈ ε(u s ) and J ≈ 1 + ∇ · u s , which give where C is a fourth order constant tensor (symmetric, positive definite), known as Hooke tensor.
In the linearized setting it is possible to reformulate the problem in terms of the (more commonly used) variable p instead of the added mass. As a result, the approximate tangent problem for the configuration at rest reads as follows: find u s , v f , p such that For simplicity, in what follows we assume θ = 0. The system (2.1) is closed with appropriate boundary conditions naturally following from the ones of the nonlinear problem. For the sake of clarity we report them here

Numerical approximation of the tangent problem
We start form the time discretization, based on a simple backward Euler approach. We will discuss later on how higher order time discretizations are also viable and the resulting discrete problem maintains its fundamental traits, such that the numerical solvers developed in what follows will still be applicable.
We consider a partition of the time interval of interest (0, T ), given by 0 = t 0 < t 1 < ... < t n < ... < t N = T with, for simplicity, constant time step size ∆t = t n − t n−1 . The temporal derivatives within the model (2.1) are approximated by finite differences We assume that besides the initial data the first time step has been already determined. From the second time step the fully dynamic linearized model can then be approximated by the Implicit Euler discretization using the above finite difference approximations: For n ≥ 2, given u n−1 Such problem must satisfy the same boundary conditions of (2.1) at each time t n , where f n , u D,n s etc. denote suitable approximations of the external problem data at time t n . In what follows we will apply the lifting technique to nonhomogeneous Dirichlet boundary data (in other words, a change of variable is introduced, by subtracting from the solution a function that is regular enough and equal to the prescribed datum on the boundary, such that the problem for the new variable is transformed into a standard homogeneous Dirichlet-type problem). In this way, all the forcing terms of the problem (volume forces and surface forces/data) will be implicitly represented in the volume term f n without significant loss of generality. Initial conditions are prescribed as in (2.2e)-(2.2h) by suitably approximating the initial data. Finally, we stress that the mass conservation equation has been divided by the constant fluid density ρ f in order to highlight an apparent symmetry between the equations. Remark 2.1 (Higher order time discretization). Applying alternative diagonally implicit Runge-Kutta schemes results in coupled systems of governing equation of similar type as (2.3). Material parameters possibly have to be scaled appropriately, and the right hand side source terms may then also include further previous data. However, we stress that the analysis of the splitting in this work does not depend on the choice of the time discretization similarly as in [38] and it could possibly be used also in the framework of space-time finite elements, used for example in [39] for the approximation of Biot poroelasticity system.
Let V , W , Q denote suitable function spaces for the solid displacement, fluid velocity, and fluid pressure, respectively, at discrete time t n , incorporating in particular homogeneous essential boundary conditions on the relevant boundaries, Remark 2.2. We note that in the weak formulation the constraints (1 − φ)u , φv ∈ H(div; Ω) are formally required for the corresponding terms to be well-defined. This is in fact a regularity condition on φ, where it is sufficient to consider that both φ and 1/φ belong to W s,r (Ω) with s > d/r and r ≥ 1. More details in [21].
Then the canonical weak formulation of (2.1) reads: The numerical discretization in space is based on the Galerkin projection of the solution Also, all the physical parameters of the tangent problem are assumed to be constant in time and uniformly bounded in space. Under these assumptions, the fully discrete version of the problem is formally equivalent to (2.4 and the test functions are taken in the same discrete spaces. Then, to avoid redundancy of notation, we will identify problem (2.4) with the fully discrete one and we will omit to specify the subindex h, unless strictly necessary. The finite element spaces will be kept generic throughout the derivation of the numerical solution algorithms, until the discussion of suitable numerical examples that will refer to precise choices of such spaces.

A two-way split inspired by alternating minimization
In the following, we introduce an iterative splitting for the semi-discrete approximation (2.4), decoupling the momentum equation for the solid phase and the remaining two equationsthe method will be directly applicable for the fully discrete approximation. The systematic construction (and later analysis) of the decoupling scheme is based on the general framework introduced in [37]. The central idea is to first equivalently rewrite the semi-discrete approximation as an auxiliary convex minimization problem, and second apply alternating minimization to derive a robust block-partitioned solver. Ultimately, reformulated in terms of (2.4), the final scheme is closely related to the undrained split for the quasi-static Biot equations [40], adding a div-div stabilization term to the momentum equation for the solid phase.
In what follows we require the following assumption, which has two modeling consequences: On one side, it rules out the possibility of considering the incompressible limit (κ s → ∞) with this approach, and on the other one it imposes that the domain cannot be composed only of fluid (φ = 1).

Problem formulation as convex minimization
We choose u n s and v n f as primary variables. Under Assumption 1, the mass conservation equation can be inverted with respect to the pressure, such that This allows to formally reduce (2.4) to a two-field formulation for the solid displacement and fluid velocity: where the momentum equation for the fluid has been scaled by ∆t, and g n s ∈ V and g n f ∈ W are defined by The symmetry and uniform coercivity of the governing equations (3.2) identify those as the optimality conditions of a block-separable convex minimization problem. Namely it holds with the energy given by The formulation (3.3)-(3.4) serves as basis for the succeeding construction of a robust split scheme for (2.4).

Robust splitting via alternating minimization
Following the approach of [37], we propose an iterative block-partitioned solver for the problem (2.4). In particular, the fundamental alternating minimization algorithm is applied to the equivalent variational formulation (3.3), cf. Alg. 1 for the definition of a single iteration with index k. By construction, the approximate solution consecutively minimizes the system energy J .
Algorithm 1: Iteration k ≥ 1 of the alternating minimization applied to (3.3) By introducing a pressure iterate, p n,k , analogously to (3.1) Alg. 1 can be equivalently reformulated in the context of the three-field formulation (2.4).
In particular, the k-th iteration of the iterative splitting scheme decouples in two steps. In the first step, a div-div stabilized momentum equation for the solid phase is solved: In the second step, the mass conservation and fluid momentum equations are solved: given In the remainder of this paper, we refer to the scheme (3.6)-(3.7) as the alternating minimization split.
In equation (3.6) the term naturally emerges to stabilize, namely control, the increment of u n s at each iteration, such that convergence is guaranteed. This will be proved in more detail in Section 5.1.

Diagonally L 2 -stabilized two-way split
Another prominent class of block-partitioned solvers for coupled problems with saddle-point structure are L 2 -stabilized splits, which have been successful especially in the context of coupled flow and mechanics. The so-called fixed-stress split for the quasi-static Biot equations [40,41], for instance, decouples solid and flow computations and employs simple L 2 -stabilization of the mass conservation equation, resulting in unconditional convergence [24,26]. It is worth mentioning that in practice the fixed-stress split often is superior to the undrained split [41], also motivating further investigation in the context of this work. Moreover, in the context of thermoporoelasticity diagonal L 2 -stabilization has been recently investigated for coupled systems consisting of more than two equations [42]. In particular, it has been observed that adding stabilization to multiple equations can be beneficial.
In the following, we present a diagonally L 2 -stabilized two-way split for (2.4). At first, we allow for stabilization of any of the three equations, introducing three stabilization parameters: β s (tensor-valued), β f (tensor-valued), β p (scalar-valued), potentially varying in space.
A single iteration of the splitting scheme is composed of two steps. Let k ≥ 1 denote the iteration index. Following the idea of the fixed-stress approach, the stabilized fluid flow problem is solved first; this is not required for convergence. The L 2 -stabilized fluid flow step reads: given The second (L 2 -stabilized solid mechanics) step reads: Finding physically motivated stabilization parameters β s , β f , β p , as for the original fixedstress split, mainly due to an increased complexity of the model (additional terms associated to dynamics) compared to the quasi-static Biot equations. For the same consideration, a simple discussion utilizing convex minimization and alternating minimization as in Section 3 is problematic, while being possible for the original fixed-stress split [37]. However, to get a better intuition, a closer look at (skew-)symmetries of the governing equations (2.4) and (partial) Schur complements is useful; in particular, the skew-symmetry of the u s − p coupling motivates positive stabilization of the mass conservation equation, and the symmetry of the u s − v f coupling suggests negative stabilization of the momentum equation of the solid phase. Since after all the u s − (v f , p) coupling is neither symmetric nor skew-symmetric, these observations merely lead to inaccurate insight. Instead, a succeeding convergence analysis in Section 5.2 is going to suggest practical, potentially vanishing values for the parameters, which eventually lead to unconditional stability.

A priori convergence analysis of the proposed two-way splits
In this section, we address the a priori convergence analysis of both the alternating minimization split (3.6)-(3.7) and the diagonally L 2 -stabilized two-way split (4.1)-(4.2), proposed in Section 3 and Section 4, respectively. The two primary goals are to (i) prove the linear convergence of the alternating minimization split, and (ii) determine ranges and specific practical values for the stabilization parameters employed within the diagonally L 2 -stabilized two-way split ensuring convergence. The two goals will be achieved using different techniques. For item (i) the interpretation of the alternating minimization split as alternating minimization applied to a strongly convex minimization problem is extensively exploited, allowing for the systematic application of sharp abstract convergence results from the literature; for item (ii) a slightly more technical approach is chosen due to the fact that the two-way split (4.1)-(4.2) does not fully conform with any (skew-)symmetry. In particular, we relax the classical (quotient) convergence criterion by means of the root convergence criterion, briefly called r-convergence, see for example [43]. More precisely, we formulate a general convergence criterion (based on relative stability) that turns out to be a sufficient condition for the r-convergence of the proposed iterative method.

Convergence analysis of the alternating minimization split for the tangent model
Guaranteed linear convergence of the alternating minimization split (3.6)-(3.7) is a direct consequence of its interpretation as alternating minimization applied to a (strongly) convex optimization problem, cf., Section 3.2 and e.g., [44]. Furthermore, using simple yet largely sharp abstract convergence results for alternating minimization in a Banach space setting, cf. [45], an upper bound of the rate of convergence can be provided. In the aforementioned work, it is showed that in each of the two steps of the alternating minimization, the energy values of the iterates are sequentially decreased with the decrease merely governed by convexity and continuity properties of the restricted minimization problems. Since the energy J is quadratic, energy differences relative to the optimum will directly translate to distances to the solution, measured in the problem-specific norm induced by the Hessian of the energy (at an arbitrary point). We In order to estimate the rate of convergence, we introduce a technical, a priori material constant γ ≥ 0, given by with κ m > 0 denoting the smallest eigenvalue of the permeability tensor k f , K dr,φ 0 ,min > 0 being a porosity dependent bulk modulus type constant given by and C Korn,1 , C Korn,2 > 0 taking on the role of generalized Korn/Poincaré constants, defined as the minimum positive numbers such that It is fair to assume that C Korn,1 and C Korn,2 are closely related to the inverse of the drained bulk modulus K dr := K dr,0,min . Finally, focusing only on the fully transient model, the linear convergence result for the alternating minimization split scheme reads as follows.
Theorem 5.1 (Linear convergence of the alternating minimization split ). Let (u n s , v n f ) ∈ V ×W , n ≥ 2, denote the solution to (3.3), and let (u n,k s , v n,k f ) ∈ V ×W , k ≥ 1, denote the corresponding approximation defined by Alg. 1. Let γ ≥ 0 be the material constant defined as in (5.1). Then, for all k ≥ 1, it holds that The convergence result is similar as for the undrained split for the quasi-static Biot equations, cf. [24]. In particular, the theoretical result suggests degenerating convergence for nearly incompressible and impermeable media. In contrast to the quasi-static Biot equations, porosity heterogeneities may also affect the performance of the splitting scheme, as the material constants γ 1 and γ 2 depend on the spatial gradients of φ 0 . However, a numerical test in Section 6.3 does only show a weak influence.
The proof of Theorem 5.1 is a plain application of the following abstract convergence result for the alternating minimization, here specifically formulated in terms of Alg. 1.
Lemma 5.2 (Convergence of the alternating minimization [45]). Let | · |, | · | s , and | · | f denote semi-norms on V × W , V , and W , respectively, such that: Let J : V × W → R be Frechét differentiable with DJ denoting its derivative such that: The energy J is strongly convex wrt. | · | with modulus σ > 0, i.e., for all u s ,ū s ∈ V and v f ,v f ∈ W it holds that The partial functional derivatives D us J and D v f J are uniformly Lipschitz continuous wrt. | · | s and | · | f with Lipschitz constants L s and L f , respectively, i.e., for all , and let (u n,k s , v n,k f ) ∈ V × W denote the corresponding approximation defined by Alg. 1. Then, for all k ≥ 1 it follows that With this, we are able to prove Thm. 5.1.
Proof of Thm. 5.1. In order to apply Lemma 5.2, we verify the conditions (A 1 )-(A 3 ). First of all, we note that the energy J is quadratic. Since | · | is induced by the Hessian of J , i.e., (for arbitrary (u s , v f ) ∈ V × W ), the convexity property (A 2 ) is satisfied with σ = 1. Similarly, by defining | · | s and | · | f on V and W , respectively, as partial Hessians of J |u | 2 It remains to examine (A 1 ). In the following, we show that one can choose β s = β f = (1+γ) −1 , i.e., it holds For both estimates, the following inequality will be of help Indeed, for T 1 , using the product rule, the Cauchy-Schwarz inequality and Young's inequality, we obtain for all ζ > 0 Further, employing the definitions of K dr,φ 0 ,min and C Korn,1 , see (5.3) and (5.4), it follows that Similarly, employing the definitions of κ m , the smallest eigenvalue of k f , and C Korn,2 , see (5.5), for T 2 it holds and thereby (5.8) follows. Finally, we show (5.7). By definition of | · | s it holds that Hence, (5.7a) follows from (5.8). By (i) definition of | · | f , (ii) suitable addition with zero, and application of the Cauchy-Schwarz and Young's inequalities, and (iii) (5.8), it holds Hence, we obtain (5.7b), and thereby (A 1 ). Ultimately, the assumptions of Lemma 5.2 are satisfied, and it follows for all k ≥ 1 that Moreover, since J is quadratic, (u n s , v n f ) is a local minimum of J , and | · | is induced by the functional Hessian of J via (5.6), we have that

Convergence analysis of the diagonally L 2 -stabilized split for the tangent model
The essence of the diagonally L 2 -stabilized split (4.1)-(4.2) is the decoupling of the mechanical displacement from the remaining variables (fluid pressure and velocity). Such a split does neither fully conform with a symmetry nor a saddle point structure of the governing equations. In view of a convergence analysis aiming at employing some contraction argument or similar, it therefore cannot be expected that all coupling terms can be simultaneously canceled by suitable testing as often done, cf., e.g., [26]. To mitigate this complication, the concept of relative stability will be exploited instead, allowing for a simpler discussion of the coupling terms. In the following, the analysis is presented in two steps: (i) a central abstract convergence result for positive real-valued sequences satisfying a relative stability property is introduced; (ii) the result is applied to the diagonally L 2 -stabilized split (4.1)-(4.2) to show a priori convergence.

Abstract convergence criterion based on relative stability
Consider a real-valued (positive) sequence {x k } k ⊂ R + satisfying the stability property: without any additional requirement for the stability constant c. We call this property the relative stability criterion for the sequence {x k } k ⊂ R + . This criterion ensures r-linear convergence for subsequences (still wrt. the original sequence), a weaker form of standard r-linear convergence, covering both contractive and certain non-contractive sequences. For the proof of Lemma 5.3, we state the following auxiliary result.
Proof. Let k ∈ N and ε > 0 be arbitrary but fixed. Assume without loss of generality that x k > 0. Then the assertion follows by contradiction: Assume it holds x k+i > εx k , for all i = 1, ..., 1 cε ; we conclude that it holds which contradicts (5.12).
Proof of Lemma 5.3. The idea of the proof is to employ Lemma 5.4 and construct a subsequence of {x k } k , which is linearly (first order) quotient converging, and then conclude r-linear convergence wrt. the original sequence. Assume without loss of generality that x 0 > 0. Let m ∈ N such that cm > 1, and let ε := 1 cm < 1, such that 1 cε = m. By Lemma 5.4 there exists some n 1 ∈ {1, ..., m} such that it holds x n 1 ≤ εx 0 . Analogously, for any i = 2, ..., there exists some n i ∈ {1, ..., m}, satisfying x i−1 j=1 n j +n i ≤ εx i−1 j=1 n j .
Next, we define {k l } l ⊂ N by setting k l := l j=1 n j for all l ∈ N. Since ε < 1 and n j ≤ m for all j, it holds that For {x k l } l , we conclude for, so far, arbitrary m ∈ N. Minimizing the right hand side wrt. m ultimately yields the assertion.

Convergence analysis of the diagonally L 2 -stabilized split based on the concept of relative stability
In the following, we establish linear convergence of the diagonally L 2 -stabilized two-way split (4.1)-(4.2). The primary aim of the analysis is to determine ranges for the stabilization parameters β s , β f and β p , which a priori guarantee convergence; in addition, we are going to suggest a practical (for simplicity of the presentation not necessarily optimally tuned) set of values. The reader interested in the analysis of optimal convergence rate is referred to analogous studies of the fixed-stress split for the quasi-static Biot equations [27].
For the convergence analysis, the concept of relative stability and r-linear convergence for subsequences introduced in the previous section is applied. Ultimately, the final result states that it is sufficient to stabilize the mass conservation equation along the lines of the fixedstress split for the quasi-static Biot equations [26,27,40,41], in order to guarantee convergence. Additional destabilization, i.e., negative stabilization, of the momentum equation for the solid phase theoretically improves the convergence speed. Fluid (de-)stabilization does not further improve the convergence rate.
To ease the presentation of the analysis, we introduce two notations: (N1) Weighted squares · 2 A , defined by ω 2 A := Aω, ω , where ω can be a tensor-, vectoror scalar-valued function on Ω, and the weight A is a (potentially non-positive definite) function on Ω with adequate dimensionality such that the above definition is well-defined.
Theorem 5.5 (Relative stability and convergence of the diagonally L 2 -stabilized two-way split).
, and d k p := p n,k − p n,k−1 denote increments for k ≥ 1. Furthermore, let K dr,φ 0 ,min and C Korn,1 as defined in (5.3) and (5.4), resp.; let δ 1 > 0 and δ 2 ∈ (0, 2) be tuning parameters; and let the stabilization parameters satisfy where A B for tensor-valued maps A and B on Ω iff. A − B is uniformly positive definite. Then the scheme (4.1)-(4.2) satisfies a relative stability criterion of type (5.12), namely where β s and β p denote augmented stabilization parameters (introduced for simpler presentation) is uniformly positive, subsequences of d k u , d k v , d k p r-linearly converge to zero, in the sense of Lemma 5.3.
Proof. The proof is organized in five steps, starting with governing equations for increments.
Increment equations By subtracting (4.1)-(4.2) at iteration k and k − 1, k ≥ 2, we obtain Testing with current increments Testing and summing (5.15) with u = d k u , v = ∆t d k v , and p = ∆t d k p , and finally summing over indices k = m + 1, ..., M , for arbitrary m < M , yields (1−φ 0 ) 2 κs + T 1 = T 2 + T 3 (5.16) (employing notation (N2) and) with We discuss the terms T 1 , T 2 and T 3 separately. For the stabilization term T 1 , we apply binomial identities of type (a − b)a = 1 2 a 2 − b 2 + (a − b) 2 and telescope sums, resulting in (5.17) (employing notation (N1)). For the coupling term T 2 we apply summation by parts, leading to For the coupling term T 3 , simple expansion and reformulation, aiming at constructing quadratic terms present on the left hand side of (5.16), and gathering those, respectively, results in We discuss the coupling terms T 4 , T 5a and T 5b separately in the two following steps. Finally, after choosing β s , β f and β p satisfying (5.13) (in particular translating to β s 0), and dropping several positive terms in (5.24), we obtain the stability result After all, relative stability in the sense of (5.12) can be deduced for any choice for δ 1 > 0 and δ 2 ∈ (0, 2), since m and M have been chosen arbitrary. By this the assertion follows.
Remark 5.6 (Incompressible media). We note that in contrast to the alternating minimization split (3.6)-(3.7), the diagonally L 2 -stabilized two-way split (4.1)-(4.2) remains well defined in the extreme case of (quasi-)incompressible solid material, i.e., (1−φ 0 ) 2 κs = 0. According to the above theory, convergence is not guaranteed anymore, yet still may be possible in practice, see also examples in Section 6.
We close this section with suggesting a practical set of stabilization parameters guided by the previous convergence analysis. We emphasize that one could optimize the effective stability constant in (5.14) wrt. δ 1 , δ 2 , β s , β f , β p ; however, theoretical optimality does not necessarily result in practical optimality, cf. [27] for an applicable discussion.
Remark 5.7 (A practical set of stabilization parameters). We assume heterogeneities of the porosity are not crucial and pretend the porosity is constant. Then it is C Korn,1 = 0 and K −1 dr,φ 0 ,min = (1−φ 0 ) 2 K dr , where K dr = K dr,0,min denotes the standard drained bulk modulus. Moreover, we choose the values δ 1 = δ 2 = 1 in order to balance similar terms on both sides of (5.14) and follow the suggestion of the stability property to choose the stabilization parameters as "small" as possible. This results in the set which leads to destabilization of the momentum equation of the solid. However, we also highlight that merely utilizing pressure stabilization and setting β s = β f = 0 does also result in guaranteed convergence, in the style of the fixed-stress split for the quasi-static Biot equations.

Numerical tests for the convergence of the proposed splitting schemes
The aim of this section is to assess the performance of the proposed splitting schemes, the alternating minimization split, cf. Section 3, and the diagonally L 2 -stabilized two-way split, cf. Section 4, and to compare it with the theoretical convergence results in Section 5. In particular, we consider three test cases and perform an extensive parametric study for various choices of model parameters and stabilization values based on the above analyses, in addition to similar ad-hoc choices motivated by the analyses or experience of the closely related splitting schemes for the Biot equations.
As test problems, we use two classic problems, the swelling [19,20,46] and footing [27,36,47] problems. In addition, we consider a perfusion-like problem as a reference for biomedical applications. We note that each problem is loaded on a different equation: the swelling on the fluid, the footing on the solid and the perfusion on the mass balance.
We first present a sensitivity study with respect to the physical parameters for both alternating minimization and L 2 -stabilized splits independently based on the swelling test. Then, we provide a detailed comparison between both methods in all the described test problems in combination with Anderson acceleration. At the end of this section, we also compare the performance of the split scheme that results most effective, with a monolithic solution approach for the linearized problem, which may be considered to be the gold standard solution strategy. This final test sheds light on the competitiveness of the proposed schemes when used for realistic scenarios.
All numerical examples have been performed using the FEniCS project [48,49], and convergence is measured in terms of the relative residual (for larger certainty absolute residuals are not considered).

Definition of the test cases
The swelling test This test consists of a 2D slab Ω = (0, L) 2 , L = 10 −2 , in absence of volume forces and simulated in the time interval (0, 1), with time step ∆t = 0.1. It is subject to an inflow φ 0 (2µ f ε(v f ) − pI) n = −p ext n, p ext (t) = 10 3 (1 − exp(4t 2 )) on the left and null stress on the right, whereas above and below it uses a no-slip boundary condition v f = 0. The boundary conditions for the solid are sliding on the bottom and left sides, whereas the rest of the boundary is of null traction type (see Figure 1 (a)). We note that these conditions are not physical because the fluid boundary pressure should act as force on the solid as well, but we keep the proposed scenario to have this test being loaded only on the fluid equation. We have indeed tested this and observed that it presents no impact on the following study. The following default parameters (from [20]) are used (unless otherwise specified): ρ f = ρ s = 1000, µ f = 0.035, λ s = 711, µ s = 4066, κ s = 10 3 , k f = 10 −7 I, φ 0 = 0.1, all in SI units; in addition Ω is discretized using 10 elements per side. Let us denote by X k h (Ω) the Lagrangian k-th order finite element space defined on a quasi-uniform mesh of Ω of characteristic size h. Unless otherwise specified, the default finite element spaces used are: first order Lagrangian elements for the solid, V h = X 1 h (Ω) and Taylor-Hood elements for the fluid-pressure system, , which satisfy the weighted inf-sup condition partially but are more useful in practice [21]. We name this choice of elements with the shorthand notation P1/P2/P1. Finally, a relative tolerance of 10 −8 was used with respect to the ∞ norm of the residual, where all sub-problems are solved using GMRES with a relative tolerance of 10 −8 as well, preconditioned with with an incomplete LU factorization with 3 levels of depth (ILU(3) [50]).

Anderson acceleration
One key aspect of both proposed schemes is that they can be interpreted as fixed point iterations.
Although they feature in general lower convergence rates than Newton methods, they have acquired higher interest recently, also due to the development of acceleration schemes. In particular, we focus on the Anderson acceleration, which can be interpreted as a multisecant scheme, or as a preconditioned GMRES iterative method [30]. As shown later on in Tables 8, 9, 10, acceleration techniques greatly improve the performance of the proposed split schemes, by increasing their robustness with respect to varying loading conditions and significantly reducing the iteration count. In practice, using Anderson acceleration is a necessary choice to effectively use the described split schemes in demanding scenarios. In general, consider a vector-valued function g : R N → R N and the sequence By defining f k = g(x k ) − x k , Anderson acceleration of order m, abbreviated by AA(m), is given as follows: For iteration k, set m k = min{m, k} and F k = (f k−m k , ..., f k ). Compute α k = (α k 0 , ..., α k m k ) that minimizes min α=(α 0 ,...,αm k ) and then compute the next element as The order m of the scheme is usually referred to as depth, due to the use of m previous iterations. We implement this method by recasting (6.1) as an unconstrained least-squares problem, and then invert its optimality conditions using the QR factorization to avoid possible ill-conditioning of the normal equations [50].

Numerical tests for the alternating minimization split
In this section we present three numerical tests on the alternating minimization split (named Alt-min in the tables), with the aim of verifying the robustness of the scheme with respect to the parameters N = κs (1−φ 0 ) 2 , k f and highly oscillatory porosities φ 0 . As test case we adopt the swelling test described above. We consider three varying parameters where L = 10 −2 is the side length, and the permeability is treated as a scalar for simplicity, namely k f = κ f I. For each parameter aside of default parameters otherwise, we present the average number of splitting iterations throughout the simulation required for convergence in Table 1. We observe that the performance of the alternating minimization split is particularly sensitive to the bulk modulus κ s , and small permeabilities make the problem much more difficult to solve. Instead, the dependence on oscillating porosity is moderate. The results are in accordance to

Numerical tests for the diagonally L 2 -stabilized split
In this section, we study the sensitivity of the performance of the diagonally L 2 -stabilized split (named L 2 S in the tables) with respect to different combinations of physical parameters. Precisely, we use the swelling test with default coefficients, and we vary the following ones Additional tests address the influence of the ratio between the elasticity and the permeability. For this, we fix the permeability and increase the drained bulk modulus K dr = λ + 2µ d , d = 2 by scaling both Lamé parameters by the same factor.
The analysis in Section 5.2 yields the interesting fact that the solid momentum equation can be destabilized. Therefore, we compare different stabilization parameters, also ones excluded by the theory in order to investigate the theoretically suggested parameter ranges. In particular, we apply the L 2 stabilized two-way split (4.1)-(4.2) using stabilization parameters of type with different scaling factorsβ s ,β f ,β p , from now on denoted by L 2 Sβ s,βf ,βp . Considered scaling factors are listed in Table 2. Splitting iterations are terminated via the tolerance tol res = 10 −8 .
As in the previous test, performance is measured in terms of the average number of splitting iterations throughout the entire simulation, with non-convergence established whenever a solver requires more than 200 iterations.
IDβ sβfβp Description Covered by Thm. 5.5 Table 2: Considered stabilization settings in the context of the diagonally L 2 -stabilized split.
Although the analysis, cf. Thm. 5.5, does not reveal any dependence on the particular discretization, it is developed under the underlying assumption that the discrete problems are uniquely solvable. To investigate potential effects of stability of the function spaces onto the stability of the splitting, we consider progressively unstable approximation spaces, namely P1/P2/P1 and P1/P1/P1 elements for displacement, velocity and pressure, respectively.

Dependence on solid bulk
In Table 3, the iteration counts for varying κ s are displayed. We observe that for P1/P1/P1 elements no set of stabilization parameters enables convergence for larger κ s ; we note that for increasing κ s , the uniform stability of the fluid-pressure system is lost. In contrast, the use of P1/P2/P1 elements adds uniform stability to the discretization and finally also uniform robustness to any of the stablized splittings. For a non-dominating u s − v f coupling, destabilization of the solid momentum equation does not make a big difference.   Non-convergence denoted by -.

Dependence on permeability
In Table 4, the iteration counts for varying k f = κ f I are displayed. Here, a maximal count of 500 splitting iterations is used for better understanding the dependence on the permeability. Lower permeability makes the problem more difficult to solve. The reasons for this are: (i) decreasing the permeability leads to ill-conditioning of the u s -v f block; and (ii) for lower permeabilities the ellipticity of the u s -v f block loses its dominance, and instead the L 2 -type contribution has a much bigger influence. Destabilization of u s seems to effectively address the first issue. In fact, it results in significantly improving the performance, compared to mere p-stabilization, which alone fails to lead to unconditional robustness. This, on the one hand, nicely verifies the theory in Thm. 5.5. On the other hand, it indicates that suitable destabilization successfully imitates approximating the Schur complement of the L 2 -type contribution of the u s -v f block; the comparison of conservative and aggressive destabilization illustrates the potential gain but also sensitivity of destabilization. Since L 2 -stabilization of the mass conservation equation does not address the L 2 -type contribution of the u s -v f block at all, unconditional robustness cannot be expected without an additional differently scaled stabilization approach, ultimately mitigating the second issue.
Comparing the results for the P1/P1/P1 and P1/P2/P1 discretizations, we note that inf-sup stability in the fluid allows for a significant improvement on the performance. Also, in contrast to the unstable case, destabilizing the u s equations greatly improves performance.  Table 4: Average iteration count of the L 2 -stabilized solvers for a varying k f in the swelling test. Non-convergence denoted by -(more than 500 iterations in this case).
We note that this method is very sensitive to low permeabilities. In previous studies, e.g. [29], Anderson acceleration has been shown to successfully increase robustness of stabilized iterative solvers. So we present the iteration counts for the same test but using Anderson acceleration with a depth of 5 in Table 5. We note that not only there is a significant decrease in the number of splitting iterations required (up to ca. 80% for very low permeabilities), but it also enables the convergence of configurations which have previously not converge, again verifying previous observations. As long as the permeability is not too low, again aggressive u s stabilization leads to the best performance.  Table 5: Average iteration count of the L 2 -stabilized solvers for a varying k f in the swelling test using Anderson acceleration with depth 5. Non-convergence denoted by -(more than 500 iterations in this case).

Dependence on densities
In Table 6, the iteration counts for varying ρ s = ρ f are displayed. We observe that for very large densities the problem starts to become more difficult to solve. To explain, increasing densities (merely) raise the second issue mentioned in Section 6.4.2; in particular, as expected, destabilizing the solid equation does not yield any improvement, in contrast to the previous test. Iteration counts are identical for P1/P1/P1 and P1/P2/P1 elements. Thus, only the former is presented.  Table 6: Average iteration count of the L 2 -stabilized solvers for a varying ρ s = ρ f in the swelling test. Non-convergence denoted by -.

Dependence on drained bulk modulus
In Table 7, the iteration counts for varying K dr (with same Poisson ratio) are displayed. We observe that lower drained bulk modulus is associated to higher iteration counts. This can be explained along the lines of the discussion of the dependence on the permeability, cf. Section 6.4.2, since a lower drained bulk modulus leads to dominance of the L 2 -type contribution of the u s -v f block. Therefore, as expected, (aggressive) destabilization is beneficial. Additionally, a lower drained bulk modulus leads to a stronger coupling strength, and in accordance to Theorem 5.5, to a deteriorating convergence rate. Again, inf-sup stability of the discretization of the fluid-pressure coupling enables slightly improved results, especially for low bulk modulus.   Table 7: Average iteration count of the L 2 -stabilized solvers for a varying K dr in the swelling test. Non-convergence denoted by -.

Comparison of the alternating minimization and L 2 -stabilized splits
The previous two sections allow for a first comparison of the two proposed schemes. In particular, two conclusions on the respective limitations can be made: (i) for increasing solid bulk modulus, the alternating minimization split quickly deteriorates, whereas the L 2 -stabilized split remains robust; and (ii) for lower permeabilities, the performance of both schemes deteriorates, but the alternating minimization split in fact better handles the limit of very low permeabilities.
In this section, we continue the comparison of the two proposed schemes, now based on all the three suggested test cases with the parameters given in their description, enjoying different problem characteristics. The focus of the following study will also be to assess the impact of actual inf-sup stability, given for a Taylor-Hood like P2/P2/P1 discretization, opposed to the previously considered P1/P2/P1 discretization. Moreover, having observed the improving effect of Anderson acceleration in Section 6.4.2, we follow this lead and also investigate the performance of the accelerated splits, this time also for the alternating-minimization. We also consider only the L 2 S −0.5,0,1 as it is the one suggested by the analysis and it exhibits an overall more robust performance.
For the swelling test, we additionally consider two bulk moduli, κ s ∈ {10 4 , 10 8 }. Results are presented in Table 8. We observe that the inf-sup stability of the displacement plays no role, and the diagonally L 2 -stabilized split proves very robust in all the tested scenarios, performing significantly better than the alternating minimization split. For the first, Anderson acceleration barely leads to improvement due to already low iteration counts; for the latter convergence can be significantly accelerated for the lower bulk modulus. For high bulk modulus, not even Anderson acceleration enables convergence.  We present the results of the footing test in Table 9. We note that in this test the alternating minimization scheme exhibits lower iteration counts. Its success can be explained by the lower bulk modulus used, and instead the initial failure of the L 2 -stabilized scheme is due to the permeability, which is very low. This case presents localized displacements at Γ foot , which are more affected by numerical locking, which justifies the increased iteration count in the case of the P2/P2/P1 discretization.  Table 9: Average iteration count for all tested scenarios in the footing test. None stands for the plain splits; AA(m) stands for additional application of Anderson acceleration with depth m.
The results of the perfusion test are presented in Table 10. The behavior of this test is similar to the swelling one, with the L 2 -stabilized split exhibiting a robust performance, which is further improved by the use of acceleration. The alternating minimization split instead presents difficulties in attaining convergence without acceleration, which can be explained by the use of a large bulk modulus. Similarly to the swelling test, the inf-sup stability of the displacement effectively plays no role.

Comparison of splitting versus monolithic approaches
In this section we present a comparison, in terms of computational time, between the proposed splitting schemes and a monolithic approach. We consider the swelling test and we choose the L 2 S −0.5,0,1 (labelled L 2 S) as it yields the best performance for this problem. The default stopping criterion for GMRES iterations is adopted for the monolithic scheme, with a relative  tolerance equal to 10 −8 . For the splitting scheme, the convergence tests for the linear system solved at each iteration is slightly relaxed, up to 10 −6 , but the (relative) tolerance of the stopping criterion for the iterative splitting scheme is also set to 10 −8 , on the ∞ norm of the residual. We compare the computational cost, measured by the average wall time per time step, calculated on a sequence of five consecutive time steps. Both formulations are solved using P1/P2/P1 finite elements, and the number of degrees of freedom is controlled by the number of nodes on each side of the domain. The results of the comparison are reported in Table 11. The iterative schemes exhibit a better scaling with respect to the number of degrees of freedom. In particular, for problems with over 10 5 degrees of freedom (given by using 100 or more elements per side on the square domain) the wall time of the split scheme is consistently lower than of the monolithic approach. Also, the ratio between both solution times decreases monotonically with respect to the degrees of freedom as shown in the last column of the table, meaning that in this test case the superiority of iterative splitting schemes increases with the discrete problem size, which makes them a competitive solution strategy for addressing realistic scenarios, especially when considering tailored, possibly scalable preconditioners for the single subproblems.

Discussion and Conclusions
In this work we have developed splitting schemes for the linearized poromechanics problem studied in [20,21], namely the alternating minimization split and the diagonally L 2 -stabilized split. As the choice of a splitting scheme strongly depends on the application of interest, due to the strong dependence that the performance of each scheme has on the parameters, we tested the proposed methods on several benchmark problems. The conclusions of this work arise from both theoretical analysis and numerical experiments. From the standpoint of theoretical convergence properties, we observe that the effectiveness of a splitting scheme hinges on the assumptions used for the convergence analysis and the corresponding stabilization, if necessary. For instance, the alternating minimization scheme requires the algebraic inversion of the pressure, so it can be expected for it to deteriorate whenever this operation is not admissible ((1 − φ)/κ s → ∞). The diagonally L 2 -stabilized split can be interpreted as an approximate Schur complement method, where the L 2 -type contributions are not considered. This implies that it can be expected for such L 2 -stabilized schemes to present difficulties converging whenever the L 2 -type contributions are dominant, meaning small permeability or large densities. The analysis also provides the interesting possibility of destabilizing the solid momentum equation in the diagonally L 2 -stabilized scheme.
Such trends are confirmed by numerical experiments. The alternating minimization scheme performs very well in compressible scenarios but its convergence rate quickly deteriorates as the bulk modulus increases. The diagonally L 2 -stabilized split instead is robust with respect to the bulk modulus, so it should be preferred in (quasi-)incompressible regimes. The numerical experiments also confirm that the destabilization of the solid momentum equation yields good improvements of the convergence rate. Neither of the schemes is capable of handling large densities or small permeabilities -enhanced splitting schemes which successfully incorporate the L 2 -type contribution in the displacement-fluid velocity block are a topic of future research. Still, an improvement for the low permeability scenario can be seen by using inf-sup stable elements for the fluid-pressure block. This is an interesting property to be investigated, as it does not emerge in the analysis.
We have strengthened our splitting schemes with Anderson acceleration, which is a general method to improve the convergence of fixed-point iterations. It does not only improve the convergence of all methods tested, but it also enables convergence in scenarios in which it previously would not converge. Another feature of Anderson acceleration, particularly relevant in this framework, is that it reduces the influence of the stabilization parameters. This is indeed a fundamental aspect, as the user-defined choice and tuning of parameters represent a drawback of the presented methods.
Finally, we have compared the diagonally L 2 -stabilized split with a monolithic approach applied to the linearized problem. This study shows that for a sufficiently large size of the discrete problem, the iterative splitting approach is a competitive choice. Such methods may then be rightfully considered as effective options for solving realistic poromechanics problems applied to soft materials. Further investigations considering practical biomedical applications will be performed in the future.