On the optimization of the fixed‐stress splitting for Biot's equations

In this work, we are interested in efficiently solving the quasi‐static, linear Biot model for poroelasticity. We consider the fixed‐stress splitting scheme, which is a popular method for iteratively solving Biot's equations. It is well known that the convergence properties of the method strongly depend on the applied stabilization/tuning parameter. We show theoretically that, in addition to depending on the mechanical properties of the porous medium and the coupling coefficient, they also depend on the fluid flow and spatial discretization properties. The type of analysis presented in this paper is not restricted to a particular spatial discretization, although it is required to be inf‐sup stable with respect to the displacement‐pressure formulation. Furthermore, we propose a way to optimize this parameter that relies on the mesh independence of the scheme's optimal stabilization parameter. Illustrative numerical examples show that using the optimized stabilization parameter can significantly reduce the number of iterations.


INTRODUCTION
There is currently a strong interest in the numerical simulation of poroelasticity, ie, fully coupled porous media flow and mechanics. This is due to its high number of societal relevant applications, such as geothermal energy extraction, life sciences, or CO 2 storage, to name a few. The most commonly used mathematical model for poroelasticity is the quasi-static, linear Biot model. It is the coupled problem arising when considering the balance of linear momentum for the porous medium allowing for only small deformations (1) and mass conservation and Darcy's law for the fluid flow (2) (see, eg, the work of Coussy 1 ): find (u, p) such that −∇ · (2 (u) + ∇ · uI) + ∇p = f , where u is the displacement; (u) ∶= 1 2 ( ∇u + ∇u ⊤ ) is the (linear) strain tensor; and are the Lamé parameters; is the Biot-Willis constant; p and are the fluid's pressure and density, respectively; 1∕M is the compressibility constant; g is the gravitational vector; and is the permeability. The source terms f and S f represent the density of applied body forces and a forced fluid extraction or injection process, respectively.
A lot of work has been done concerning the discretization of Biot's equations (1) and (2). Various spatial discretizations, combined with the backward Euler method as temporal discretization, have been proposed and analyzed. We mention cell-centered finite volumes, 2 continuous Galerkin for the mechanics and mixed finite elements for the flow, 3-6 mixed finite elements for flow and mechanics, 4,7 nonconforming finite elements, 8 the MINI element, 9 continuous or discontinuous Galerkin, [10][11][12] or multiscale methods. [13][14][15] Continuous and discontinuous higher-order Galerkin space-time finite elements were proposed in the work of Bause et al. 16 Adaptive computations were considered, for example, in the work of Ern and Meunier. 17 A Monte Carlo approach was proposed in the work of Rahrah and Vermolen. 18 For a discussion on the stability of different spatial discretizations, we refer to the recent papers. 19,20 Independently of the chosen discretization, there are two popular alternatives for solving Biot's equations: monolithically or by using an iterative splitting algorithm. The former has the advantage of being unconditionally stable, whereas a splitting method is much easier to implement, typically building on already available, tailored, separate numerical codes for porous media flow and for mechanics. However, a naive splitting of Biot's equations will lead to an unstable scheme. 21 To overcome this, one adds a stabilization term in either the mechanics equation (the so-called undrained splitting scheme 22 ) or the flow equation (the fixed-stress splitting scheme). 23 The splitting methods have very good convergence properties, making them a valuable alternative to monolithic solvers for simulation of the linear Biot model (see, eg, the works of Both et al, 5 Kim et al, 21 Settari and Mourits, 23 and Mikelić and Wheeler 24 ). In the present work, we will discuss the fixed-stress splitting scheme. For other splitting schemes, see, for example, the works of Turska and Schrefler 25 and Turska et al. 26 After applying the backward Euler method in time to (1) and (2) and discretizing in space (using finite elements or finite volumes), one has to solve a fully coupled, discrete system at each time step. The fixed-stress splitting scheme is an iterative splitting scheme to solve this system. Let i denote the iteration index, and look for a pair (u i , p i ) to converge to the solution (u, p), when i → +∞. Algorithmically, one first solves the flow equation (2) using the displacement from the previous iteration, and then, one solves the mechanics equation (1) with the updated pressure and iterates until convergence is achieved. To ensure convergence, 5,21,24 one needs to add a stabilizing term L(p i − p i−1 ) to the flow equation (2). The free-to-be-chosen parameter L ≥ 0 is called the stabilization or tuning parameter. Choosing the value of this parameter is of major importance to the performance of the algorithm, because the number of iterations strongly depends on its value (see the works of Both et al, 5 Bause et al, 16 Both and Köcher, 27 Mikelić et al, 28 and Dana et al 29 ). Moreover, a too small or too big L will lead to slow or no convergence.
The initial derivation of the fixed-stress splitting scheme had a physical motivation 21,23 : one "fixes the (volumetric) stress," ie, imposes K dr ∇ · u i − p i = K dr ∇ · u i−1 − p i−1 and uses this to replace ∇ · u i in the flow equation. Here, K dr is the physical drained bulk modulus. The resulting stabilization parameter L, called from now on the physical stabilization parameter, is L phys = 2 K dr (depending on the mechanical properties and the Biot coefficient). In 2013, a rigorous mathematical analysis of the fixed-stress splitting scheme was performed for the first time in the work of Mikelić and Wheeler. 24 The authors show that the scheme is a contraction for any stabilization parameter L ≥ L phys 2 . This analysis was confirmed in the work of Both et al 5 for heterogeneous media using a simpler technique, and the same result was obtained for both continuous and discontinuous Galerkin higher-order space-time finite elements in the works of Bause et al 16 and Bause,30 implying that the value of the stabilization parameter does not depend on the order of the spatial discretization. The question of which stabilization parameter is the optimal one (in the sense that it requires the least number of iterations to converge) arises, and the aim of this paper is to answer this open question.
In a recent study, 27 the authors studied the convergence of the fixed-stress splitting scheme for different test cases with varying material parameters. They determined numerically the optimal stabilization parameter for each considered case. This study, together with the previous results presented in the works of Mikelić et al 28 and Both et al,5 suggests that the optimal parameter actually is a value in the interval [ L phys 2 , L phys ], depending on the data. In particular, the optimal parameter depends on the problem's boundary conditions and flow parameters, and not only on its mechanical properties and coupling coefficient. Nevertheless, to the best of our knowledge, there exists no theoretical evidence for this in the literature so far.
In this paper, we propose for the first time that the optimal stabilization parameter for the fixed-stress splitting scheme lies in the interval [ 2 , L phys ) and depends also on the fluid flow properties and stability properties of the spatial discretization. This is achieved through refining the proof techniques in the work of Both et al 5 to obtain an improved linear rate of convergence; minimizing this rate with respect to the stabilization parameter gives the "theoretical" optimal choice. Although the trends for the practical and the proposed theoretically optimal stabilization parameter are sound for varying material parameters, the theoretically calculated one does not show great practical promise in terms of being optimal (see the work of Storvik et al 31 for a supplementary numerical study). This is due to harsh bounds that have been used in the proof. Therefore, we propose a brute-force approach for optimizing the stabilization parameter, utilizing the newly found interval [ In contrast to previous works, the spatial discretization is required to be inf-sup stable, which essentially allows for the control of errors in the pressure by those in the stress. A novel consequence of our theoretical result is that under the use of an inf-sup-stable discretization, the fixed-stress splitting scheme also converges robustly in the limit case of incompressible fluids and impermeable porous media.
In Section 4, numerical experiments are performed, which show the soundness and efficiency of the proposed optimization technique. In particular, we show that the optimized stabilization parameter can be far superior to a naive choice among the classical stabilization parameters, L phys or L phys 2 .
To summarize, the main contributions of this work are as follows: • an improved, theoretical convergence result for the fixed-stress splitting scheme under the assumption of an inf-sup-stable discretization; • the derivation of an explicit interval for the optimal stabilization parameter, depending solely on the material parameters; • a brute-force approach for optimizing the stabilization parameter, relying on a nearly mesh-independent performance of the fixed-stress splitting.
We mention that the fixed-stress splitting scheme also can be applied to more involved extensions of Biot's equations, for example, including nonlinear water compressibility, 32 unsaturated poroelasticity, 33,34 the multiple-network poroelasticity theory, 35,36 finite-strain poroplasticity, 37 fractured porous media, 38 and fracture propagation. 39,40 For nonlinear problems, one combines a linearization technique, eg, the L-scheme, 41,42 with the splitting algorithm; the convergence of the resulting scheme can be proved rigorously. 32,33 Finally, we would like to mention some valuable variants of the fixed-stress splitting scheme: the multirate fixed-stress method, 43 the multiscale fixed-stress method, 29 and the parallel-in-time fixed-stress method. 44 This paper is structured as follows. The notation, the discretization, and the fixed-stress splitting scheme are presented in Section 2. The theoretical analysis of the convergence and the optimization technique are the subject of Section 3. In Section 4, numerical experiments that test the optimization technique are presented. Finally, conclusions are given in Section 5.

THE NUMERICAL SCHEME FOR SOLVING BIOT'S MODEL
In this paper, we use common notations in functional analysis. Let Ω ⊂ ℝ d be a Lipschitz domain where d is the spatial dimension. The space L 2 (Ω) is the Hilbert space of Lebesgue-measurable, square-integrable functions on Ω, and H 1 (Ω) is the Hilbert space of functions in L 2 (Ω) with derivatives (in the weak sense) in L 2 (Ω). The inner product and its associated norm in L 2 (Ω) are denoted by ⟨·, ·⟩ and || · ||, respectively, and || · || H 1 (Ω) is the standard H 1 (Ω)-norm. Vectors and tensors are written bold, and, sometimes, the scalar product and the norm will be taken for vectors and tensors. Vectorial functions are written bold-italic. T will denote the final time.
Biot equations (1) and (2) are solved in the domain Ω × (0, T) together with (for simplicity) homogeneous Dirichlet boundary conditions and a given initial condition. In time, the backward Euler method is applied with a constant time-step size ∶= T N , N ∈ ℕ. Throughout this work, the index n will refer to the time level. For the spatial discretization, a two-field Galerkin finite element formulation is considered, and two generic discrete spaces V h and Q h , associated with displacements and pressures, are introduced. Later, we require V h ×Q h to be inf-sup stable with respect to the divergence operator; the most prominent inf-sup-stable example is the Taylor-Hood element, ie, P2-P1 for displacement and pressure. 45 Nevertheless, the analysis below can be extended without difficulties to a three-field formulation as, for example, in the works of Phillips and Wheeler, 3 Both et al, 5 and Berger et al. 6 In this way, the fully discrete, weak problem reads: let n ≥ 1 and assume (u n−1 h , are obtained by using the initial condition. The fixed-stress splitting scheme 5,21,23,28 is now introduced. Denote by i the iteration index. Iterate until convergence. For i ≥ 1, given a stabilization parameter The initial guess for the iterations is chosen to be the solution at the last time step, ie, Notice that the mechanics and flow problems decouple, allowing for the use of separate simulators for both subproblems.

CONVERGENCE ANALYSIS AND OPTIMIZATION
In this section, the convergence of the scheme (5)-(6) is analyzed. We are particularly interested in finding an optimal stabilization parameter L, in the sense that the scheme requires the least amount of iterations, ie, has the smallest possible convergence rate. Before we proceed with the main result, we need some preliminaries. Definition 1. The mathematical bulk modulus, K ⋆ dr > 0, is defined as the largest constant such that By the Cauchy-Schwarz inequality, we get that the physical drained bulk modulus K dr = 2 d + is a lower bound for K ⋆ dr . However, for effectively lower-dimensional situations, eg, a one-dimensional-like compression, d can be replaced by a value closer to 1. Lemma 1 below guarantees an upper bound for K ⋆ dr . Nevertheless, there is a strong indication (based on numerical experiments; see, eg, Section 4 and the work of Both and Köcher 27 We remark that the exact value, depending on the physical situation, can be computed as a generalized eigenvalue.
Throughout this paper, we make use of the following two assumptions. Assumption 1. The constants , , , and are strictly positive, the constants 1∕M and are nonnegative, and the vector g is constant.
Proof. Consider Corollary 1. Let the continuous linear function The first statement of Corollary 1 is a characterization of an inf-sup-stable discretization Assumption 2, with inf-sup constant . Hence, the second statement of Corollary 1 holds; there exists a linear function Additionally, the following chain of inequalities holds true: where the first inequality follows from Young's inequality with C depending only on the Lamé parameters, and the second inequality results from the operator norm, ie, We obtain our desired inequality, as follows:  (5) and (6), and (u n h , p n h ) is a solution to (3) and (4). The fixed-stress splitting scheme (5)- (6) converges linearly for any L ≥ through the error inequalities where C Ω is the Poincaré constant and is the constant from (8).
Proof. Subtract (5) and (6) from (3) and (4), respectively, to obtain the error equations holding for all To prove (11), test (12)(i) with v h = e n,i u , and apply the Cauchy-Schwarz inequality and Young's inequality to the pressure term to obtain We now get (11) by applying (7). In order to prove (10), test (12) Using now Equation (12) By applying Young's inequality in (14), we obtain that, for any > 0, there holds To take care of the last term in (15), consider Equation (12) By using (7), (16) implies Inserting (17) into (16) yields By rearranging terms and inserting (18) into (15), we immediately get Using that L ≥ 2 K ⋆ dr and the Poincaré inequality, we obtain from the above The result, (19), already implies that we have convergence of the scheme. In previous works, particularly that of Both et al 5 (where the proof so far is very similar), the conclusion at this point is that L = 2 2K ⋆ dr is the optimal parameter. However, this does not consider the influence of the first term in (19). By Lemma 1, we get that there exists v h ∈ V h such that e n,i p = ∇ · v h in a weak sense and By testing now (12) From (20) and (21) and the Cauchy-Schwarz inequality, we immediately obtain which, together with (19), implies This gives the following rate of convergence, for ∈ (0, 2] and L ≥ Remark 2. Assumptions 1 and 2 are valid in various relevant physical situations. Therefore, our analysis has a wide range of applications. One can easily extend the result to heterogeneous media, ie, = (x) as long as is bounded from below by m ≥ 0. Moreover, any of the other parameters can be chosen spatially dependent as long as they are bounded from below by appropriate constants satisfying Assumption 1.

Optimality
Consider the rate obtained in (9). As rate(L, ) is an increasing function of L, it follows that, for all ∈ (0, 2], its minimum is obtained at L = 2 K ⋆ dr , giving the rate Minimizing (23) with respect to corresponds to maximizing .
It is easily seen that the maximum of (A − B) is attained at = A 2B . Therefore, the minimizer of rate( ) is since A ≥ 2B. This suggests that the theoretical optimal choice of L is .

(25)
Remark 3 (Consequence for low-compressible fluids and low-permeable porous media). Previous convergence results in the literature for the fixed-stress splitting scheme have not predicted or guaranteed any robust convergence in the limit cases M → ∞ and → 0 (for a fixed time-step size ). Now, by Theorem 1, for inf-sup-stable discretizations, robust convergence of the fixed-stress splitting scheme is guaranteed, even in the limit case. This was studied numerically in the work of Storvik et al. 31 Convergence was showed to be robust with respect to material parameters for P2-P1 elements and deteriorating for P1-P1.

Brute-force optimization of the stabilization parameter
The rate obtained in Theorem 1 is not necessarily sharp, and it is rather viewed as theoretical evidence that the optimal stabilization parameter resides in the interval [

NUMERICAL EXAMPLES
In this section, we demonstrate the effectiveness of the proposed brute-force method for optimizing the stabilization parameter for the fixed-stress splitting scheme. In particular, we show for several numerical test cases that the optimal stabilization parameter is close to being mesh independent and that the method for choosing it optimally, as described in For the implementation of the numerical examples, we use modules from the DUNE project, 47 particularly dune-functions. 48,49 If not mentioned otherwise, the inf-sup-stable Taylor-Hood pair P2-P1 is utilized as spatial discretization. As stopping criteria, we have applied relative L 2 -norms for the pressure, ie, iterations stop when h ||, consistent with Theorem 1. Constant material and fluid parameters are applied and given for each individual test case.

Notations
During the numerical experiments, we apply some specific choices of stabilization parameters several times. Therefore, we give them names here. Recall the definition of the physical drained bulk modulus K dr = 2 d + . The original stabilization parameter will be called the physical one due to the fixed-stress splitting scheme's physical origin, ie, L phys = 2 K dr . The other classical choice of stabilization parameter will be named after Mikelić and Wheeler due to their paper, 24  . The stabilization parameter obtained by the brute-force method described in Section 3.2 will be called L opt . The final parameter is the one that is proposed to be the smallest possible choice in Section 3.1, ie, L min = 2 4 +2 (see Table 1).

Dependence on boundary conditions-the unit square
We consider two test cases differing solely in the applied boundary conditions. Common for both, the domain is the unit square discretized by structured triangles, and the constant material parameters from Table 2 are considered. Moreover, we employ source terms corresponding to the analytical solution of the continuous problem (1)- (2). The pressure, p, is scaled by p ref = 10 11 Pa in order to balance the magnitude of the mechanical and fluid stresses for the chosen physical parameters. Regarding the different sets of boundary conditions, we consider the following.
• BC1: homogeneous Dirichlet data on the entire boundary for displacement and pressure.
• BC2: homogeneous Dirichlet data for the pressure; homogeneous Neumann data on top in the mechanics equation and homogeneous Dirichlet data everywhere else for the displacement.
Solutions after 10 time steps using a mesh size of h = 1∕128 are displayed in Figures 1 and 2.
To motivate the brute-force approach from Section 3.2, the performance of the fixed-stress splitting scheme has been measured for a variety of stabilization parameters and mesh sizes (see Figure 3). We observe that the numbers of iterations vary significantly for different stabilization parameters but that the optimal choice is within our proposed interval [L min , L phys ). Additionally, for fixed stabilization parameters, we observe that the numbers of iterations are close to constant with respect to the mesh size. Now, we test the brute-force approach of Section 3.2. In order to calculate L opt , we start by applying the fixed-stress splitting scheme for 11 equidistant stabilization parameters in [L min , L phys ] while only computing one time step for a mesh   Table 2. The largest value of  corresponds to L min , whereas the smallest value of  corresponds to L phys . Recall that L opt is calculated using only one time step, and therefore, there is a slight deviation between L opt and the actual optimal choice. A, BC1; B, BC2 size of h = 1∕16. Then, using the stabilization parameter that needed the least amount of iterations to converge, we apply the fixed-stress splitting scheme for the full problem using a mesh size of h = 1∕512. In Figure 3, the average numbers of iterations over 10 time steps are displayed for this "optimal" stabilization parameter, for the two classical choices L phys and L MW , and for the stabilization parameter that we consider to be the smallest possible choice, ie, L min . We see that the optimized stabilization parameter requires the least amount of iterations for both boundary conditions. It is also worth noticing that the optimal choice differs considerably for the two sets of boundary conditions.

Dependence on Poisson's ratio-L-shaped domain
To further analyze the proposed brute-force optimization of the stabilization parameter for the fixed-stress splitting scheme, we test it on an L-shaped domain as well. The L-shaped domain is considered as a subdomain of the unit square domain where the top-right quarter square has been removed, ie, L = [0, 1] 2 ∖(0.5, 1] 2 . The material and implementation parameters from Table 3 are applied, whereas the right-hand side is the same as for the unit square test case. Zero Dirichlet boundary conditions are applied everywhere, but at the top boundary ([0, 0.5] × {1}) for the mechanics equation where zero Neumann conditions are considered. A solution to this problem after 10 time steps with = 0 and mesh size 1∕h = 128 is given in Figure 4.
Given Young's modulus E and Poisson's ratio , the corresponding Lamé parameters have been determined by   Again, as for the unit square test case, we test the brute-force optimization technique that is described in Section 3.2, but now for three different Poisson's ratios. In Figure 5, the fixed-stress splitting scheme is applied to a variety of mesh sizes and with a variety of stabilization parameters to three problems with different Poisson's ratios. There are several key observations to make. First, the scheme is close to being mesh independent for all mesh sizes, stabilization parameters, and Poisson's ratios. Second, we see that the optimal stabilization parameter is in the proposed interval [L min , L phys ) for all Poisson's ratios and all mesh sizes. The final observation is that when the Poisson's ratio increases, the choice of stabilization parameter becomes less important. This is due to the fact that an increase in the Poisson's ratio can be seen as an effective decrease in the coupling strength.
To calculate the optimal stabilization parameter, we follow the recipe of Section 3.2. We apply 11 equidistant stabilization parameters in the interval [L min , L phys ] for the fixed-stress splitting scheme on a coarse mesh (1∕h = 16) for only one time step. Counting the numbers of iterations it takes to reach convergence, we choose the parameter that corresponds to the smallest number and use this for the finer mesh (1∕h = 512) and more time steps (10). We see that the parameter that is the optimal choice for the coarse mesh is also the optimal one for the finer mesh for all Poisson's ratios.

Mandel's problem
Here, we consider Mandel's problem, a relevant two-dimensional problem with a known analytical solution that is often used as a benchmark problem for discretizations. The analytical solution is derived in the works of Coussy 1 and Abousleiman et al, 50 and its expressions for pressure and displacement are given by where n , n ∈ ℕ, correspond to the positive solutions of the equation and u , F, B, c f , and a are input parameters, partially depending on the physical problem parameters. Here, we apply the values listed in Table 4. For a thorough explanation of the problem and the coefficients in (27)- (29), we refer to the works of Coussy 1 and Phillips and Wheeler. 3 We consider the domain, Ω = (0, 100) × (0, 10), discretized by a regular triangular mesh. An equidistant partition of the time interval is applied with time-step size = 10 from t 0 = 0 to T = 100. Initial conditions are inherited from the analytic solutions (27)- (29). As boundary conditions, we apply exact Dirichlet boundary conditions for the normal displacement on the top, left, and bottom boundaries. For pressure, we apply homogeneous boundary conditions on the right boundary. On the remaining boundaries, homogeneous natural boundary conditions are applied. The tolerance r is set to 10 −6 . The solution after 10 time steps with 80 vertical and horizontal nodes is displayed in Figure 6.
Similar to the unit square and L-shaped domain test cases, we test the mesh independence and the brute-force optimization technique for Mandel's problem. This time, the parameters from Table 4 are applied. In Figure 7, the mesh dependence of the fixed-stress splitting scheme is tested, and it is clear that the performance of the scheme is independent of this choice. At the same time, we confirm that the optimal stabilization parameters actually are in the proposed interval [L min , L phys ).   Table 4. The largest value of  corresponds to L min , whereas the smallest value of  corresponds to L phys To calculate the optimal stabilization parameter, we have applied the optimization technique of Section 3.2. First, the fixed-stress splitting scheme is applied for one time step using a coarse mesh with 10 horizontal and 10 vertical nodes for 11 different stabilization parameters in the interval [L min , L phys ]. Choosing the parameter that yields the lowest number of iterations, we apply the scheme for finer meshes and count the number of iterations. As for the other test cases, we see that the optimal parameter indeed is optimal. Moreover, a poor choice of stabilization parameter can result in a huge number of iterations.

3D footing problem
The numerical section is concluded with a three-dimensional example, ie, a footing problem similar to a test case studied in the work of Adler et al. 46 We consider a unit cube subject to normal compression, ramped in time n (t) = t·10 10 N·m 2 ∕s, applied to a part of the top boundary Γ N ∶= [0.25, 0.75] × [0.25, 0.75] × {1}. The bottom is fixed in all directions, and the remaining boundary is considered to be stress free. A no-flow boundary condition is applied at the compression zone Γ N , and zero pressure is enforced on the remaining boundary. Furthermore, zero body forces are applied. The medium is considered isotropic with the same material parameters as used in Section 4.2 (cf Table 2). For the numerical discretization, we consider a set of four meshes with mesh size h ∈ {1∕8, 1∕16, 1∕32, 1∕64} and employ the inf-sup-stable MINI element. 51 The simulation result for the final time step is visualized in Figure 8.
Due to high computational cost, optimizing the stabilization parameter of the fixed-stress splitting becomes tedious for fine meshes in 3D. Motivated by the previous results, the optimal stabilization parameter is assumed to be nearly mesh independent. This allows for a brute-force search for the optimal, practical stabilization parameter utilizing the coarsest grid (cf Section 3.2). For validation of the optimization strategy, the performance of the splitting scheme is measured in the range [L min , L phys ] suggested by Theorem 1; for the finest mesh, we restrict the validation only to a neighborhood of the optimized stabilization parameter. The performance measured in terms of the number of iterations is presented in Figure 9. A large contrast in the performance can be observed for different stabilization parameters, emphasizing the need for a suitable stabilization parameter. Finally, as before, we observe that, indeed, the optimal, practical stabilization parameter is only slightly mesh dependent; it is close to the physical bulk modulus K dr = 2 d + . All in all, the brute-force  Table 2. The largest value of  corresponds to L min , whereas the smallest value of  corresponds to L phys search strategy from Section 3.2 has, again, been confirmed to be a suitable method to obtain a satisfactory stabilization parameter for finer meshes.

CONCLUSIONS
In this work, we have considered the quasi-static, linear Biot model for poroelasticity and studied theoretically and numerically the convergence of the fixed-stress splitting scheme. An improved convergence result has been proved, indicating the nontrivial dependence of the optimal stabilization parameters on not only mechanical properties but also fluid flow properties and discretization properties. We observe numerically that the fixed-stress splitting scheme is close to being mesh independent and determine a novel domain in which the optimal stabilization/tuning parameter is found, ie, [ . On the basis of these observations, we propose a brute-force method with low cost for choosing the optimal stabilization parameter, ie, the parameter that corresponds to the smallest amount of fixed-stress iterations. Through numerical experiments, we have showed that this optimization method results in a much faster fixed-stress splitting scheme than those obtained by choosing the classical stabilization parameters L =