Non-Iterative Partitioned Methods for Uncoupling Evolutionary Groundwater–Surface Water Flows

: We present an overview of a modern, efficient approach for uncoupling groundwater–surface water ﬂows governed by the fully evolutionary Stokes–Darcy equations. Referred to as non-iterative partitioned methods, these algorithms treat the coupling terms explicitly and at each time level require only one Stokes and one Darcy sub-physics solve, thus taking advantage of existing solvers optimized for each sub-ﬂow. This strategy often results in a time-step condition for stability. Furthermore, small problem parameters, speciﬁcally those related to the physical characteristics of the porous media domain, can render certain time-step conditions impractical. Despite these obstacles, researchers have made signiﬁcant progress towards efﬁcient, stable, and accurate partitioned methods. Herein, we provide a comprehensive survey and comparison of recent developments utilizing these non-iterative numerical schemes.


Introduction
Access to the clean freshwater is absolutely imperative for the continued survival of humankind.As a necessity for our agricultural, industrial and domestic practices, water constitutes an integral part of all civilizations.However, only 2.5% of the water present on Earth is freshwater, and the majority of this amount is either frozen or inaccessible.Furthermore, 96% of accessible freshwater comes from aquifers underground.Because of the scarcity of this resource, we must prioritize the protection and conservation of groundwater sources.Too often, human and natural processes threaten groundwater quality, sometimes irreversibly.For example, in hydro-fracturing, companies inject a mixture of water with sand and chemicals at high pressure into a well to create fractures to allow for the collection of shale gas.Companies do not recover the majority of the chemicals in this mixture and many fear that eventually these pollutants will leave the well to contaminate the local groundwater supply.Pesticide application in agriculture can have devastating effects on surrounding freshwater resources due to chemical run-off into nearby rivers, lakes, and streams, and seepage deep into the soil.Furthermore, many storage facilities for radioactive materials exist underground for both safety and convenience.Over time, as storage containers become compromised, nuclear waste can migrate into nearby freshwater aquifers.Even natural processes may result in contaminated freshwater, as evident in the devastation of forests growing above coastal aquifers from salt-water intrusion.
Tracking these contaminants necessitates accurate numerical models for this coupled flow.Scientists have thoroughly studied the individual groundwater and surface water flows (see, for example, Pinder and Celia [1], Watson and Burnett [2], or Bear [3] for an extensive study on subsurface flows, and Kundu, Cohan and Dowling [4] for surface water flows).As a result, many accurate and efficient solvers for the independent flow processes exist.Modeling the interaction of groundwater and surface water, however, presents additional difficulties as we must preserve the physical processes in each sub-flow while accurately describing the activity occurring along the interface.
An attractive and practical strategy, which is the main focus of this survey paper, is to make use of the existing solvers for separate fluid and porous media flow by investigating methods that uncouple the flow equations in time so that the individual flow problems may be solved separately.Called partitioned methods (also domain decomposition methods), these methods allow us to utilize, in a black-box manner, solvers already optimized for the separate flow problems.It is important that the partitioned methods maintain stability and accuracy along the interface where the two flows meet.In addition, potentially small physical parameters create an additional challenge for stability.We are concerned with methods that are efficient for the time-dependent models, in particular, the ones that are stable over long-time intervals, since groundwater moves slowly and numerical simulations may span long-time periods.Along these same lines, we want methods that converge within a reasonable amount of time to be of practical use, making higher-order convergence a desirable property.
In recent years, several researchers have made substantial progress in the development of non-iterative, partitioned methods applied to the evolutionary groundwater-surface water flow problem.Based on an implicit discretization in time for each subproblem, these methods, however, make use of results from previous time steps to predict the values on the interface at the current time step, thus requiring only one solve for groundwater and one solve for surface water flow at each time level (thus non-iterative).In this work, we will review and discuss several such methods so as to illustrate the current status of this important problem.The modeling of coupled fluid-porous media flow begins with the coupling of the Stokes or Navier-Stokes equations describing the flow in the fluid region, along with the Darcy or Brinkman equations for the flow in the aquifer, or porous media region containing the groundwater.This survey focuses on the Stokes-Darcy coupling that is suitable for slow moving flows over large domains.
Studies on the continuum surface water-groundwater model have been performed in [5][6][7][8].The literature on numerical analysis of methods for the coupled Stokes-Darcy problem has grown extensively since [9,10] (see, for example, [11][12][13] for analysis of the steady-state problem).There exist many effective and efficient domain decomposition techniques for decoupling the Stokes-Darcy system in the stationary case [14][15][16][17][18][19][20][21][22][23][24].To solve the fully evolutionary Stokes-Darcy problem, one approach is monolithic discretization by an implicit scheme (see, e.g., [25,26]).These schemes can also be solved by an iterative domain decomposition method at each time step.In general, any decoupling technique for stationary Stokes-Darcy (many cited above) may be applied to find the solution at each time level in the time-dependent case.
Non-iterative partitioned methods, an alternative approach, are advantageous in that they allow uncoupling into only one (SPD) Stokes and one (SPD) Darcy system per time step.Mu and Zhu presented the first non-iterative partitioned scheme in [27], proposing employing Backward Euler discretization for each subproblem while treating the coupling term explicitly by Forward Euler.Layton, Tran and Trenchea revisited this method in [28], with an improved analysis showing long-time stability.In that work, the authors also developed and tested for efficiency a second first-order scheme, Backward Euler-Leap Frog.Following these methods, others proposed several other implicit-explicit (IMEX) methods of high order, such as Crank-Nicolson-Leap Frog [29], second-order backward-differentiation with Gear's extrapolation [30], and Adam-Moulton-Bashforth [30,31].Although these methods use explicit discretizations for the coupling terms, all are now known to be long-time stable and optimally convergent uniformly in time (possibly under a small time-step constraint).With the addition of suitable stabilization terms, it is possible to further enhance the stability property, for instance, a stabilized Crank-Nicolson-Leap Frog, developed in [32,33], requires no time-step restriction for the long-time stability and convergence.Another way for uncoupling groundwater-surface water systems is using splitting schemes.Unlike the aforementioned IMEX schemes that solve for separate sub-flows in parallel, splitting methods require sequential sub-problem solves at each time step.In [34], the authors proposed four first and second-order splitting schemes.Theoretical and numerical evidence provided therein suggests that these methods are stable for larger time steps than the first order IMEX schemes and, in particular, a good option in case of small physical parameters.Finally, asynchronous (aka, multiple-time-step, multi-rate) partitioned methods allow for different time steps in the two subregions, motivated by the observation that the flow in fluid region occurs with higher velocities compared to flow in porous media region.Such methods may be more efficient, as we may apply two different time steps to separately solve the fast and slow flows.Developed in [35,36], these asynchronous techniques utilize the Backward Euler-Forward Euler time discretization, with long-time stability acquired in the latter work.
We organize this paper as follows.Section 2 reviews the preliminaries of the Stokes-Darcy equation, including interface conditions, variational formulation and semi-discrete approximations.We briefly discuss the implicit time discretization, together with the iterative domain-decomposition approach.Section 3 focuses on first-order partitioned methods.We will survey several different approaches including first-order IMEX schemes and splitting schemes.We review high-order methods in Section 4 and asynchronous partitioned techniques in Section 5. Finally, we provide some conclusions and outlooks in Section 6.

The Stokes-Darcy Equation
Let the fluid region be denoted by Ω f and the porous media region by Ω p .Assume both domains are bounded and regular.Let I represent the interface between the two domains.We assume the time-dependent Stokes flow in Ω f and the time-dependent groundwater flow along with Darcy's law in Ω p .The Stokes-Darcy equation, describing the fluid velocity field u = u(x, t) and pressure p = p(x, t) on Ω f and the porous media hydraulic head φ = φ(x, t) on Ω p , can be written as follows: u(x, t) = 0, in ∂Ω f \I and φ(x, t) = 0, in ∂Ω p \I, + coupling conditions across I. (1) Here, f f denotes the body force in the fluid region, f p is the sink or source in the porous media region, ν > 0 is the kinematic viscosity of the fluid, S 0 is the specific mass storativity coefficient and K is the hydraulic conductivity tensor, assumed to be symmetric, positive definite with spectrum(K) ∈ [k min , k max ].
It is worth noting that values of S 0 and the smallest eigenvector k min of K can be very small (see Tables 1 and 2 for the values of S 0 and k min for different materials).As we shall see, this poses a major challenge in designing partitioned methods with good stability.Indeed, partitioning often induces time-step restrictions for long-time stability, which may become severe in the case of small system parameters.

Interface Conditions
To close the coupled problem formulation, a set of conditions has to be defined on the interface.Let n f /p denote the indicated, outward pointing, unit normal vector on I.The first two coupling conditions involve the conservation of mass and balance of forces on I: In addition, we need a tangential condition on the fluid region's velocity along the interface.In [5], Beavers and Joseph proposed the following slip-flow condition, expressing that slip velocity along I is proportional to the shear stresses along where α BJ is a dimensionless constant depending solely on the porous media properties and ranges from 0.01 to 5, g is the gravitational acceleration constant, and u p is the average velocity in the porous media region.The validity of Beavers-Joseph interface condition has been supported by abundant empirical evidence; however, one challenge in adopting this condition is that the bilinear form in the weak formulation is not coercive.Several simplifications have been considered.In [6], Saffman proposed a modification to the Beavers-Joseph coupling condition by dropping the porous media averaged velocity u p , based on observations that the term u p • τ i is negligible compared to the fluid velocity u • τ i .This simplified condition was mathematically justified in [39] and has been shown satisfactory for many fluid-porous media systems.Known as Beavers-Joseph-Saffman(-Jones) coupling condition, this is the third and final condition we use in this article: , on I for any τ i tangent vector on I.
For the analysis and numerical methods for Stokes-Darcy systems with Beavers-Joseph condition, we refer to [21,26,40].

Variational Formulation and Semi-Discrete Approximations Using Finite Element Method
We denote the L 2 (I) norm by • I and the L 2 (Ω f /p ) norms by • f /p , respectively, and the corresponding inner products are denoted by (•, •) f /p .In addition, define the H div (Ω f ) and the functional spaces and the bilinear forms Note that, setting v = u, ψ = φ and adding, the coupling terms exactly cancel out in the monolithic sum yielding the energy estimate for the coupled system.
To discretize the Stokes-Darcy problem in space by the finite element method (FEM), we select finite element spaces velocity: based on a conforming FEM triangulation with maximum triangle diameter denoted by "h".We do not assume mesh compatibility or interdomain continuity at the interface I between the FEM meshes in the two subdomains.The Stokes velocity-pressure FEM spaces are assumed to satisfy the usual discrete inf-sup condition for stability of the discrete pressure (see, e.g., [41]) The semi-discretization for the time-dependent Stokes-Darcy problem is as follows: find It is worth noting that the coupling between the Stokes and the Darcy sub-problems is exactly skew symmetric.

Fully-Discrete Approximations with Fully Implicit Temporal Schemes
Let t n := n∆t and w n := w(x, t n ) for any function w(x, t).For V being a Banach space with norm • V , we denote the following discrete norms where N = T/∆t and T can be ∞.
, on the other hand, depends on the time step ∆t.However, for functions smooth in time, this norm converges to the continuous norm • L 2 (0,T;V) as ∆t → 0. In those cases, one can reasonably assume the uniform bound of | • | L 2 (0,T;V) , independent of ∆t.
The most natural time discretization for the Stokes-Darcy equation is perhaps the first-order backward Euler scheme, which, in combination with the aforementioned finite element Galerkin method for the spatial discretization, leads to the following fully implicit, coupled problem.
Stability and convergence analysis of this scheme were conducted in [25][26][27], for both Beavers-Joseph and Beavers-Joseph-Saffman-Jones interface conditions.Higher order fully implicit schemes, such as the Crank-Nicolson, can also be considered.In general, fully implicit methods possess superior stability compared to IMEX or splitting temporal schemes.The major concern here is that this approach must solve a coupled problem at each time level.Partitioning the coupled problem at each time step is possible, but involves an iterative procedure with additional cost.In principle, any decoupled methods developed for the stationary model can be used in iteration at each time level.

First Order Partitioned Schemes
An attractive alternative to fully implicit, fully coupled discretization is exploiting information obtained in previous time steps to construct a non-iterative uncoupling scheme, which only need a single Stokes solve and a single Darcy solve at each time step.This approach allows the use of legacy subproblems' codes and obviously requires less programming effort (compared to solving coupled Stokes-Darcy system directly) as well as less computation cost (compared to iterative domain decomposition approach).As the interface values are obtained in an explicit manner, the main challenge here is how to obtain optimal accuracy and good stability properties.Many non-iterative partitioned methods have been developed in the literature recently [27][28][29][30][31][32][33][34][35][36]42], whose stability and accuracy have been proved (over a long time or without time-step condition) and numerically tested.Several of them maintain good performance even in the case of small parameters.The rest of this paper represents an overview of these developments.Our discussion will be divided into three parts: in Section 3, we survey first order schemes; in Section 4, high order schemes will be discussed; Section 5 is devoted to asynchronous partitioned methods.Unless otherwise stated, C denotes a generic positive constant whose value may be different from place to place but which is independent of mesh size, time step, and final time.For all the methods surveyed, approximations are needed at the first few (one or more) time steps to begin, and we always assume these are computed to sufficient accuracy.

Backward Euler-Forward Euler
The first non-iterative uncoupling scheme is Backward Euler-Forward Euler (BEFE), proposed by Mu and Zhu in [27] (and referred to as DBES therein).This method applies Backward Euler discretization for the subproblems and treats the coupling terms by explicit Forward Euler: A stability analysis for BEFE was given in [27].These results only apply for bounded time intervals [0, T] with T < ∞, as the estimates include e cT multipliers and thus grow exponentially with T. The long-time stability of BEFE was established in [28].An important feature of this proof, also of other long-time results coming next, is that no form of Gronwall's inequality was used.This result can be stated as follows.
Then, the following hold: BEFE is first order in time.A convergence analysis of this scheme can be found in [27], with a very recent improvement in [43].

Backward Euler-Leap Frog
Backward Euler-Leap Frog (BELF) is another IMEX scheme, first proposed in [28].This method is a combination of the three level implicit method with the coupling terms treated by the explicit Leap-Frog method.Approximations are needed at the first two time steps to begin.The stability region of the usual Leap-Frog time discretization for y = λy is exactly the interval of the imaginary axis −1 ≤ Im(∆tλ) ≤ +1.Thus, LF is unstable for every problem except for ones that are exactly skew symmetric such as the coupling herein.
The Backward Euler-Leap Frog scheme can be formulated as follows: As with any explicit scheme, BELF inherits a time-step restriction for the stability.The following long-time stability result was established in [28].Proposition 2 (Long-time stability of BELF, [28]).Consider the scheme BELF.Assume that the following time-step condition is satisfied ∆t min{ νk min , S 0 νk min , νk 2  min , S 0 ν 2 k min }; then, BELF possesses the same stability properties as those for BEFE in Proposition 1.More precisely, It was also proved that BELF achieves the optimal convergence rate uniformly in time, as shown below.
Proposition 3 (Error estimate of BELF, [28]).Consider the scheme BELF.Assume the following time-step condition is satisfied ∆t min{ νk min , S 0 νk min , νk then the solution of BELF satisfies the uniform in time error estimates Numerical tests illustrating the theoretical stability and convergence properties of BEFE and BELF were presented in [28].In particular, the stability of these two methods is compared to that of the fully implicit method in the case of small k min for a Stokes-Darcy flow on Ω f = (0, 1) × (1, 2) and Ω p = (0, 1) × (0, 1) with the interface I = (0, 1) × {1}.Given the source terms f f ≡ 0, f p ≡ 0, the initial condition set all the physical parameters (except for ν and k min ) to 1. Letting h = 1 10 , ν = 1 10 , the evolution of the energy p with k min = 10 −6 is shown in Figure 1.Since the true solution decays as t → ∞, any growth in E n indicates instability.The plot reveals that while not unconditionally stable like the fully implicit method, BEFE and BELF only require mild constraints on ∆t for their stability.Indeed, BELF is already stable for ∆t     In summary, the first-order IMEX methods BEFE and BELF allow a parallel, non-iterative uncoupling of the Stokes-Darcy system at each time step and, on the other hand, enjoy the desirable strong stability and convergence properties.One disadvantage, as shown in their time-step restrictions, is that they may become highly unstable when one of the parameters S 0 and k min is small.In the next subsection, this constraint is relaxed with another type of first-order partitioned methods, splitting schemes.

First Order Sequential Splitting Schemes
In [34], several methods for non-iterative, sub-physics uncoupling the evolutionary Stokes-Darcy problem were proposed, using ideas from splitting methods.The estimates and tests therein suggest that these methods are stable for larger timesteps than the IMEX based partitioned methods BEFE and BELF, and in particular, a very good option when either k min or S 0 (but not both) is small.Here, the Stokes and Darcy systems are uncoupled, but, unlike the aforementioned IMEX schemes, sequentially solved.
In the first Backward Euler time-split (BEsplit1) scheme, the coupling term in the φ equation is evaluated at the newly computed value u n+1 h so we compute Algorithm 4 Backward Euler time-split 1 (BEsplit1) The long-time stability of this scheme can be stated as follows.
Proposition 4 (Long-time stability of BEsplit1, [34]).Consider the scheme BEsplit1.Assume the following time-step condition is satisfied: The second Backward Euler time-split (BEsplit2) is the previous method in the opposite order, i.e., computing u n h → φ n+1 h → u n+1 h .The analysis in [34] revealed that control was needed for a term u n+1 h − u n h div .This led to the insertion of the grad-div stabilization term (∇ • (u n+1 h − u n h )/∆t, ∇ • v h ) acting on the time discretization of u t .This term is exactly zero for the continuous problem so it does not increase the method's consistency error.
Algorithm 5 Backward Euler time-split 2 (BEsplit2) The stability result of BEsplit2 is presented below.
Propositions 4 and 5 impose two slightly different conditions on ∆t, both of which are mild when one of S 0 and k min is small.In those cases, BEsplit1 and BEsplit2 are preferable choices to first-order IMEX schemes, with the small price of solving the uncoupled subproblems sequentially, instead of in parallel.However, it is worth remarking that these methods may become unstable if both parameters are small.Finally, BEsplit1 and BEsplit2 can be shown to be optimally convergent under the same time-step conditions for stability.For a thorough analysis, we refer to [44].
Two other splitting methods were proposed in [34], whose details are omitted here for brevity.SDsplit is a first order scheme, long-time stable under the condition ∆t min {S 0 , k min } h, and thus seems less favorable than BEsplit1 and BEsplit2 in theory.CNsplit, on the other hand, is stable with ∆t √ S 0 h and a very good option in case of small k min .This scheme is second order.Numerical tests checking and comparing the largest time step for which the four methods are stable over long-time intervals were also performed in [34].Taking the initial condition as in (3), the body sources to be 0, the system parameters (except k min and S 0 ) to be 1.0 and mesh size h = 1  10 , the authors computed the system energy E N at final time T = 10 with different time-step sizes.Since the true solution decays as t → ∞, large E N indicates instability.The performance of presented splitting methods was plotted for three cases: (i) O(1) k min and small S 0 , (ii) small k min and O(1) S 0 , and (iii) small k min and small S 0 (see Figures 2 and 3).These plots show that for small parameter k min or S 0 , BEsplit1 and BEsplit2 are stable for large time steps.The performance of SDsplit is close to those of BEsplit1 and BEsplit2, suggesting that its theoretical condition was not optimized.These three first order splitting methods display superior stability to IMEX methods in our previous tests.The second order CNsplit in general requires a much smaller time step, but still possesses strong stability in the case of small k min and large S 0 .

High Order Partitioned Schemes
In this section, we review recent developments in high-order partitioned schemes for uncoupling the Stokes-Darcy system.Thus far, most of the proposed methods are of IMEX type (except for CNsplit, which was mentioned above).The IMEX schemes we will discuss in detail here include Crank-Nicolson-Leap Frog [29,32,33], second-order backward-differentiation with Gear's extrapolation [30], and Adam-Moulton-Bashforth [30,31].

Crank-Nicolson-Leap Frog
The Crank-Nicolson-Leap Frog (CNLF) method is a second-order scheme that employs the implicit Crank-Nicolson discretization of subdomain terms and treats the interface terms explicitly with Leap Frog.CNLF was developed for uncoupling systems of evolutionary equations in [45][46][47].This method was first applied to uncouple Stokes-Darcy system in [29].
CNLF is a three-level method.The first terms, (u 0 h , p 0 h , φ 0 h ), arise from the initial conditions of the problem.To obtain (u 1 h , p 1 h , φ 1 h ), one must use another numerical method.Note that approximations in this first step will affect the overall convergence rate of the method, and as usual, we assume them to be sufficiently accurate.This scheme can be stated as follows: CNLF provably possesses a strong stability and convergence properties, as shown in [29] and presented below.Specifically, the time-step condition for CNLF does not depend on k min , making this method a very good choice for fluid-porous media systems with small k min .Proposition 6 (Long-time stability of CNLF, [29]).Consider the scheme CNLF.Assume the following time-step condition is satisfied: Proposition 7 (Error estimate of CNLF, [29]).Consider the scheme CNLF.Assume the following time-step condition is satisfied: ∆t max{min{h 2 , S 0 }, min{h, S 0 h}}, as in Proposition 6.If the solution of the Stokes-Darcy problem (1) is regular in the sense that and the discrete norms of u, p and φ in L 2 (0, ∞; While independent of k min , the conditional stability of CNLF may still be restrictive when faced with small S 0 .To tackle this difficulty, the authors of [32,33] proposed a strategy to improve the stability property of CNLF by adding appropriate stabilization terms to both the Stokes as well as the groundwater flow equation.
The resulting numerical scheme, denoted CNLFstab, is unconditionally, uniformly in time stable, as well as second-order convergent.More specifically, it was shown in [33] that the added stabilization terms,

∇ •
eliminate the time-step restriction without affecting the second-order accuracy of the method.Indeed, among all the methods we have discussed so far, CNLFstab exhibits the best stability and convergence properties.Here, C f ,p is a constant satisfying In the special case of a flat interface I, with Ω f and Ω p being arbitrary domains, C f ,p equals 1 (see [48, Lemmas 3.1 and 3.2]).We make the above discussion rigorous with the following results.
Proposition 9 (Error estimate of CNLFstab, [33]).Consider the scheme CNLFstab.Assume u, p and φ satisfy the same regularity condition as in Proposition 7, then the solution of CNLFstab satisfies the unconditional, uniform in time error estimate The stability and convergence properties of CNLF and CNLFstab were illustrated and compared via numerical tests in [33,44].Set the body sources to be 0, the system parameters (except k min and S 0 ) to be 1.0 and mesh size h = 1 10 , the authors computed the system energy E N at final time T = 10 with different time-step sizes.Since the true solution decays as t → ∞, large E N indicates instability.The performance of CNLF and CNLFstab was plotted for four cases: (i) small k min and small S 0 , (ii) small S 0 , (iii) small k min , and (iv) k min = S 0 = 1.0.These plots show that, in all situations, CNLFstab is stable for large time steps, regardless of the size of k min and S 0 .This is a vast improvement over CNLF, which, for ∆t ≤ 1/80, was only stable when k min = S 0 = 1.0, and unstable in all other cases (see Figure 4).

Second-Order Backward-Differentiation Formula with Gear's Extrapolation
In [30], the authors introduced a second-order scheme, which discretizes in time via a second-order backward-differentiation formula (BDF2), and treats the interface term via a second-order explicit Gear's extrapolation formula.
Algorithm 8 Second-order backward-differentiation (BDF2) BDF2 provably possesses a strong stability and convergence properties, as stated below.Specifically, an inspection of the analysis argument in [30] shows that the time-step condition for BDF2 does not depend on S 0 , making this method a very good choice for fluid-porous media systems with small S 0 (and moderate k min ).Proposition 10 (Long-time stability of BDF2, [30]).Consider the scheme BDF2.Assume the following time-step restriction is satisfied: We note that the solution of BDF2 was also proved to be long-time stable in H 1 norm, i.e., ∇u n h 2 f + ∇φ n h 2 p ≤ C, ∀n ≥ 0, in [30], with a more restrictive time-step condition.A strong convergence property of BDF2 is stated below.
Proposition 11 (Error estimates of BDF2, [30]).Consider the scheme BDF2.Assume ∆t is sufficiently small (independent of mesh size and final time).If the solution of the Stokes-Darcy problem (1) is long-time regular in the sense that then the solution of BDF2 satisfies the uniform in time error estimates This result showed not only the error estimate of the velocity with respect to the L 2 norm, but also with respect to the H 1 norm, as well as the error estimate of the pressure.The two latter estimates are not second order in time; however, these were also observed in the numerical experiments [30].Finally, the authors suggested that the stabilization terms may be added to the Stokes and Darcy solves correspondingly, with parameters γ f , γ p ≥ 0. While the analysis does not take advantage of the stabilization term, the numerical experiments demonstrate the benefit of this strategy in the sense that the presence of the stabilization term relaxes the time-step restriction.

Adam-Moulton-Bashforth
The second-order Adam-Moulton-Bashforth method (AMB2), studied in [30], combines the second-order implicit Adams-Moulton treatment of the symmetric terms and the second-order explicit Adams-Bashforth treatment of the interface term.
Algorithm 9 Second-order Adam-Moulton-Bashforth method (AMB2) The third-order Adams-Moulton-Bashforth method (AMB3), studied in [31], is a combination of the third-order explicit Adams-Bashforth treatment for the coupling term and the third-order Adams-Moulton method for the remaining terms.To the best of our knowledge, this has been the only third order IMEX scheme that was applied to uncouple the Stokes-Darcy equation so far.
Algorithm 10 Third-order Adam-Moulton-Bashforth method (AMB3) These methods proved to be long-time stable under small time-step restrictions.The explicit dependence of these conditions on the system parameters can be elaborated from [30,31] as follows.
Similar to CNLF and BDF2, appropriate stabilization terms can be added for both algorithms AMB2 and AMB3, which were numerically shown to significantly relax the time-step constraint (see [30,31]).

Asynchronous Schemes
In surface water-groundwater models, the flow in fluid regions is often associated with higher velocities, compared to flow in porous media region.In such cases, it may be desirable to apply an asynchronous scheme (aka, multiple-time-step scheme, multi-rate scheme), which computes fast solutions using a small time step and consider a larger time step for slow solutions.The first partitioned scheme that allows different time steps in the fluid and porous region for the nonstationary Stokes-Darcy problem was probably proposed and analyzed in [35].In that work, the decoupling is based on lagging the interfacial coupling terms following the BEFE method; thus, we will refer to the scheme as asynchronous BEFE or BEFE-as1.Let ∆s be the (small) time-step size in the fluid region Ω f and ∆t be the (large) time-step size in the porous region Ω p such that ∆t = r∆s.In addition, define n k := kr and let N = T/∆s, the number of small time step, and M = T/∆t = N/r, the number of large time step.The algorithm in [35] reads as follows.
Algorithm 11 Asynchronous Backward Euler-Forward Euler 1 (BEFE-as1) For k = 0 to T ∆t − 1, do the following: Find φ n k+1 ∈ X h p such that for all ψ h ∈ X h p : The stability of BEFE-as1 was proved in [35] over a bounded interval.In particular, one has the following: Proposition 13 (Stability of BEFE-as1).Consider the scheme BEFE-as1.Let T > 0 be any fixed time.Assume the following condition on the small time step ∆s S 0 νk min .
The study on multi-rate schemes continues with [36], in which the second asynchronous strategy based on BEFE was proposed.This method in computing the hydraulic head φ, instead of using free flow velocity averaged over multiple previous steps as in BEFE-as1, only uses the free flow velocity value at the immediately previous time level.As such, the long-time stability was acquired, with the time-step restriction depending not only on the model parameters, but also including the ratio between the time steps applied in the free flow and porous medium domains.A remarkable property of this method is that it conserves mass across the interface, which does not seem possible with BEFE-as1.To be precise, we state the method as follows.
Algorithm 12 Asynchronous Backward Euler-Forward Euler 2 (BEFE-as2) For k = 0 to T ∆t − 1, follow the same procedure as in BEFE-as1, except for Step 2, where it is replaced by: 2. set U n k := u n k+1 h .The stability results of BEFE-as2 can be established as follows.
Proposition 14 (Long-time stability of BEFE-as2, [36]).Consider the scheme BEFE-as2.Assume following condition on the time-step condition is satisfied: For the error analysis and numerical experiments illustrating the convergence rate and mass conservation, we refer to [36].

Conclusions
In solving the coupled Stokes-Darcy equations, the non-iterative partitioned approach is an attractive alternative to fully implicit, monolithic discretization (combining with either direct coupled problem solve or iterative domain decomposition methods).First, these uncoupling schemes allow the use of legacy sub-problems' codes, in which the spatial mesh, time step and numerical method may be optimized according to each subprocess.Second, this approach, by exploiting the interface information obtained in previous time steps, only needs a single Stokes solve and a single Darcy solve per time level (some splitting methods may require a couple solves), and is therefore very cost effective.Since the coupling terms are treated in an explicit manner, obtaining optimal accuracy and good stability properties is a major concern with these methods.In recent years, many proposed partitioned schemes have surpassed this challenge, with proven long-time stability and optimal convergence properties.Further improvements in efficiency include high-order discretizations, stabilization strategies, and asynchronous schemes.Table 3 summarizes and compares the stability and convergence properties for all numerical schemes surveyed herein.There are, however, several important questions that remain open, in our point of view.

1.
For the long-time stability, most of the current methods require time-step conditions sensitive to the sizes of system parameters.These conditions may become restrictive for the fluid-porous media coupling with small parameters, particularly S 0 and k min .See Table 3 for the stability dependence on problem parameters of all methods discussed here.While there are a few methods achieving some successes in this case, e.g., CNLFstab, in our opinion, long-time stable and accurate schemes in case of small parameters, in particular, when both S 0 and k min are small, are worth further study.

2.
Most of the existing methods have not accounted for the dependence of the time step on the domain size.This is an important problem, especially for domains with large aspect ratios.3.
To our knowledge, there are no adaptations of the asynchronous approach beyond first-order schemes.High order asynchronous methods are desirable and the next logical step.4.
The primary motivation for modeling the fully evolutionary Stokes-Darcy flow is transport contaminant tracking, a major concern in several modern environmental problems.However, the problem of coupling numerical methods for the time-dependent Stokes-Darcy equation, in particular non-iterative partitioned methods discussed herein, with the transport equation to simulate the path of chemicals is largely open.

1 30 ,
followed by BEFE at ∆t 1 50 .These conditions are much weaker than those predicted by the theory.

Figure 1 .
Figure 1.The evolution of system energy with k min = 10 −6 for different choices of time step [28].Copyright c 2013 Society for Industrial and Applied Mathematics.Reprinted with permission.All rights reserved.
K∇φ, ∇ψ) p , and c I (u, φ) = g 2 min , S 0 ν 2 k min }, as in Proposition 2. If the solution of the Stokes-Darcy problem (1) is long-time regular in the sense that u

min Third-Order in Time Method Type Stability Condition AMB3
parallel ∆t min {ν, k min }