Operator splitting for two-dimensional incompressible fluid equations

We analyze splitting algorithms for a class of two-dimensional fluid equations, which includes the incompressible Navier-Stokes equations and the surface quasi-geostrophic equation. Our main result is that the Godunov and Strang splitting methods converge with the expected rates provided the initial data are sufficiently regular.

The quasi-geostrophic equation has recently received considerable attention from the mathematical community. In particular, since the two-dimensional Navier-Stokes equations admit smooth solutions, it has been a question if the quasigeostrophic equation exhibits similar behavior. In this respect, it is common to distinguish between three cases of α for the geostrophic equation. When α ∈ (1, 2), the dissipation term Λ α provides enough regularization to guarantee the existence of smooth solutions [6]. When α ∈ (0, 1), the global properties of solutions are still open. The remaining case (α = 1) is known as the critical case, and the global behavior of solutions was settled only recently (cf. [1,14,8] and the references therein).
In terms of the physical applicability of (1.1), it seems like the most relevant models are the critical geostrophic equation and of course the Navier-Stokes equations. In particular, the critical geostrophic equation has been proposed as a simplified model for strongly rotating atmospheric flow. We refer the reader to [19] for more on the physical aspects of the model.
We now turn to the main topic of the present paper, namely operator splitting algorithms for computing approximate solutions to (1.1). Generally speaking, the label "operator splitting" alludes to the well-known idea of constructing numerical methods for an intricate partial differential equation by reducing the original equation to a succession of simpler equations, each of which can be handled by some efficient and tailor-made numerical method. The operator splitting approach has been comprehensively described in a large number of articles and books. We do not survey the literature here, referring the reader to the bibliography in [9].
Regarding the Navier-Stokes equations, which is a special case of (1.1), operator splitting (viscous splitting) has been analyzed and applied in a great number of works, see, e.g., the book by Majda and Bertozzi [17]. Indeed, viscous splitting has been frequently utilized as a design principle for numerical methods for the Navier-Stokes equations, including vortex or particle methods and transport-diffusion or characteristic-Galerkin methods. Error estimates for Godunov and Strang viscous splitting algorithms have been established, e.g., in [17], utilizing arguments that are different from ours and which rely on the well-known fact that there exist unique, global smooth solutions to the two dimensional Euler and Navier-Stokes equations.
In this paper we apply operator splitting to separate the effects in (1.1) of the transport term u = curl Λ −β θ · ∇θ and the fractional diffusion term Λ α θ. This type of splitting is reasonable as it allows for specialized hyperbolic methods to be applied in the transport step and specialized "Fourier space" methods in the fractional diffusion step. The interested reader can consult [9] for further information on operator splitting.
Our main contribution is that we contribute rigorous proofs of the expected convergence rates for operator splitting applied to the general active scalar equation (1.1). Our approach is inspired by the recent paper [10] (see also [11]) on splitting algorithms for the KdV equation, and for that reason our results apply under the standing assumption that there exists a smooth solution to (1.1). This assumption is verified for the Navier-Stokes equations and the quasi-geostrophic equation with α ≥ 1. It is also reasonable to expect the existence of unique smooth solution in the regime α ∈ [1, 2] and β ∈ [1, 2], cf. [4,12,13,18] for results in that direction.
Let us now discuss our splitting methods in more details. For this purpose, we write (1.1) in the form: We will need the solution operators Φ A (t, θ 0 ) and Φ B (t, θ 0 ), defined as the solutions to the abstract differential equations: For local-in-time existence results for the inviscid equation θ t + u · ∇θ = 0 (i.e., existence of the Φ B operator) when the initial data belong to Sobolev or Triebel-Lizorkin spaces and β = 1, see [5,7,2]; although for β = 1 such results cannot be found in the literature, they can be proved by properly adapting the arguments in [5,7,2]. Regarding the Φ A operator, the fractional diffusion equation has a solution given by v(t) = G α (t) v 0 , where G α (t, x) is the fundamental solution in R 2 which can be expressed in terms of the Fourier transform G α (t, ξ) = e −t|ξ| α .
If G α (x) denotes the inverse Fourier transform of e −|ξ| α , then We refer to [16,20] for further information.
The first operator splitting method we will study is known in the literature as Godunov splitting. The method is defined as follows: Set θ 0 = θ 0 and sequentially determine approximate solutions θ n , n = 1, . . . , M , satisfying Formally, one can show that θ(n∆t)−θ n = O(∆t) in an appropriate spatial norm in the limit ∆t → 0 and n∆t = t. In Section 2.2, we will rigorously prove this linear convergence rate. Specifically, we show that (for small ∆t) The second method we will consider is known as Strang splitting: Let θ 0 = θ 0 and determine sequentially Formally, one can show that the method is second order in ∆t. That is, θ(n∆t) − θ n = O(∆t 2 ) in an appropriate spatial norm in the limit ∆t → 0 and n∆t = t. In Section 2.2, we show that for any sufficiently small time step ∆t > 0 and for all k ≥ max{5, 3α}.
To build fully discrete numerical methods for the fluid equation (1.1), we have to replace the exact solutions operators Φ A and Φ B by appropriate numerical methods. However, we will not discuss that here.
The paper is organized as follows: Section 2 is of an introductory nature and collect some results to be used later on. The convergence rate result for the Godunov splitting is proved in Section 3, while the Strang splitting is analyzed in Section 4.

Preliminary material
2.1. Existence and regularity results. Presently there is no complete existence theory for the equation (1.1). In this paper, we will assume the existence of a unique solution with the same regularity as the initial data. In the literature, one can find results confirming this assumption for some specific cases of α and β.
In what follows we shall also need a local-in-time existence results for Sobolev regular solutions to the inviscid version of (1.1). However, as mentioned in the introduction, the currently available results [5,7,2] apply only to the inviscid quasi-geostrophic equation (β = 1) with initial data θ| t=0 = θ 0 belonging to some Sobolev space H k . Throughout this paper we will simply make the standing assumption that the inviscid generalized quasi-geostrophic equation (θ t + u · ∇θ = 0, u = curl Λ −β θ for β ∈ [1, 2]) possesses such a sufficiently regular solution.
2.2. Fractional calculus. The fractional Laplace operator Λ α occurring in (1.1) is defined using Fourier transform, namely Our normalization of the Fourier transform reads In the upcoming analysis, we will need a different representation of Λ α . Since we require α ∈ (0, 2), [7] provides the identity , for all f in the Schwartz class (in particular all f ∈ C ∞ 0 ). Here, P.V. · · · dz denotes the principal value integral, and we have introduced the notation We will make use of the following Leibniz-like formula "with remainder": By adding and subtracting, we see that Multiplying (2.1) with |z| −2−β and integrating over z provides the following representation of G β .

Lemma 2.3 (Leibniz formula
). If f and g are sufficiently smooth functions the following identity holds pointwise, for α = 2.
In our analysis, we will make several applications of the above Leibniz rule. The following proposition provides an L p estimate for the term G α . It is a variation of a result due to Constantin [3].
Proposition 2.4. Let 1 < p < ∞. For each fixed α ∈ (0, 2], there is a constant C depending on α such that Proof. The result follows directly from the Hölder inequality when α = 2. We may thus assume that α ∈ (0, 2). Let us write G α (f, g) as the sum of two parts: We commence by estimating the first term: 1 qi = 1, we can apply the generalized Hölder inequality to obtain In order to work efficiently in our Sobolev spaces we will need to provide some standard definitions, mostly to fix the notation. Let l denote a two-dimensional multi-index, i.e., l = (l 1 , l 2 ), l j ∈ N 0 . Then we write We will be working the Sobolev spaces (where S denotes the set of tempered distributions). If k is a natural number, H k is the standard Sobolev space with inner product and norm given by In our analysis, we will apply the following corollary of the previous proposition.
Our strategy is to prove the desired estimate for each of terms in the sum separately. For this purpose, we let s = 0, . . . , k be arbitrary and consider an arbitrary component of ∇ s G α (f, g). Using Lemma 2.3 and the standard Leibniz rule, we deduce ∂ s ∂x l ∂y s−l G α (f, g) For m, n such that 2 ≤ m + n ≤ k, Proposition 2.4 can applied to conclude Conversely, if m,n is such that 0 ≤ m + n ≤ 1, Proposition 2.4 allow us deduce since k ≥ 3. Now, by taking the L 2 norm on both sides of (2.2) and applying the previous calculations, we gather Hence, any given component of ∇ s G α (f, g) satisfies the desired bound. Thus, which completes our proof.
2.3. Two auxiliary lemmas. We will make heavily use of the following two lemmas throughout the paper. Their proofs are technical and tedious, but straightforward. For this reason, proofs are deferred to the appendix.
Lemma 2.7. Let k ≥ 4 be an integer. The following estimates hold (2.5)

Godunov splitting
In this section we prove rigorously the expected linear rate of convergence for Godunov splitting. As part of the proof we also show that this splitting method is well-defined and that it produces regular approximations. The results are valid under a condition on the length of the time step ∆t.
We begin by precisely defining Godunov splitting for (1.1). For this purpose, write (1.1) as Using the operators A and B, we define the solution operators Φ A and Φ B as the solutions to The Godunov splitting method is classically defined as follows: For ∆t > 0 given, construct a sequence {θ n , θ n+1/2 } T /∆t n=0 of approximate solutions to (1.1) by the following procedure: Let θ 0 = θ 0 and determine inductively To facilitate the convergence analysis, we will need a different definition of the Godunov method. Our definition can be seen as an extension of the splitting solution {θ n , θ n+1/2 } n to all of [0, T ]. The most used method of extension is to let "time run twice as fast" in each of the sub-intervals [t n , t n+1/2 ] and [t n+1/2 , t n+1 ], where as usual t r = r∆t for r ∈ [0, ∞), thus obtaining Although it appears to be a natural extension, it does not seem to be appropriate for our purpose. Instead, we will follow the approach taken in the paper [10], and introduce two time variables instead of one.
Our Godunov splitting method is given by the following definition.
The time-continuus Godunov splitting solution ϑ : Ω ∆t → R is defined as the solution to cf. Figure 1.
To measure the error, we will use the function where θ is the (smooth) solution of (1.1). It is not trivial trivial to obtain the existence of a splitting solution ϑ in the sense of Definition 3.1. Since the best available existence result we have for the hyperbolic step is local-in-time, it is unclear if we can iterate the steps up to any given large time T . In fact, well-posedness of the method is one of our main results.
The following theorem is our main result in this section. 2]. Then, for ∆t sufficiently small, we have the following: Theorem 3.2 will be an outcome of the results proved in the subsections below (Subsection 3.3 will bring the pieces together).
3.1. Evolution equations for the error. To prove Theorem 3.2, we will analyze a set of evolution equations satisfied by the error e, see [10], which we now derive.
We shall need the following Taylor expansion satisfied by an operator E: Using the definition of ϑ and the above Taylor formula, we deduce where we have introduced the "forcing" term By direct calculation, where we have defined the commutator For the fluid equation (1.1), we have that and The Leibniz formula (Lemma 2.3) yields Hence, for the fluid equation (1.1), equations (3.2)-(3.3) read: The equations (3.4) and (3.5) constitute our main tool for proving Theorem 3.2.

3.2.
Estimates on the error. Lemma 3.5 below provides an estimate on the error e that will be a key ingredient in the proof of Theorem 3.2. The estimate depends on two auxiliary results (Lemmas 3.3 and 3.4) that we first prove. The first auxiliary results gives perhaps the most fundamental property of our splitting solution. It states that if the splitting solution is in H k , k ≥ 6, then it is actually in H k+2 .
The following lemma is our second auxiliary result.
Lemma 3.4. Let k, t, and τ be as in the previous lemma. Then, (3.4), multiplying componentwise with ∇ s F , summing over s = 0, . . . , k − 2, and integrating over the domain gives (3.9) An application of Corollary 2.5 to (3.9) gives (here we do not make use of the regularizing effect of Λ −β ) In view of Lemma 3.3, this means that which concludes the proof.
The next lemma will be the key ingredient in the proof of Theorem 3.2.
Let H(t, τ ) denote the statement ϑ(s, σ) H k−2 ≤ γ, (s, σ) ∈ Ω t,τ ∆t , and C(t, τ ) the statement ∈ Ω t,τ ∆t , where γ is some value which will specified below. Let us for a moment assume that assertions (1)-(4) of Lemma 3.6 are true for H and C. In this case Lemma 3.6 tells us that C(t, τ ) holds for all (t, τ ) ∈ Ω ∆t . Hence, the H k−2 norm of the splitting solution never blows up and thus the existence part of Theorem 3.2 follows readily (a local-in-time solution can be extended to all of (0, T )).
Let us now verify assertions (1)-(4) of Lemma 3.6 for our choice of H(t, τ ) and C(t, τ ). First, we observe that assertions (2) and (3) clearly hold. For assertion (4) to be true, we assume that γ satisfies (3.13) Then it remains to verify assertion (1). For this purpose, let us assume that H(t, τ ) is true for some (t, τ ) ∈ Ω ∆t . Lemma 3.5 can then be applied to obtain the bound (3.14) We need to compare ϑ(s, σ) with the corresponding value on the diagonal, viz. ϑ(s, s). Two cases need to be considered: (i) If σ > s, we observe that (using that α ≤ 2 to make sure the expression is finite) which implies that Using this inequality and applying Lemma 3.3 we find Thus in both cases we conclude that ϑ(s, σ) H k−2 ≤ ϑ(s, s) H k−2 + ∆t e Cγs θ 0 H k , |s − σ| ≤ ∆t.
By applying the previous inequality, adding and subtracting θ, and involving (3.14), we estimate (3.15) Now, we fix γ and ∆t according to and note that this is not conflict with (3.13). Then, (3.15) gives ∈ Ω t,τ ∆t , which verifies (1) in Lemma 3.6.
At this point we have verified assertions (1)-(4) of Lemma 3.6 for our choice of H(t, τ ) and C(t, τ ). Consequently, Lemma 3.6 tells us that C(t, τ ) is true for all times (t, τ ) ∈ Ω ∆t . In other words,

Strang splitting
In this section we prove the expected second-order convergence rate of Strang splitting applied to (1.1). We also prove that the method produces regular solutions up to any given finite final time T > 0. The results are valid under a condition on the size of each time step ∆t.
We will continue to use the notation and definitions introduced in the previous section, unless explicitly stating otherwise. The Strang method we consider in this section reads as follows: For ∆t > 0 given, construct a sequence {θ n } T /∆t n=0 of approximate solutions to (1.1) by the following procedure: Let θ 0 = θ 0 and determine sequentially As with the Godunov splitting, the convergence analysis will require a timecontinuous interpolation of the Strang approximations {θ n } n . In contrast to the Godunov case, we will now introduce three time variables instead of two. This approach is different from the one taken in [10], and indeed appears more natural. In [10] the authors stick to two time variables and interpret Strang splitting in terms of two "∆t/2" Godunov splittings, and in alternating order.
The time-continuous Strang splitting solution ϑ : Ω ∆t → R is defined as cf. Figure 2.
In each box tn , it is of particular interest to consider the function ϑ restricted to the diagonal, i.e., the function ϑ t 2 , t, t 2 for t ∈ [t n , t n+1 ]. Observe that at each point on the diagonal, ϑ is a Strang splitting solution with a specific time step. Specifically, ϑ t 2 , t, t 2 is a Strang splitting solution with time step t − t n . Consequently, ϑ t n 2 , t n , t n 2 = θ n , n = 0, . . . , T /∆t , and hence that ϑ t 2 , t, t 2 can be seen as an extension of {θ n } n to all of [0, T ]. To measure the error between the splitting approximation and the exact solution, we will use the function where θ is the (smooth) solution of (1.1). We state two main results in this section.
for some constant C = C ( θ 0 H k ). 4.1. The error evolution equations. To prove Lemma 4.2 we will use the same approach as for the Godunov method. However, since we are now using three time variables instead of two, we must derive a new set of evolution equations governing the error.
By a direct calculation, the error e satisfies the time-evolution where F (t) = F t 2 , t, t 2 and , τ, ω)).
Let Ω t,τ,ω ∆t denote the set of times prior to (t, τ, ω): (4.5) We begin by observing that Lemma 3.3 can be easily adapted to yield the following result.
Lemma 4.4. Let 6 ≤ k ∈ N, and assume the existence of (t, τ, ω) ∈ Ω ∆t such that Then Next we prove a Strang version of Lemma 3.4.
Note that while F = 0 on the line (s, t n , tn 2 ), there holds 1 2 (ϑ t − B(ϑ)) = 0 on (s, t n , tn 2 ). Hence, we can repeat the arguments of Lemma 3.4 to conclude that ∈ Ω t,τ,ω ∆t . (4.8) We fix n such that (t, τ, ω) . Now, to estimate F H k−2 at an arbitrary point (s, σ, ζ) ∈ Ω t,τ,ω ∆t , we will first integrate in the ω direction to the plane given by ω = tn 2 and then apply the estimate (4.8). To perform this integration, we apply ∇ s to (4.3), multiply the result with ∇ s F , sum over s = 0, . . . , k − 2, and integrate by parts, to obtain 9) and the goal is to bound the terms I 1 , I 2 , and I 3 .
By first estimating as in (3.9), then using Corollary 2.5, and finally applying Lemma 4.4, we obtain To bound the two other terms I 2 and I 3 , we first apply Lemma 2.7 to find

An application of Lemma 4.4 then yields
Inserting these bounds back into (4.9), Applying the Gronwall lemma we obtain Finally, we apply (4.8) to conclude which completes the proof. Recall that Ω t,τ,ω ∆t denotes the set of all times prior to (t, τ, ω), cf. (4.5). To prove the existence part of Lemma 4.2 we will utilize the same strategy as with did for the Godunov method in the previous section. That is, to apply the bootstrap lemma (Lemma 3.6) with H(t, τ, ω) as the statement ϑ(s, σ, ζ) H k−2 ≤ γ, (s, σ, ζ) ∈ Ω t,τ,ω ∆t , (4.10) and C(t, τ, ω) the statement where γ is some value to be specified below. Note that there is no problem with extending Lemma 3.6 so that it allows three time variables instead of two (i.e., hypothesis and statement depending on three variables).
By definition, where the estimate follows from Lemma 2.6, and where we have used (4.10). The fundamental theorem of calculus and (4.13) provides us with the estimate and note that this is not conflict with (4.11). Then, (4.16) gives Since (s, σ, ζ) ∈ Ω t,τ,ω ∆t was arbitrarily chosen, this verifies (1) in Lemma 3.6. We have now verified assertions (1)-(4) of Lemma 3.6 for H(t, τ, ω) and C(t, τ, ω). Consequently, Lemma 3.6 tells us that C(t, τ, ω) is true for all times (t, τ, ω) ∈ Ω ∆t . In other words, By applying the definition of Θ q+1 i,j, , the definition ϑ ω (cf. (4.1)), and the Leibniz rule, we deduce : ∇ s Θ q+1 i,j, (t , τ , s) dxds. Inserting this expression into the previous equality yields : ∇ s Θ q+1 i,j, (t , τ , s) dxds. Let us consider three separate cases of m + n + r in the quadruple sum above. (i) If m + n + r = q + 1, the corresponding term in the above reads Here, the first inequality is an application of Lemma 2.7 and the last inequality is an application of Lemma 4.4.
(ii) If m + n + r = 0, we can also apply Lemma 2.7 to conclude (iii) For the remaining cases (1 ≤ m + n + r ≤ q), we first apply the Leibniz rule to the spatial derivatives, then we overestimate the terms using the Hölder inequality.
The following lemma is almost a corollary of the previous lemma.  1). Then, for any k such that θ 0 H k+3α ≤ C, we have ∈ Ω ∆t . Proof. Let i and j denote any one of t, τ , or ω. An arbitrary component of ∇ 2 t F can then be written F ij := ∂ i ∂ j F . By definition, we have that By applying the Hölder inequality together with the previous lemma, we estimate where we have grossly overestimated most of the terms. We have also applied Lemma 4.7 with l = 1 for the H k+2 norm and with l = 2 for the H k norm. The constant C depends on θ 0 H k+3α . Theorem 4.3 is a consequence of the following (remarkable) fact: Proof. We will prove Lemma 4.9 by direct calculation. Let us begin by estimating F ω tn 2 , t n , tn 2 . Since F (t, t n , tn 2 ) = 0, (4.3) tells us that Similarly, we see that (4.2) yields It remains to estimate F t tn 2 , t n , tn 2 . However, as F (t, t n , tn 2 ) = 0, we must have F t tn 2 , t n , tn 2 = 0. This, together with (4.26) and (4.27), concludes the proof.
Using the previous lemma, we can now prove that the error produced along the diagonal (t/2, t, t/2) is second order in ∆t.  where the last equality is an application of Lemma 4.9. By taking the H k norm on both sides of (4.29) and applying the previous lemma, we gather which concludes the proof.

4.4.
Proof of Theorem 4.3. We have now gathered all the ingredients needed to prove the sought after second-order error estimate.
Performing the same calculations as in (3.10), (3.11), (3.12), and then applying Lemma 4.10, yields Since e(0) = 0, an application of the Gronwall inequality to the previous inequality gives which concludes the proof of Theorem 4.3.
Appendix A. Proof of Lemmas 2.6 and 2.7 The purpose of this appendix is to provide proofs of Lemmas 2.6 and 2.7. Both lemmas have been crucial to our convergence analysis. In particular, Lemmas 3.3, 3.5, 4.5, and 4.7, all rely on their validity. Proof. Let us confine to the three-dimensional case (N = 3) as the other cases are almost identical. By applying the Leibniz rule (with multiindex notation α = (α 1 , α 2 , α 3 )), we obtain the following expression Let us now consider four separate cases of i 1 + i 2 + i 3 in the above quadruple sum.
Lemma A.2. Let k ≥ 4 be an integer. The following estimates hold Proof. The proof of (2.4) is easily obtained by the calculations of the previous proof. To prove (2.5), it is only step (i) of the previous proof which is no longer true. However, this is also the reason for the k + 1 on g. That is, step (i) is replaced by a simpler Hölder inequality.