# An efficient algorithm for simulating ensembles of parameterized flow problems

An efficient algorithm for simulating ensembles of parameterized flow problems Abstract Many applications of computational fluid dynamics require multiple simulations of a flow under different input conditions. In this paper, a numerical algorithm is developed to efficiently determine a set of such simulations in which the individually independent members of the set are subject to different viscosity coefficients, initial conditions and/or body forces. The proposed scheme, when applied to the flow ensemble, needs to solve a single linear system with multiple right-hand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and long-term stable under certain conditions on the time-step size and a parameter deviation ratio. A rigorous numerical error estimate shows the scheme is of first-order accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results. 1. Introduction Numerical simulations of incompressible viscous flows have important applications in engineering and science. In this paper we consider settings in which one wishes to obtain solutions for several different values of the physical parameters and several different choices for the forcing functions appearing in the partial differential equation (PDE) model. For example, in building low-dimensional surrogates for the PDE solution such as sparse-grid interpolants or proper orthogonal decomposition approximations, one has to first determine expensive approximation of solutions corresponding to several values of the parameters. Sensitivity analyses of solutions often need to determine approximate solutions for several parameter values and/or forcing functions. An important third example is quantifying the uncertainties of outputs from the model equations. Mathematical models should take into account the uncertainties invariably present in the specification of physical parameters and/or forcing functions appearing in the model equations. For flow problems, because the viscosity of the liquid or gas often depends on the temperature, an inaccurate measurement of the temperature would introduce some uncertainty into the viscosity of the flow. Direct measurements of the viscosity using flow meters and measurements of the state of the system are also prone to uncertainties. Of course, forcing functions, e.g., initial condition data, can and usually are also subject to uncertainty. In such cases, due to the lack of exact information, stochastic modeling is used to describe flows subject to a random viscosity coefficient and/or random forcing. Subsequently, numerical methods are employed to quantify the uncertainties in system output. It is known that uncertainty quantification, when a random sampling method such as a Monte Carlo method is used, could be computationally expensive for large-scale problems because each individual realization requires a large-scale computation but on the other hand, many realizations may be needed in order to obtain accurate statistical information about the outputs of interest. Therefore, for all the examples discussed and for many others, how to design efficient algorithms for performing multiple numerical simulations becomes a matter of great interest. The ensemble method that forms the basis for our approach was proposed in Jiang & Layton (2014); there, a set of J solutions of the Navier–Stokes equations (NSE) with distinct initial conditions and forcing terms is considered. All solutions are found, at each time step, by solving a linear system with one shared coefficient matrix and J right-hand sides (RHS), reducing both the storage requirements and computational costs of the solution process. The algorithm of Jiang & Layton (2014) is first-order accurate in time; it is extended to higher-order accurate schemes in Jiang (2015, 2017). Ensemble regularization methods are developed in the studies by Jiang (2015), Jiang & Layton (2015) and Takhirov et al. (2016) for high Reynolds number flows, and a turbulence model based on ensemble averaging is developed in Jiang et al. (2015). The ensemble algorithm has also been extended to magnetohydrodynamics flows in Mohebujjaman & Rebholz (2017), to natural convection problems in Fiordilino (2017a) and to parametrized flow problems in Gunzburger et al. (2017b). Ensemble algorithms incorporating reduced-order modeling techniques are studied in Gunzburger et al. (2017a, 2018). Recently, the ensemble method has been introduced in Fiordilino (2017b); Luo & Wang (2018a,b) for uncertainty quantification problems on random linear parabolic equations. In this paper, we develop a numerical scheme for simulating ensembles of the NSE flow problems in which not only the initial data and body force function, but also the viscosity coefficient, may vary from one ensemble member to another. Specifically, we consider a set of J NSE simulations on a bounded domain subject to no-slip boundary conditions in which the jth individual member solves the system   $$\left\{ \begin{array}{rcll} u_{j, t}+u_{j}\cdot\nabla u_{j}-\nabla \cdot (\nu_j \nabla u_{j})+\nabla p_{j} &=& f_{j}(x,t) \quad &\textrm{in}\ \varOmega\times [0, \infty), \\ \nabla\cdot u_{j} &=& 0 \quad &\textrm{in}\ \varOmega \times [0, \infty), \\ u_{j} &=& 0 \quad &\textrm{on}\ \partial\varOmega \times [0, \infty),\\ u_{j}(x,0) &=& u_{j}^{0}(x) \quad & \textrm{in} \ \varOmega, \end{array}\right.$$ (1.1)which, for each j, corresponds to a different variable kinematic viscosity $$\nu _j= \nu _j(x)$$ and/or distinct initial data $$u_{j}^{0}$$ and/or body forces $$f_{j}$$. In the sequel, it is assumed that $$\nu _j(x) \in L^\infty (\varOmega )$$ and $$\nu _j(x)\geqslant \nu _{j, \textrm{min}}>0$$. Due to the nonlinear convection term, implicit and semiimplicit schemes are invariably used for time integration. For a semiimplicit scheme the associated discrete linear systems would be different for each individually independent simulation, i.e., for each j. As a result, at each time step, J linear systems need to be solved to determine the ensemble, resulting in a huge computational effort. For a fully implicit scheme, the situation is even worse because one would have to solve many more linear systems due to the nonlinear solver iteration. To tackle this issue we propose a novel discretization scheme that results, at each time step, in a common coefficient matrix for all the ensemble members. 1.1 The ensemble-based semiimplicit scheme To focus on the main idea we temporarily ignore the spatial discretization and consider only the ensemble-based implicit–explicit temporal integration scheme   \left\{ \begin{aligned} \frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}+\overline{u}^{n}\cdot\nabla u_{j}^{n+1} +(u_{j}^{n}-\overline{u}^{n})\cdot\nabla u_{j}^{n} -\nabla \cdot (\overline{\nu} \nabla u_{j}^{n+1}) -\nabla \cdot \left((\nu_j-\overline{\nu})\nabla u_{j}^{n}\right) +\nabla p_{j}^{n+1} &=f_{j}^{n+1},\\ \nabla\cdot u_{j}^{n+1}&=0, \end{aligned}\right. (1.2)where $$\overline{u}^{n}$$ and $$\overline{\nu }$$ are the ensemble means of the velocity and viscosity coefficients, respectively, defined as   $$\overline{u}^{n}:=\frac{1}{J}\sum_{j=1}^{J}u_{j}^{n} \qquad \textrm{and}\qquad \overline{\nu}:=\frac{1}{J}\sum_{j=1}^{J}\nu_{j}.$$We also define $$\overline{\nu }_{\textrm{min}}:= \frac{1}{J}\sum _{j=1}^{J}\nu _{j, \textrm{min}}$$. After rearranging the system, we have, at time $$t_{n+1}$$,   \left\{ \begin{aligned} \frac{1}{\varDelta t}u_{j}^{n+1}\!+\overline{u}^{n}\cdot\nabla u_{j}^{n+1} \!-\! \nabla\!\cdot\! (\overline{\nu} \nabla u_{j}^{n+1}) \!+\!\nabla p_{j}^{n+1} &= f_{j}^{n+1}\!+\! \frac{1}{\varDelta t}u_{j}^{n} \!-(u_{j}^{n}\!-\overline{u}^{n})\cdot\nabla u_{j}^{n} \!+\!\nabla \!\cdot\! \left( (\nu_j\!-\!\overline{\nu}) \nabla u_{j}^{n} \right), \\ \nabla\cdot u_{j}^{n+1}&=0. \end{aligned}\right. (1.3)It is clear that the coefficient matrix of the resulting linear system will be independent of j. Thus, for the flow ensemble, to advance all members of the ensemble by one time step, we need to solve only a single linear system with J RHS. Compared with solving J individually independent simulations, this approach used with either a single LU factorization for small-scale problems or a block Krylov subspace method (Gutknecht, 2007; Parks et al., 2016) for large-scale problems is computationally more efficient and significantly reduces the required storage. When the size of the ensemble becomes huge, it can be subdivided into p sub-ensembles so as to balance memory, communication and computational costs and then (1.2) can be applied to each sub-ensemble. The rest of this section is devoted to establishing notation and to providing other preliminary information. Then in Section 2 we prove a conditional stability result for a fully discrete finite element discretization of (1.2). In Section 3 we derive an error estimate for the fully discrete approximation. Results of the preliminary numerical simulations that illustrate the theoretical results are given in Section 4, and Section 5 provides some concluding remarks. 1.2 Notation and preliminaries Let $$\varOmega$$ denote an open, regular domain in $$\mathbb{R}^{d}$$ for d = 2 or 3 having boundary denoted by $$\partial \varOmega$$. The $$L^{2}(\varOmega )$$ norm and inner product are denoted by ∥ ⋅ ∥ and (⋅, ⋅), respectively. The $$L^{p}(\varOmega )$$ norms and the Sobolev $$W^{k}_{p}(\varOmega )$$ norms are denoted by $$\|\cdot \|_{L^{p}}$$ and $$\|\cdot \|_{W_{p}^{k}}$$, respectively. The Sobolev space $$W_{2}^{k}(\varOmega )$$ is simply denoted by $$H^{k}(\varOmega )$$ and its norm by $$\|\cdot \|_{k}$$. For functions v(x, t) defined on (0, T), we define, for $$1\leqslant m <\infty$$,   $$\| v \|_{\infty,k} \textrm{ }:=\textrm{EssSup}_{[0,T]}\| v(\cdot, t)\|_{k}\qquad \textrm{and}\qquad \|v\|_{m,k} \textrm{ }:= \left( \int_{0}^{T}\|v(\cdot, t)\|_{k}^{m}\, \textrm{d}t \right) ^{1/m}.$$Given a time step $$\varDelta t$$, associated discrete norms are defined as   $${\left\vert \! \left\vert \! \left\vert v \right\vert\! \right\vert \! \right\vert}_{\infty,k}=\max\limits_{0\leqslant n\leqslant N}\Vert v^{n}\Vert_{k} \qquad \textrm{and}\qquad{\left\vert \! \left\vert \! \left\vert v \right\vert\! \right\vert \! \right\vert}_{m,k}:= \left(\sum_{n=0}^{N}\|v^{n}\|_{k}^{m}\varDelta t \right)^{1/m},$$where $$v^{n}=v(t_{n})$$ and $$t_{n}=n\varDelta t$$. Denote by $$H^{-1}(\varOmega )$$ the dual space of bounded linear functions on $$H^{1}_{0}(\varOmega )=\{v\in H^{1}\,:\, v=0 \ \textrm{on}\ \partial \varOmega \}$$; a norm on $$H^{-1}(\varOmega )$$ is given by   $$\|\,f\|_{-1}=\sup_{0\neq v\in H^{1}_{0}(\varOmega)}\frac{(\,f,v)}{\Vert\nabla v\Vert}.$$ The velocity space X and pressure space Q are given by   $$X:=[H_{0}^{1}(\varOmega)]^{d} \qquad \textrm{and}\qquad Q:=L_{0}^{2}(\varOmega)=\{q\in L^2(\varOmega)\,\,:\,\, \textstyle\int_\varOmega q\, \textrm{d}x =0\},$$respectively. The space of weakly divergence-free functions is   $$V\textrm{ }:=\{v\in X\,\,:\,\, (\nabla\cdot v,q)=0,\quad \forall\, q\in Q\} .$$ A weak formulation of (1.1) reads as follows: for j = 1, … , J, find $$u_j: \,[0,T]\rightarrow X$$ and $$p_j:\, [0,T]\rightarrow Q$$ for a.e. t ∈ (0, T] satisfying   \begin{equation*} \left\{ \begin{aligned} (u_{j,t},v)+(u_{j}\cdot\nabla u_{j},v)+(\nu_j\nabla u_{j},\nabla v) -(p_{j},\nabla\cdot v) &= (f_{j},v) &\forall\, v\in X, \\ (\nabla\cdot u_{j},q) &= 0 &\forall\, q\in Q, \end{aligned}\right. \end{equation*}with $$u_{j}(x,0)=u_{j}^{0}(x)$$. Our analysis is based on a finite element method (FEM) for spatial discretization. However, the results also extend, without much difficulty, to other variational discretization methods. Let $$X_{h}\subset X$$ and $$Q_{h}\subset Q$$ denote families of conforming velocity and pressure finite element spaces on a regular subdivision of $$\varOmega$$ into simplices; the family is parameterized by the maximum diameter h of any of the simplices. Assume that the pair of spaces $$(X_h,Q_h)$$ satisfy the discrete inf-sup (or $$\textrm{LBB}_h$$) condition required for the stability of the finite element approximation and that the finite element spaces satisfy the approximation properties   $$\inf_{v_h\in X_h}\| v- v_h \| \leqslant C h^{k+1}\Vert u \Vert_{k+1} \quad \quad\quad\quad \forall\, v\in [H^{k+1}(\varOmega)]^d,$$ (1.4)  $$\inf_{v_h\in X_h}\| \nabla ( v- v_h )\| \leqslant C h^k \Vert v\Vert_{k+1}\quad \;\;\,\quad\quad\quad \forall\, v\in [H^{k+1}(\varOmega)]^d,$$ (1.5)  $$\inf_{q_h\in Q_h}\| q- q_h \| \leqslant C h^{s+1}\Vert p\Vert_{s+1} \quad \quad\quad\quad\quad \forall\, q\in H^{s+1}(\varOmega),$$ (1.6)where the generic constant C > 0 is independent of mesh size h. An example for which the $$\textrm{LBB}_h$$ stability condition and the approximation properties are satisfied is the family of Taylor–Hood $$P^{s+1}$$–$$P^{s}$$, $$s\geqslant 1$$ element pairs. For details concerning FEMs see, e.g., Ciarlet (2002) and Girault & Raviart (1979, 1986); Gunzburger (1989); Layton (2008) for FEMs for the NSE. The discretely divergence-free subspace of $$X_{h}$$ is defined as   $$V_{h}\ :=\{v_{h}\in X_{h}\,\,:\,\,(\nabla\cdot v_{h},q_{h})=0\quad \forall\, q_{h}\in Q_{h}\}.$$Note that, in general, $$V_h\not \subset V$$. We assume the mesh and finite element spaces satisfy the standard inverse inequality   $$h\Vert\nabla v_{h}\Vert \leqslant C_{(\textrm{inv})}\Vert v_{h}\Vert \quad\forall\, v_{h}\in X_{h},$$ (1.7)which is known to hold for standard finite element spaces with locally quasi-uniform meshes (Brenner & Scott, 2008). We also define the standard explicitly skew-symmetric trilinear form   $$b^{\ast}(u,v,w):=\tfrac{1}{2}(u\cdot\nabla v,w)-\tfrac{1}{2}(u\cdot\nabla w,v),$$which satisfies the bounds (Layton, 2008)   $$b^{\ast}(u,v,w)\leqslant C \left(\Vert \nabla u\Vert\Vert u\Vert\right)^{1/2}\Vert\nabla v\Vert\Vert\nabla w \Vert \quad \forall\, u, v, w \in X,$$ (1.8)  $$b^{\ast}(u,v,w)\leqslant C \Vert \nabla u\Vert\Vert\nabla v\Vert\left(\Vert\nabla w \Vert\Vert w\Vert\right)^{1/2} \quad \forall\, u, v, w \in X .$$ (1.9)We also denote the exact and approximate solutions at $$t=t^n$$ as $$u_{j}^{n}$$ and $$u_{j, h}^{n}$$, respectively. 2. Stability analysis The fully discrete finite element discretization of (1.2) is given as follows. Given $$u_{j,h}^{0}\in X_{h}$$, for n = 0, 1, … , N − 1, find $$u_{j,h}^{n+1}\in X_{h}$$ and $$p_{j,h}^{n+1}\in Q_{h}$$ satisfying   $$\begin{cases} \left(\frac{u_{j,h}^{n+1}-u_{j,h}^{n}}{\varDelta t},v_{h}\right)+b^{\ast}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h})+b^{\ast}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n}, v_{h})-(p_{j,h}^{n+1},\nabla\cdot v_{h})\\ \qquad\qquad\quad\;\ \ \ \ + \, (\overline{\nu} \nabla u_{j,h}^{n+1},\nabla v_{h}) +\left( (\nu_j-\overline{\nu}) \nabla u_{j,h}^{n},\nabla v_{h} \right) =(f_{j}^{n+1},v_{h})\quad \forall\, v_{h}\in X_{h},\\ \big(\nabla\cdot u_{j,h}^{n+1},q_{h}\big)=0 \quad\forall\, q_{h}\in Q_{h}. \end{cases}$$ (2.1)We begin by proving the conditional, nonlinear, long-time stability of scheme (2.1) under a time-step condition and a parameter deviation condition. Theorem 2.1 (Stability). For all j = 1, … , J, if for some $$\mu$$, $$0\leqslant \mu <1$$, and some $$\varepsilon$$, $$0< \varepsilon \leqslant 2-2\sqrt{\mu }$$, the following time-step condition and parameter deviation condition both hold,   $$C\frac{\varDelta t}{ h \overline{\nu}_{\textrm{min}}}\left\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\right\Vert{}^{2} \leqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)},$$ (2.2)  $$\frac{\|\nu_j -\overline{\nu}\|_{\infty} }{\overline{\nu}_{\textrm{min}}} \leqslant \sqrt{\mu},\quad$$ (2.3)then scheme (2.1) is nonlinearly, long-time stable. In particular, for j = 1, … , J and for any $$N\geqslant 1$$, we have   $$\begin{split} &\quad\frac{1}{2}\|u_{j,h}^{N}\|^{2}+\frac{1}{4}\sum_{n=0}^{N-1}\|u_{j,h} ^{n+1}-u_{j,h}^{n}\|^{2}+\sum_{n=0}^{N-1}\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\\ &\qquad+\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_\textrm{{min}}}\right) \|\nabla u_{j,h}^{N}\|^{2}\\ &\qquad\qquad\leqslant\sum_{n=0}^{N-1}\frac{2\varDelta t}{\overline{\nu}_{\textrm{min}}}\|\,f_{j}^{n+1}\|_{-1}^{2}+ \frac{1} {2}\|u_{j,h}^{0}\|^{2} +\overline{\nu}_\textrm{{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\|\nabla u_{j,h}^{0}\|^{2}. \end{split}$$ (2.4) Proof. The proof is given in Appendix A. Remark 2.2 It is seen from (2.2) that the upper bound in the time-step condition increases as $$\varepsilon$$ decreases. As $$\varepsilon \rightarrow 0$$ the bound approaches $$1-\sqrt{\mu }$$. Because the relative deviation of the viscosity coefficient in (2.3) is bounded by $$\sqrt{\mu }$$, the two stability conditions are oppositional to each other. Remark 2.3 Noting that condition (2.2) depends only on known quantities such as the solution at $$t_n$$ and that scheme (2.1) is a one-step method, (2.2) can be used to adapt $$\varDelta t$$ in order to guarantee stability for the ensemble simulations. Remark 2.4 When the viscosity $$\nu _j$$ is a constant instead of being variable, the relative viscosity coefficient deviation ratio required for stability is still bounded by $$\sqrt{\mu }$$, and condition (2.3) becomes $$\max_{j}\frac{\left |\nu _j-\overline{\nu }\right |}{\overline{\nu }}\leqslant \sqrt{\mu }$$. 3. Error analysis In this section we give a detailed error analysis of the proposed method under the same type of time-step condition (with possibly a different constant C on the left-hand side of the inequality) and the same parameter deviation condition. Assuming that $$X_{h}$$ and $$Q_{h}$$ satisfy the $$\textrm{LBB}_{h}$$ condition, scheme (2.1) is equivalent to the following: given $$u_{j,h}^{0}\in V_{h}$$, for n = 0, 1, … , N − 1 find $$u_{j,h} ^{n+1}\in V_{h}$$ such that   \begin{align} \left(\frac{u_{j,h}^{n+1}-u_{j,h}^{n}}{\varDelta t},v_{h}\right)&+b^{\ast}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h})+b^{\ast}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h} ^{n},v_{h})\nonumber\\ &+\overline{\nu}(\nabla u_{j,h}^{n+1},\nabla v_{h})+((\nu_j-\overline{\nu})\nabla u_{j,h}^{n},\nabla v_{h})=(\,f_{j}^{n+1},v_{h})\quad\forall\, v_{h}\in V_{h}. \end{align} (3.1)To analyse the rate of convergence of the approximation we assume the following regularity for the exact solutions:   \begin{align*} u_{j} \in L^{\infty}(0,T;H^{k+1}(\varOmega))\cap H^{1}(0,T;H^{k+1}(\varOmega))\cap H^{2}(0,T;L^{2}(\varOmega)),\\ p_{j} \in L^{2}(0,T;H^{s+1}(\varOmega))\quad \textrm{and}\quad f_{j} \in L^{2} (0,T;L^{2}(\varOmega)). \end{align*}Let $$e_{j}^{n}=u_{j}^{n}-u_{j,h}^{n}$$ denote the approximation error of the jth simulation at the time instance $$t_n$$. We then have the following error estimates. Theorem 3.1 (Convergence of scheme (2.1)). For all j = 1, … , J, if for some $$\mu$$, $$0\leqslant \mu <1$$ and some $$\varepsilon$$, $$0< \varepsilon \leqslant 2-2\sqrt{\mu }$$ the following time-step condition and parameter deviation condition both hold,   $$C\frac{\varDelta t}{\overline{\nu}_{\textrm{min}} h}\left\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\right\Vert{}^{2} \leqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)},$$ (3.2)  $$\frac{\|\nu_j -\overline{\nu}\|_{\infty} }{\overline{\nu}_{\textrm{min}}} \leqslant \sqrt{\mu},\quad$$ (3.3)then there exists a positive constant C independent of the time step such that   \begin{align*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} &+\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left( 1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \\ \leqslant&\, e^{\frac{CT}{\overline{\nu}_{\textrm{min}}^{3}}} \left\{\frac{1}{2}\Vert e_{j}^{0}\Vert^{2} +\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla e_{j}^{0}\Vert^{2}\right. \\ &\qquad\qquad +\frac{1}{2}h^{2k+2}\Vert u_j^0 \Vert_{ k+1}^2+\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\overline{\nu}_{\textrm{min}}\varDelta t h^{2k} \Vert u_{j}^{0}\Vert_{k+1}^{2}\\ &\qquad\qquad+C\varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla u_{j,t}\Vert^2_{2,0}+C\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert \! \left\vert u_j \right\vert\!\right\vert \! \right\vert}^2_{2,k+1} +C\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \\ &\qquad\qquad+Ch^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t \Vert \nabla u_{j,t}\Vert^2_{2,0} + C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \\ &\qquad\qquad+C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2\Vert \nabla u_{j,t}\Vert^2_{2,0}+ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \\ &\qquad\qquad+C\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2}+C\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}\Vert u_{j,t}\Vert_{2,k+1}^{2} +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}\Vert u_{j,tt}\Vert_{2,0}^{2}\Bigg\} \\ &+\frac{1}{2}h^{2k+2}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}_{\infty, k+1}^2+\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left( 1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}. \end{align*} Proof. The proof is given in Appendix B. In particular, when Taylor–Hood elements (k = 2, s = 1) are used, i.e., the $$C^{0}$$ piecewise-quadratic velocity space $$X_{h}$$ and the $$C^{0}$$ piecewise-linear pressure space $$Q_{h}$$, we have the following estimate. Corollary 3.2 Assume that $$\Vert e_{j}^{0}\Vert$$ and $$\Vert \nabla e_j^0\Vert$$ are both $$\mathcal{O}$$(h) accurate or better. Then if $$(X_{h},Q_{h})$$ is chosen as the $$(P_2, P_1)$$ Taylor–Hood element pair, we have   \begin{equation*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} +\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\Big( 1-\frac{\sqrt{\mu}}{2}\Big)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \leqslant C (h^{2} + \varDelta t^2 +h \varDelta t). \end{equation*} 4. Numerical experiments In this section we illustrate the effectiveness of our proposed method (1.2) and the associated theoretical analyses in Sections 2 and 3 by considering two examples: a Green–Taylor vortex problem and a flow between two offset cylinders. The first problem has a known exact solution that is used to illustrate the error analysis. The second example does not have an analytic solution but has complex flow structures; it is used to check the stability analysis and demonstrate the necessity and efficiency of the proposed method. 4.1 The Green–Taylor vortex problem The Green–Taylor vortex flow is commonly used for testing convergence rates, e.g., see Chorin (1968), Tafti (1996), John & Layton (2002), Berselli (2005), Barbato et al. (2007) and Jiang & Tran (2015). The Green–Taylor vortex solution given by   \begin{align} & u(x,y,t)=-\cos(\omega\pi x)\sin(\omega\pi y)e^{-2\omega^{2}\pi^{2}t/\tau }, \nonumber\\ & v(x,y,t)=\sin(\omega\pi x)\cos(\omega\pi y)e^{-2\omega^{2}\pi^{2}t/\tau }, \\ & p(x,y,t)=-\tfrac{1}{4}(\cos(2\omega\pi x)+\cos(2\omega\pi y))e^{-4\omega ^{2}\pi^{2}t/\tau} \nonumber \end{align} (4.1)satisfies the NSE in $$\varOmega =(0,1)^{2}$$ for $$\tau =\textrm{Re}$$ and initial condition   $$u^{0}=\big(\!\!-\cos(\omega\pi x)\sin(\omega\pi y), \sin(\omega\pi x)\cos(\omega\pi y)\big)^{\textrm{T}}.$$The solution consists of an $$\omega \times \omega$$ array of oppositely signed vortices that decay as $$t\rightarrow \infty$$. In the following numerical tests we take $$\omega =1$$, $$\nu =1/\textrm{Re}$$, T = 1, h = 1/m and $$\varDelta t/h=2/5$$. The boundary condition is assumed to be inhomogeneous Dirichlet, that is, the boundary values match that of the exact solution. We consider an ensemble of two members, $$u_{1}$$ and $$u_{2}$$, corresponding to two incompressible NSE simulations with different viscosity coefficients $$\nu _j$$ and initial conditions $$u_{j, 0}$$. We investigate the ensemble simulations and compare them with independent simulations. For j = 1, 2 we define by $$\mathcal{E}^{\, E}_ {\, j} = u_{j}-u_{j, h}$$ the approximation error of the jth member of the ensemble simulation and by $$\mathcal{E}^{\, S}_{\, j} = u_{j}-u_{j, h}$$ the approximation error of the jth independently determined simulation. Here, the superscript ‘E’ stands for ‘ensemble’ whereas ‘S’ stands for ‘independent’. Case 1. We set the viscosity coefficient $$\nu _1=0.2$$ and initial condition $$u_{1, 0}= (1+\varepsilon )u^0$$ for the first member $$u_1$$ and $$\nu _2=0.3$$ and $$u_{2, 0}= (1-\varepsilon )u^0$$ for the second member $$u_2$$, where $$\varepsilon = 10^{-3}$$. For this choice of parameters, we have $$\vert \nu _j -\overline{\nu }\vert /\overline{\nu } =\frac{1}{5}$$ for both j = 1 and j = 2 so that condition (2.3) is satisfied. We first apply the ensemble algorithm; results are shown in Table 1. It is seen that the convergence rate for $$u_{1}$$ and $$u_{2}$$ is first order. Table 1 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors for ensemble simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  View Large Table 1 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors for ensemble simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  View Large We next compare the ensemble simulations with independent simulations. To this end, we perform each NSE simulation independently using the same discretization setup. The associated approximation errors are listed in Table 2. Comparing with Table 1, we observe that the ensemble simulation is able to achieve accuracies close to those of the independent, more costly simulations. Table 2 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors in independent simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  View Large Table 2 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors in independent simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  View Large Case 2. We now set $$\nu _1= 0.01$$ and $$\nu _2= 0.49$$ while keeping the same initial conditions as for Case 1. With this choice of parameters, $$\vert \nu _j -\overline{\nu }\vert /\overline{\nu } =\frac{24}{25}$$ for both j = 1 and j = 2, which still satisfies (2.3) but is closer to the upper limit. The ensemble simulation errors are listed in Table 3, which shows that the rate of convergence for the second member is nearly 1 and for the first member is approaching 1. Table 3 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in ensemble simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  View Large Table 3 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in ensemble simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  View Large The approximation errors for two independent simulations under the same discretization setup are listed in Table 4. Comparing the ensemble simulation results in Table 3 with the independent simulations, we find that the accuracy of the first member in the ensemble simulation degrades slightly whereas that of the second member in the ensemble simulation improves a bit. Overall, the ensemble simulation is able to achieve the same order of accuracy as the independent simulations. Table 4 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in independent simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  View Large Table 4 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in independent simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  View Large 4.2 Flow between two offset cylinders Next we check the stability of our algorithm by considering the problem of a flow between two offset circles (see, e.g., Jiang & Layton, 2014, 2015; Jiang, 2015; Jiang et al., 2015). The domain is a disk with a smaller off-center obstacle inside. Letting $$r_{1}=1$$, $$r_{2}=0.1$$ and $$c=(c_{1},c_{2})=\big(\frac{1}{2} ,0\big)$$, the domain is given by   $$\varOmega=\big\{(x,y)\,\,:\,\,x^{2}+y^{2}\leqslant r_{1}^{2} \,\,\textrm{and}\,\, (x-c_{1})^{2} +(y-c_{2})^{2}\geqslant r_{2}^{2}\big\}.$$ The flow is driven by a counterclockwise rotational body force   $$f(x,y,t)=\big(-6y(1-x^{2}-y^{2}), 6x(1-x^{2}-y^{2})\big)^{\textrm{T}}$$with no-slip boundary conditions imposed on both circles. The flow between the two circles shows interesting structures interacting with the inner circle. A von K$$\acute{\textrm{a}}$$rm$$\acute{\textrm{a}}$$n vortex street is formed behind the inner circle and then re-interacts with that circle and with itself, generating complex flow patterns. We consider multiple numerical simulations of the flow with different viscosity coefficients using the ensemble-based algorithm (2.1). For spatial discretization, we apply the Taylor–Hood element pair on a triangular mesh that is generated by Delaunay triangulation with 80 mesh points on the outer circle and 60 mesh points on the inner circle and with refinement near the inner circle, resulting in 18,638 degrees of freedom; see Fig. 1. In order to illustrate the stability analysis, we select two different sets of viscosity coefficients for Case 1:$$\nu _1=0.005$$, $$\nu _2=0.039$$, $$\nu _3=0.016$$; Case 2:$$\nu _1=0.005$$, $$\nu _2=0.041$$, $$\nu _3=0.014$$. Note that the viscosity coefficients in Case 2 are obtained by slightly perturbing those in Case 1. The average of the viscosity coefficients is $$\overline{\nu }= 0.02$$ for both cases. However, the stability condition (2.3) is satisfied in the first case but is not satisfied in the second one, i.e., we have Case 1:$$\frac{\vert \nu _1-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{4}<1$$, $$\frac{\vert \nu _2-\overline{\nu }\vert }{\overline{\nu }}=\frac{19}{20}<1$$, $$\frac{\vert \nu _3-\overline{\nu }\vert }{\overline{\nu }}=\frac{1}{5}<1$$; Case 2:$$\frac{\vert \nu _1-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{4}<1$$, $$\frac{\vert \nu _2-\overline{\nu }\vert }{\overline{\nu }}=\frac{21}{20}>1$$, $$\frac{\vert \nu _3-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{10}<1$$. The second member of Case 2 has a perturbation ratio greater than 1. Simulations of both cases are subject to the same initial condition and body forces for all ensemble members. In particular, the initial condition is generated by solving the steady Stokes problem with viscosity $$\nu =0.02$$ and the same body force f(x, y, t). All the simulations are run over the time interval [0, 5] with a time-step size $$\varDelta t= 0.01$$. This time-step size is smaller than the one ensuring a stable ensemble simulation in Case 1, thus condition (2.2) holds. Furthermore, the same $$\varDelta t$$ is used in Case 2. Because the viscosity coefficients in Case 2 are slightly perturbed from those in Case 1, and $$\varDelta t$$ is chosen small, we believe condition (2.2) still holds. But condition (2.3) no longer holds. Therefore, we expect the ensemble simulation to be unstable even when using the small time-step size $$\varDelta t=0.01$$. For testing the stability, we use the kinetic energy as a criterion and compare the ensemble simulation results with independent simulations using the same mesh and time-step size. Fig. 1. View largeDownload slide Mesh for the flow between two offset cylinders example. Fig. 1. View largeDownload slide Mesh for the flow between two offset cylinders example. Fig. 2. View largeDownload slide For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 2. View largeDownload slide For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 3. View largeDownload slide For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 3. View largeDownload slide For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). The comparison of the energy evolution of ensemble-based simulations with the corresponding independent simulations is shown in Figs 2 and 3. It is seen that for Case 1 the ensemble simulation is stable, but for Case 2 it becomes unstable. This phenomenon coincides with our stability analysis since condition (2.3) holds for all members of Case 1, but does not hold for the second member of Case 2. Indeed, it is observed from Fig. 3 that the energy of the second member in Case 2 blows up after t = 3.7, then affecting the other two members, and this results in their energy dramatically increasing after t = 4.7. Next we use this test example to illustrate the necessity of ensemble simulations. Indeed, the ensemble calculation is much needed when flow problems are under parameter uncertainty. This is because even though one could obtain a ‘best estimate’ of the parameter that is close enough to the true parameter value, the corresponding individual solution may vary significantly from the true solution due to the nonlinearity of the problems. However, the ensemble mean of the model problems at several selected parameter samples, which are drawn from a probability distribution (usually experimentally determined), tends to smooth out possible errors and gives a better forecast than the single run using the ‘best estimate’ parameter value. Here we consider a simple computational setting in which the true viscosity $$\nu = 0.01$$. Due to the errors in measurements, the ‘best estimate’ of the viscosity is 0.00995. Suppose the viscosity value empirically follows a uniform distribution on [0.009, 0.011]. We first run the simulation with viscosity $$\nu =0.01$$ on the domain partitioned by a Delaunay triangular mesh with 80 mesh points on the outer circle and 40 mesh points on the inner circle, take the uniform time-step size to be $$\varDelta t=0.002$$ and label the output solution ‘true’ as it is the numerical solution associated with the true datum. We then run the simulation associated to the ‘best estimate’ $$\nu =0.00995$$ and label its solution ‘best estimate’, which is to be compared with the ensemble forecast. For the latter, we draw J (= 16, 32, 64) independent, uniformly distributed samples of $$\nu$$ from [0.009, 0.011], and seek the mean of these J simulation results by using the proposed ensemble algorithm. Time evolution of the kinetic energy for the single ‘best estimate’ ensemble means of J simulations are shown in Fig. 4, together with that of the ‘true’ solution. It is observed that although the ‘best estimate’ of $$\nu$$ is very close to the true value, the difference between these two solutions is relatively large. However, the ensemble mean is able to yield a more accurate forecast. The greater the number of realizations, the better forecast that ensemble mean could achieve. Fig. 4. View largeDownload slide For the flow between two offset cylinders, comparison of ‘true’ solution, ‘best estimate’ solution and ensemble forecast. Fig. 4. View largeDownload slide For the flow between two offset cylinders, comparison of ‘true’ solution, ‘best estimate’ solution and ensemble forecast. Finally, we illustrate the efficiency of the proposed ensemble algorithm. To this end, we solve a group of J flow problems using a sequence of individual simulations and one ensemble simulation, respectively. These problems are subject to fixed initial and boundary conditions and a forcing function, but different viscosity parameter values. We then compare the efficiency of these two approaches by considering the CPU time. Assume the viscosity parameter $$\nu _j$$ to be a random variable that distributes uniformly on the interval [0.4, 0.5]. We generate a set of J (= 2, 4, 8, 16) random samples using MATLAB, partition the domain by a Delaunay triangular mesh with 80 mesh points on the outer circle and 40 mesh points on the inner circle and take the uniform time-step size to be $$\varDelta t=0.002$$. Both individual and ensemble simulations are implemented by using FreeFem++ with UMFPACK, which solves the linear system with multifrontal LU factorization. Numerical results together with the CPU time for simulations (measured in seconds) are listed in Table 5, where $$\Vert \overline{u}^{E}(T)\Vert$$ is the $$L_2$$ norm of the mean of the velocity obtained by the ensemble algorithm (2.1) at the final simulation time T = 0.1 and $$\Vert \overline{u}^{S}(T)\Vert$$ is that of individual simulations. Table 5 Efficiency comparison between the ensemble simulations and individual simulations on a set of J parametrized NSE flow between two offset cylinders problems at T = 0.1. It is observed that the ensemble algorithm outperforms the individual simulations as the former spends less CPU times than the latter but keeps the same accuracy J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  View Large Table 5 Efficiency comparison between the ensemble simulations and individual simulations on a set of J parametrized NSE flow between two offset cylinders problems at T = 0.1. It is observed that the ensemble algorithm outperforms the individual simulations as the former spends less CPU times than the latter but keeps the same accuracy J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  View Large Define the speedup factor of the ensemble algorithm to be the ratio of execution time for the sequence of individual simulations to that of the ensemble simulation. It is observed that for this simple test, the computational time for completing the sequence of individual simulations increases almost linearly as the size of the ensemble increases, which is expected as all the work including assembly of matrices and LU decomposition needs to be done repeatedly for each ensemble member. However, the ensemble algorithm takes significantly less CPU time while maintaining the same accuracy. It is seen that for J = 2 the speedup factor is 2.30, while for J = 16 the speedup factor increases to 4.81. This saving rate will keep increasing as the size of the ensemble increases. 5. Conclusions In this paper, we consider a set of Navier–Stokes simulations in which each member may be subject to a distinct viscosity coefficient, initial conditions and/or body forces. An ensemble algorithm is developed for the group by which all the flow ensemble members, after discretization, share a common linear system with different right-hand-side vectors. This leads to a great saving in both storage requirements and computational costs. The stability and accuracy of the ensemble method are analysed. Two numerical experiments are presented. The first is for Green–Taylor flow and serves to illustrate the first-order accuracy in time of the ensemble-based scheme. The second is for a flow between two offset cylinders and serves to show that our stability analysis is sharp and the proposed method is effective. As a next step, we will investigate higher-order accurate schemes for the flow ensemble simulations. Funding US Air Force Office of Scientific Research (FA9550-15-1-0001); US Department of Energy Office of Science (DE-SC0009324, DE-SC0016591, DE-SC0016540); US National Science Foundation (DMS-1522672, DMS-1720001); University of Missouri Research Board. References Barbato, D., Berselli, L. & Grisanti, C. ( 2007) Analytical and numerical results for the rational large eddy simulation model. J. Math. Fluid Mech. , 9, 44– 74. Google Scholar CrossRef Search ADS   Berselli, L. ( 2005) On the large eddy simulation of the Taylor-Green vortex. J. Math. Fluid Mech. , 7, S164– S191. Google Scholar CrossRef Search ADS   Brenner, S. & Scott, R. ( 2008) The Mathematical Theory of Finite Element Methods , 3rd edn. New York: Springer-Verlag Google Scholar CrossRef Search ADS   Chorin, A. ( 1968) Numerical solution for the Navier–Stokes equations. Math. Comp. , 22, 745– 762. Google Scholar CrossRef Search ADS   Ciarlet, P. ( 2002) The Finite Element Method for Elliptic Problems . Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Fiordilino, J. ( 2017a) A second order ensemble timestepping algorithm for natural convection. Preprint, arXiv:1708.00488. Fiordilino, J. ( 2017b) Ensemble timestepping algorithms for the heat equation with uncertain conductivity. Preprint, arXiv:1708.00893. Girault, V. & Raviart, P.-A. ( 1979) Finite Element Approximation of the Navier-Stokes Equations . Lecture Notes in Mathematics , vol. 749. Berlin: Springer. Google Scholar CrossRef Search ADS   Girault, V. & Raviart, P.-A. ( 1986) Finite Element Methods for Navier-Stokes Equations - Theory and Algorithms . Berlin: Springer. Google Scholar CrossRef Search ADS   Gunzburger, M. ( 1989) Finite Element Methods for Viscous Incompressible Flows - A Guide to Theory, Practice, and Algorithms.  London: Academic Press. Gunzburger, M., Jiang, N. & Schneier, M. ( 2017a) An ensemble-proper orthogonal decomposition method for the nonstationary Navier-Stokes equations. SIAM J. Numer. Anal. , 55, 286– 304. Google Scholar CrossRef Search ADS   Gunzburger, M., Jiang, N. & Schneier, M. ( 2018) A higher-order ensemble/proper orthogonal decomposition method for the nonstationary Navier-Stokes equations. Int. J. Numer. Anal. Model. , 15, 608– 627. Gunzburger, M., Jiang, N. & Wang, Z. ( 2017b) A second-order time-stepping scheme for simulating ensembles of parameterized flow problems. Comput. Meth. Appl. Math.  (in press). https://doi.org/10.1515/cmam-2017-0051. Gutknecht, M. ( 2007) Block Krylov space methods for linear systems with multiple right-hand sides: an introduction. Modern Mathematical Models, Methods and Algorithms for Real World Systems (A. H. Siddiqi, I. S. Duff and O. Christensen eds). New Delhi: Anamaya Publishers, pp. 420– 447. Jiang, N. ( 2015) A higher order ensemble simulation algorithm for fluid flows. J. Sci. Comput. , 64, 264– 288. Google Scholar CrossRef Search ADS   Jiang, N. ( 2017) A second-order ensemble method based on a blended backward differentiation formula timestepping scheme for time-dependent Navier-Stokes equations. Numer. Meth. Partial. Diff. Eqs. , 33, 34– 61. Google Scholar CrossRef Search ADS   Jiang, N. & Layton, W. ( 2014) An algorithm for fast calculation of flow ensembles. Int. J. Uncertain. Quantif. , 4, 273– 301. Google Scholar CrossRef Search ADS   Jiang, N. & Layton, W. ( 2015) Numerical analysis of two ensemble eddy viscosity numerical regularizations of fluid motion. Numer. Meth. Partial Diff. Eqs. , 31, 630– 651. Google Scholar CrossRef Search ADS   Jiang, N. & Tran, H. ( 2015) Analysis of a stabilized CNLF method with fast slow wave splittings for flow problems. Comput. Meth. Appl. Math. , 15, 307– 330. Jiang, N., Kaya, S. & Layton, W. ( 2015) Analysis of model variance for ensemble based turbulence modeling. Comput. Meth. Appl. Math. , 15, 173– 188. John, V. & Layton, W. ( 2002) Analysis of numerical errors in large eddy simulation. SIAM J. Numer. Anal. , 40, 995– 1020. Google Scholar CrossRef Search ADS   Layton, W. ( 2008) Introduction to the Numerical Analysis of Incompressible Viscous Flows . Philadelphia: Society for Industrial and Applied Mathematics (SIAM) . Google Scholar CrossRef Search ADS   Luo, Y. & Wang, Z. ( 2018a) An ensemble algorithm for numerical solutions to deterministic and random parabolic PDEs. SIAM J. Numer. Anal. , 56, 859– 876. Google Scholar CrossRef Search ADS   Luo, Y. & Wang, Z. ( 2018b) A multilevel Monte Carlo ensemble scheme for solving random parabolic PDEs. Preprint, arXiv:1802.05743. Mohebujjaman, M. & Rebholz, L. ( 2017) An efficient algorithm for computation of MHD flow ensembles. Comput. Meth. Appl. Math. , 17, 121– 138. Google Scholar CrossRef Search ADS   Parks, M., Soodhalter, K. M. & Szyld, D. B. ( 2016) A block recycled GMRES method with investigations into aspects of solver performance. Preprint, arXiv:1604.01713. Tafti, D. ( 1996) Comparison of some upwind-biased high-order formulations with a second order central-difference scheme for time integration of the incompressible Navier-Stokes equations. Comput. & Fluids , 25, 647– 665. Google Scholar CrossRef Search ADS   Takhirov, A., Neda, M. & Waters, J. ( 2016) Time relaxation algorithm for flow ensembles. Numer. Meth. Partial Diff. Eqs. , 32, 757– 777. Google Scholar CrossRef Search ADS   Thomée, V. ( 1997) Galerkin Finite Element Methods for Parabolic Problems . Berlin, Heidelberg: Springer. Google Scholar CrossRef Search ADS   Appendix A. Proof of Theorem 2.1 Setting $$v_{h}=u_{j,h}^{n+1}$$ and $$q_h = p_{j, h}^{n+1}$$ in (2.1) and then adding two equations we obtain   \begin{align*} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert ^{2}+\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2}+ \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n}, u_{j,h}^{n},u_{j,h}^{n+1}) \\ &+ \varDelta t \Vert \overline{\nu}^{\frac{1}{2}} \nabla u_{j,h}^{n+1}\Vert^{2}= \varDelta t (f_{j}^{n+1}, u_{j,h}^{n+1}) - \varDelta t \left( \left( \nu_j - \overline{\nu}\right)\nabla u_{j,h}^n, \nabla u_{j,h}^{n+1}\right). \end{align*}Note that $$\overline{\nu }\geqslant \overline{\nu }_{\textrm{min}}>0$$. We apply Young’s inequality to the terms on the RHS and have, for all $$\alpha , \beta>0$$,   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} + \overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla u_{j,h}^{n+1}\Vert^{2} + \varDelta t b^{*} \left(u_{j,h}^{n} -\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1} -u_{j,h}^{n} \right)\nonumber\\ \leqslant &\frac{\alpha\overline{\nu}_{\textrm{min}}\varDelta t}{8} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} + \frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2} +\frac{\beta\overline{\nu}_{\textrm{min}}\varDelta t}{4} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} +\frac{\|\nu_j-\overline{\nu}\|^2_{\infty}\varDelta t}{\beta\overline{\nu}_{\textrm{min}}}\Vert \nabla u_{j,h}^n\Vert^2. \end{align} (A1)Both $$\frac{\beta \overline{\nu }_{\textrm{min}}\varDelta t}{4} \Vert \nabla u_{j,h}^{n+1} \Vert ^{2}$$ and $$\frac{\|\nu _j-\overline{\nu }\|_{\infty }^2\varDelta t}{\beta \overline{\nu }_{\textrm{min}}}\Vert \nabla u_{j,h}^n\Vert ^2$$ on the RHS of (A1) need to be absorbed into $$\overline{\nu }_{\textrm{min}}\varDelta t \Vert \nabla u_{j,h}^{n+1}\Vert ^{2}$$ on the LHS. To this end, we minimize $$\frac{\beta \overline{\nu }_{\textrm{min}}\varDelta t}{4}+\frac{\|\nu _j-\overline{\nu }\|_{\infty }^2\varDelta t}{\beta \overline{\nu }_{\textrm{min}}}$$ by selecting $$\beta = \frac{2 \| \nu _j-\overline{\nu }\|_{\infty } }{\overline{\nu }_{\textrm{min}}}$$ so that (A1) becomes   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} + \overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla u_{j,h}^{n+1}\Vert^{2} + \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})\nonumber\\ \leqslant&\, \frac{\alpha\overline{\nu}_{\textrm{min}}\varDelta t}{8} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} + \frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert f_{j}^{n+1}\Vert_{-1}^{2} +\frac{\| \nu_j-\overline{\nu}\|_{\infty}\varDelta t}{2} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} +\frac{\| \nu_j-\overline{\nu}\|_{\infty} \varDelta t}{2}\Vert \nabla u_{j,h}^n\Vert^2. \end{align} (A2)Next we bound the trilinear term using the inequality (1.9) and the inverse inequality (1.7), obtaining   \begin{align*} - \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})&\leqslant C\varDelta t\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert\Vert\nabla u_{j,h} ^{n}\Vert\left[\Vert \nabla ( u_{j,h}^{n+1}-u_{j,h}^{n})\Vert\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert\right]^{1/2}\\ &\leqslant C\varDelta t\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert\Vert\nabla u_{j,h} ^{n}\Vert Ch^{-\frac{1}{2}} \Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert. \end{align*}Using Young’s inequality again gives   $$- \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})\leqslant C \frac{\varDelta t^{2}}{h} \Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2} \Vert\nabla u_{j,h}^{n}\Vert^{2}+ \frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h} ^{n}\Vert^{2}.$$ (A3)Substituting (A3) into (A2) and combining like terms, we have   \begin{align*} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\\&+\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2} \leqslant\frac{2\varDelta t}{\alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2} + C \frac{\varDelta t^{2}}{h} \Vert\nabla( u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2} \Vert\nabla u_{j,h}^{n}\Vert^{2} \\&+\frac{\| \nu_j-\overline{\nu}\|_{\infty} \varDelta t}{2}\Vert \nabla u_{j,h}^n\Vert^2. \end{align*}For any $$0<\sigma <1$$,   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+ \overline{\nu}_{\textrm{min}}\varDelta t \Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)\left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert^{2}\right)\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \left[(1-\sigma)\Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla u_{j,h}^{n}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}} \varDelta t \left[\sigma\Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla u_{j,h}^n\Vert^2 \leqslant\frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert f_{j}^{n+1}\Vert_{-1}^{2}\text{ .} \end{align} (A4)Select $$\alpha =4-\frac{2(\sigma +1)}{\sigma }\sqrt{\mu }$$. Since $$\alpha$$ is supposed to be greater than 0, we have   $$\sigma>\frac{\sqrt{\mu}}{2-\sqrt{\mu}} \in (0,1).$$ (A5)Now taking $$\sigma = \frac{\sqrt{\mu }+\varepsilon }{2-\sqrt{\mu }}$$, where $$\varepsilon \in (0, 2-2\sqrt{\mu })$$, (A4) becomes   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+ \overline{\nu}_{\textrm{min}}\varDelta t \Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)\left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \bigg[\Big( (1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu}-(1-\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big) -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\bigg] \Vert\nabla u_{j,h}^{n}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}} \varDelta t \Big[ \frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\Big]\Vert \nabla u_{j,h}^n\Vert^2\leqslant\frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert _{-1}^{2}\text{ .} \end{align} (A6)Stability follows if the following conditions hold:   $$(1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu} -(1-\sigma)\frac{1}{2}\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{\overline{\nu}_{\textrm{min}}} -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \geqslant0,$$ (A7)  $$\frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\geqslant 0.$$ (A8)Using assumption (2.3), we have   \begin{equation*} \frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}} = \frac{2+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}-\frac{\|\nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\right) \geqslant 0 \end{equation*}so that (A7) holds. Together with assumption (2.2), we then have   \begin{align*} &(1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu} -(1-\sigma)\frac{1}{2}\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{\overline{\nu}_{\textrm{min}}} -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \\ &\qquad\geqslant (1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu}-(1-\sigma)\frac{1}{2}\sqrt{\mu}-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \\ &\qquad=\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}=0, \end{align*}so (A8) holds. Therefore, assuming that (2.2) and (2.3) hold, (A6) reduces to   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert^{2}\right) \leqslant \frac{2\varDelta t}{ \overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2}\text{ .} \end{align} (A9)Summing (A9) from n = 0 to n = N − 1 results in   \begin{align} \frac{1}{2}\|u_{j,h}^{N}\|^{2} &+\frac{1}{4}\sum_{n=0}^{N-1}\|u_{j,h}^{n+1}-u_{j,h}^{n}\|^{2} +\sum_{n=0}^{N-1}\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ & +\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon} -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \|\nabla u_{j,h}^{N}\|^{2}\nonumber\\ \leqslant& \sum_{n=0}^{N-1}\frac{2\varDelta t}{\overline{\nu}_{\textrm{min}}}\|f_{j}^{n+1}\|_{-1}^{2}+ \frac{1} {2}\|u_{j,h}^{0}\|^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon} -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\|\nabla u_{j,h}^{0}\|^{2}. \end{align} (A10)This completes the proof of stability. Appendix B. Proof of Theorem 3.1 The weak solution of the NSE $$u_{j}$$ satisfies   \begin{align} \left(\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}, v_{h}\right)&+ b^{*}(u_{j}^{n+1}, u_{j}^{n+1}, v_{h})+ (\nu_j\nabla u_{j}^{n+1}, \nabla v_{h})\nonumber\\ &- (p_{j} ^{n+1},\nabla\cdot v_{h})=(\,f_{j}^{n+1}, v_{h}) + \textrm{Intp}(u_{j}^{n+1};v_{h})\qquad\forall\, v_{h}\in V_{h}, \end{align} (B1)where $$\textrm{Intp}(u_{j}^{n+1};v_{h})= \left (\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t} (t^{n+1}),v_{h}\right)$$. Let   $$e_{j}^{n}=u_{j}^{n}-u_{j,h}^{n}=(u_{j}^{n}-R_{h} u_{j}^{n})+(R_{h} u_{j} ^{n}-u_{j,h}^{n})=\eta_{j}^{n}+\xi_{j,h}^{n},$$where $$R_{h} u_{j}^{n} \in V_{h}$$ is the Ritz projection of $$u_{j}^{n}$$ onto $$V_{h}$$ defined by $$\big (\overline{\nu }~\nabla (u_{j}^{n}-R_{h} u_{j}^{n}), \nabla v_h\big ) = 0$$. The associated Ritz projection error satisfies (see, e.g., Thomée, 1997)   $$\|R_h u_j^n - u_j^n\| + h \|\nabla (R_h u_j^n - u_j^n)\| \leqslant C h^{k+1} \|u_j^n\|_{k+1}.$$ (B2)Subtracting (3.1) from (B1) gives   \begin{align*} \left(\frac{\xi_{j,h}^{n+1}-\xi_{j,h}^{n}}{ \varDelta t},v_{h}\right) &+(\overline{\nu}\nabla\xi_{j,h}^{n+1},\nabla v_{h}) +((\nu_j-\overline{\nu})\nabla (u_{j}^{n+1}-u_{j}^n),\nabla v_{h}) +((\nu_j-\overline{\nu})\nabla \xi_{j,h}^n,\nabla v_h) \nonumber\\ & +b^{*}(u_{j}^{n+1},u_{j}^{n+1},v_{h})-b^{*}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h}) -b^{*}(u_{j,h}^{n}-\overline{u}_{h}^n,u_{j,h} ^{n},v_{h})-(p_{j}^{n+1},\nabla\cdot v_{h}) \nonumber\\ =&-\left(\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{\varDelta t},v_{h}\right) -(\overline{\nu}\nabla\eta_{j}^{n+1},\nabla v_{h})-((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla v_{h}) +\textrm{Intp}(u_j^{n+1};v_{h}).\nonumber \end{align*}Setting $$v_{h}=\xi _{j,h}^{n+1}\in V_{h}$$ and rearranging the nonlinear terms, we have   $$\begin{split} &\frac{1}{\varDelta t}\left(\frac{1}{2}\|\xi_{j,h}^{n+1}\|^{2}-\frac{1}{2}\|\xi _{j,h}^{n}\|^{2}+\frac{1}{2}\|\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\|^{2}\right)+ \|\overline{\nu}^{\frac{1}{2}}\nabla\xi_{j,h}^{n+1}\|^{2}\\ &\qquad=-\,((\nu_j-\overline{\nu})\nabla (u_{j}^{n+1}-u_{j}^n),\nabla \xi_{j,h}^{n+1}) -((\nu_j-\overline{\nu})\nabla \xi_{j,h}^n,\nabla \xi_{j,h}^{n+1}) -((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla \xi_{j,h}^{n+1}) \\ &\qquad\quad-b^{*}(u_{j,h}^n-\overline{u}_h^n,u_{j,h}^{n+1}-u_{j,h}^{n},\xi_{j,h}^{n+1}) -b^{*}(u_{j}^{n+1},u_{j}^{n+1},\xi_{j,h}^{n+1})+b^{*}(u_{j,h}^{n}, u_{j,h}^{n+1},\xi_{j,h}^{n+1})\\ &\qquad\quad+(p_{j}^{n+1},\nabla\cdot\xi_{j,h}^{n+1})-\left(\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t},\xi_{j,h}^{n+1}\right) +\textrm{Intp}(u_j^{n+1};\xi_{j,h}^{n+1}). \end{split}$$ (B3)We first bound the viscous terms on the RHS of (B3):   \begin{align} -((\nu_j-\overline{\nu}) \nabla (u_{j}^{n+1}-u_j^n), \nabla\xi_{j,h}^{n+1}) &\leqslant \| \nu_j-\overline{\nu}\|_{\infty} \|\nabla (u_j^{n+1}-u_j^n)\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\|\nabla(u_{j}^{n+1}-u_j^n)\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber \\ &\leqslant \frac{\varDelta t}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \left(\int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2}\, \textrm{d}t\right)+ C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B4)  \begin{align} -((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla\xi_{j,h}^{n+1}) &\leqslant\| \nu_j-\overline{\nu}\|_{\infty} \|\nabla\eta_{j} ^{n}\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B5)  \begin{align} -(\nu_j-\overline{\nu})(\nabla\xi_{j,h}^{n},\nabla\xi_{j,h}^{n+1}) &\leqslant\| \nu_j-\overline{\nu}\|_{\infty} \|\nabla\xi_{j,h}^{n}\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_1}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \|\nabla\xi_{j,h}^{n}\|^{2} + C_1\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber \\ &\leqslant \frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2} \|\nabla\xi_{j,h}^{n}\|^{2}+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2} \|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B6)in which we note that both terms on the RHS of (B6) need to be hidden in the LHS of the error equation, thus $$C_1= \frac{\|\nu _j-\overline{\nu }\|_{\infty }}{2\overline{\nu }_{\textrm{min}}}$$ is selected to minimize the summation. Next we analyse the nonlinear terms on the RHS of (B3) one by one. For the first nonlinear term we have   \begin{align}\!\! -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j,h}^{n+1}-u_{j,h}^n, \xi_{j,h}^{n+1}) =&-\,b^{\ast}(u_{j,h}^n-\overline{u}_h^n, e_{j}^{n+1}-e_{j}^n, \xi_{j,h}^{n+1}) +b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j}^{n+1}-u_{j}^n, \xi_{j,h}^{n+1})\nonumber\\ =&-\,b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n+1}, \xi_{j,h}^{n+1})+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n}, \xi_{j,h}^{n+1})\nonumber\\ &+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \xi_{j}^{n}, \xi_{j,h}^{n+1})+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j}^{n+1}-u_{j}^n, \xi_{j,h}^{n+1})\,. \end{align} (B7)Using inequality (1.8) and Young’s inequality, we have the estimates   \begin{align} -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla\eta_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla\eta_{j}^{n+1}\|^{2} +C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B8)and   \begin{align} -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla\eta_{j}^{n}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1} \|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B9)Because $$b^{\ast }(\cdot ,\cdot ,\cdot )$$ is skew-symmetric we have   \begin{align*} b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n},\xi_{j,h}^{n+1})&=b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi _{j,h}^{n}-\xi_{j,h}^{n+1},\xi_{j,h}^{n+1}) \\& =b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n+1},\xi_{j,h}^{n+1}-\xi_{j,h}^n)\,. \end{align*}Then by inequality (1.8) we obtain   \begin{align} b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n},\xi_{j,h}^{n+1}) &\leqslant \Vert \nabla (u_{j,h}^n-\overline{u}_h^n)\Vert\Vert \nabla \xi_{j,h}^{n}\Vert \Vert \nabla (\xi_{j,h}^{n+1}-\xi_{j,h}^n)\Vert^{1/2}\Vert \xi_{j,h}^{n+1}-\xi_{j,h}^n\Vert^{1/2}\nonumber\\ &\leqslant C\Vert \nabla (u_{j,h}^n-\overline{u}_h^n)\Vert\Vert \nabla \xi_{j,h}^{n}\Vert h^{-1/2}\Vert \xi_{j,h}^{n+1}-\xi_{j,h}^n\Vert\nonumber\\ &\leqslant \frac{1}{4\varDelta t}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2}+\left[ C\frac{\varDelta t}{h}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}. \end{align} (B10)For the last term of (B7) we have   \begin{align} b^{*}(u_{j,h}^n-\overline{u}_h^n,u_{j}^{n+1}-u_{j}^{n},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla(u_{j}^{n+1}-u_{j}^{n})\|\|\nabla\xi_{j,h}^{n+1} \|\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla (u_{j}^{n+1}-u_{j}^{n})\|^{2} +C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &\leqslant \frac{C^2}{4C_0}\varDelta t \overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\left(\int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \, \mathrm{d}t\right) +C_0\overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} .\nonumber \end{align} (B11)Next we bound the last two nonlinear terms on the RHS of (B3) as   \begin{align} &-b^{*}(u_{j}^{n+1}, u_{j}^{n+1},\xi_{j,h}^{n+1}) +b^{*}(u_{j,h}^{n},u_{j,h}^{n+1},\xi_{j,h}^{n+1})\nonumber\\ &\qquad= -b^{*}(e_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u_{j,h}^{n},e_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u^{n+1}_{j}-u_j^{n}, u_{j}^{n+1}, \xi_{j,h}^{n+1}) \nonumber\\ &\qquad= -b^{*}(\eta^{n}_j,u_{j}^{n+1},\xi_{j,h}^{n+1})-b^{*}(\xi_{j,h}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u^{n}_{j,h}, \eta_{j}^{n+1}, \xi_{j,h}^{n+1})-b^{*}(u^{n+1}_{j}-u^{n}_j, u_{j}^{n+1}, \xi_{j,h}^{n+1}), \end{align} (B12)where, with the assumption $$u_j^{n+1}\in L^{\infty }(0,T; H^1(\varOmega ))$$, we have   \begin{align} -b^{*}(\eta_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla \eta_{j}^{n}\|\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B13)Using inequality (1.8), Young’s inequality and $$u_j^{n+1}\in L^{\infty }(0,T; H^1(\varOmega ))$$ we get   \begin{align} -b^{*}(\xi^{n}_{j,h}, u_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\| \nabla\xi ^{n}_{j,h} \|^{1/2} \Vert\xi^{n}_{j,h}\Vert^{1/2}\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\leqslant C\| \nabla\xi^{n}_{j,h} \|^{1/2} \Vert\xi^{n}_{j,h}\Vert ^{1/2}\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant C\Big(\frac{1}{4\alpha}\|\nabla\xi^{n}_{j,h} \|\|\xi^{n}_{j,h} \| + \alpha\|\nabla\xi_{j,h}^{n+1}\|^{2}\Big)\\ &\leqslant C\Big[\frac{1}{4\alpha} \big(\frac{\varDelta}{2}\|\nabla\xi^{n}_{j,h} \|^{2}+\frac{1}{2\varDelta}\|\xi^{n}_{j,h} \|^2 \big) + \alpha\|\nabla\xi_{j,h}^{n+1}\|^{2} \Big]\nonumber\\ &\leqslant C_0\overline{\nu}_{\textrm{min}} \|\nabla \xi^{n}_{j,h} \|^{2}+\frac{C^4}{64C_0^3\overline{\nu}_{\textrm{min}}^{3}}\|\xi^{n}_{j,h} \|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2}\,,\nonumber \end{align} (B14)where we set $$\alpha = \frac{C_0\overline{\nu }_{\textrm{min}}}{C}$$ and $$\varDelta = \frac{8C_0^2\overline{\nu }_{\textrm{min}}^2}{C^2}$$. By Young’s inequality, (1.8) and the result (A10) from the stability analysis, i.e., $$\Vert u_{j,h}^n \Vert ^2 \leqslant C$$, we also have   \begin{align} b^{*}(u^{n}_{j,h}, \eta_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla u^{n}_{j,h}\|^{1/2}\Vert u_{j,h}^n\Vert^{1/2}\|\nabla \eta_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla u_{j,h}^{n}\|\|\nabla\eta^{n+1}_{j}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B15)and   \begin{align} b^{*}(u_{j}^{n+1}-u_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j}^{n+1}-u_{j}^{n})\|\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h} ^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla(u_{j}^{n+1}-u_{j}^{n})\|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &= \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1} \varDelta t^2 \left\|\frac{\nabla u_{j}^{n+1}- \nabla u_{j}^{n}}{\varDelta t}\right\|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} \\ &= \frac{C^2}{4C_0} \overline{\nu}_{\min}^{-1} \varDelta t^2 \int_{\varOmega}\left(\frac{1}{\varDelta t}\int_{t^{n}}^{t^{n+1}} \nabla u_{j,t} \, \mathrm{d}t\right)^{2} \textrm{d}\varOmega+ C_0 \overline{\nu} _{\min}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t} \|^{2}\, \mathrm{d}t + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2}. \nonumber \end{align} (B16)For the pressure term in (B3), because $$\xi _{j,h}^{n+1}\in V_{h}$$ we have, for any $$q_{j, h}^{n+1}\in Q_h$$,   \begin{align} (p_{j}^{n+1},\nabla\cdot\xi_{j,h}^{n+1})&=(p_{j}^{n+1}-q_{j,h}^{n+1}, \nabla\cdot\xi_{j,h}^{n+1})\leqslant\sqrt{d}\,\|p_{j}^{n+1}-q_{j,h}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\\ &\leqslant \frac{1}{4d\, C_0}\overline{\nu}_{\textrm{min}}^{-1} \|p_{j}^{n+1}-q_{j,h}^{n+1}\|^{2} + C_0\,\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \text{ .}\nonumber \end{align} (B17) The other terms are bounded as   \begin{align} \Big( \frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t},\xi_{j,h}^{n+1} \Big) &\leqslant C \Big\|\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t} \Big\| \|\nabla\xi_{j,h} ^{n+1}\|\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t} \right\|^{2} +C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{1}{\varDelta t}\int_{t^{n}}^{t^{n+1}} \eta_{j,t} \textrm{ } \mathrm{d}t \right \|^2+C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\| \eta_{j,t}\|^{2}\textrm{ } \mathrm{d}t+C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B18)and   \begin{align} \textrm{Intp}(u_{j}^{n+1};\xi_{j,h}^{n+1}) &=\left(\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1}),\xi_{j,h}^{n+1}\right)\leqslant C\left\|\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1})\right\| \|\nabla \xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1}) \right\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber \\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\|u_{j,tt}\|^{2}\, \mathrm{d}t + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B19)Combining (B4)–(B19), we have   \begin{align}&\frac{1}{\varDelta t}\left(\frac{1}{2}||\xi_{j,h}^{n+1}||^{2}-\frac{1}{2}||\xi _{j,h}^{n}||^{2}+\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \right)\nonumber\\ &\quad+ C_0 \overline{\nu}_{\textrm{min}} \Vert \nabla \xi_{j,h}^{n+1}\Vert^2+ \overline{\nu}_{\textrm{min}} \left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\left(\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2}-\Vert\nabla \xi_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[(1-\sigma)\left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\sigma\left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla \xi_{j,h}^n\Vert^2\nonumber\\ &\leqslant C \bigg[ \frac{1}{\overline{\nu}_{\textrm{min}}^{3}}\Vert\xi_{j,h}^{n}\Vert^{2} +\varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \mathrm{d}t \\ &\quad+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2} \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1} }\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla u_{j,h}^{n}\Vert\Vert\nabla\eta_{j}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\mathrm{d}t \nonumber \bigg]. \end{align} (B20)Note that the generic constant C independent of $$\varDelta t$$ is used on the RHS. It, however, depends on the geometry and mesh due to the use of an inverse inequality in (B10). Similarly to the stability analysis we take $$C_0= \frac{1}{14}(1-\frac{\sigma +1}{2\sigma }\sqrt{\mu }) = \frac{1}{14} \frac{\varepsilon }{\sqrt{\mu }+\varepsilon }( 1-\frac{\sqrt{\mu }}{2})$$ with $$\sigma =\frac{\sqrt{\mu }+\varepsilon }{2-\sqrt{\mu }}$$. Then (B20) becomes   \begin{align}& \frac{1}{\varDelta t}\left(\frac{1}{2}||\xi_{j,h}^{n+1}||^{2}-\frac{1}{2}||\xi _{j,h}^{n}||^{2}+\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \right)\nonumber\\&\quad+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right) \overline{\nu}_{\textrm{min}} \Vert \nabla \xi_{j,h}^{n+1}\Vert^2 \nonumber\\&\quad+ \overline{\nu}_{\textrm{min}} \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\left(\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2}-\Vert\nabla \xi_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla \xi_{j,h}^n\Vert^2\nonumber\\ \leqslant&\ C\bigg[ \frac{1}{\overline{\nu}_{\textrm{min}}^{3}}\Vert\xi_{j,h}^{n}\Vert^{2} +\varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \textrm{ }\mathrm{d}t +\overline{\nu}_{\textrm{min}}\Vert\nabla\eta_{j}^{n+1}\Vert^{2} \\ \quad&\quad+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2} \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1} }\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla u_{j,h}^{n}\Vert\Vert\nabla\eta_{j}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\mathrm{d}t\nonumber \bigg]. \end{align} (B21)By the convergence condition (3.3) we have   \begin{equation*} \frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\geqslant \frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\sqrt{\mu}}{2}> 0 \end{equation*}and   \begin{align*}& \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\\ &\qquad\geqslant \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\sqrt{\mu}}{2}\right)-\frac{\sqrt{\mu}}{2}\\ &\qquad\geqslant \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2-\sqrt{\mu}}{\sqrt{\mu}+\varepsilon}\right)-\frac{\sqrt{\mu}}{2} \geqslant \frac{\sqrt{\mu}}{2}-\frac{\sqrt{\mu}}{2}=0 \end{align*}and by the convergence conditions (3.2) and (3.3) we have   \begin{align*} &\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}- \overline{u}_{h}^{n})\Vert^{2}\\ &\qquad\geqslant\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2-\sqrt{\mu}}{\sqrt{\mu}+\varepsilon}\right)-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}\\ &\qquad\geqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}\geq 0. \end{align*}Summing (B20) from n = 1 to N − 1 and multiplying both sides by $$\varDelta t$$ gives   \begin{align*}& \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}+\frac{1}{4}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n+1} -\xi_{j,h}^{n}\Vert^{2}\nonumber\\&\quad+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2} \\ &\quad+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{N}\|^{2}\\ \leqslant&\ \frac{1}{2}\|\xi_{j,h}^{0}\|^{2}+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} +\frac{C\varDelta t}{\overline{\nu}_{\textrm{min}}^{3}}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n}\Vert^{2}\\ &\quad+C\varDelta t\sum_{n=0}^{N-1} \bigg[ \varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \,\mathrm{d}t +\overline{\nu}_{\textrm{min}}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\\ &\qquad\qquad\qquad+\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2}\\ &\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \\ \quad&\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\,\mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert \\ \quad&\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} + \overline{\nu}_{\textrm{min}}^{-1} \varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\, \mathrm{d}t \bigg] \text{ .} \end{align*}Using the Ritz projection error (B2) and the result (2.4) from the stability analysis, i.e., $$\varDelta t \sum _{n=0}^{N-1}\Vert \nabla u_{j,h}^{n} \Vert ^2 \leqslant C$$, we have   $$C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t\sum_{n=0}^{N-1}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert \leqslant C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}\varDelta t \sum_{n=0}^{N-1}\Vert u_j^{n+1}\Vert^2_{k+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert$$ (B22)  \begin{align} \leqslant &\ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}\left( \varDelta t \sum_{n=0}^{N-1}\Vert u_{j}^{n+1}\Vert_{k+1}^4+\varDelta t \sum_{n=0}^{N-1}\Vert\nabla u_{j,h}^{n}\Vert^2 \right)\\ \leqslant &\ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}.\nonumber \end{align} (B23)By applying again the Ritz projection error (B2) and the interpolation error (1.6) we have   \begin{align} \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}&+ \overline{\nu}_{\textrm{min}}\varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}_{\textrm{min}}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{N}\|^{2} +\sum_{n=0}^{N-1}\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \nonumber\\ & +\frac{1}{15}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ \leqslant\ &\frac{1}{2}\|\xi_{j,h}^{0}\|^{2}+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} +\frac{C\varDelta t}{\overline{\nu}_{\textrm{min}}^{3}}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n}\Vert^{2} \\ &+ C\varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla u_{j,t}\Vert^2_{2,0} +C\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber \\ &+Ch^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t \Vert\nabla u_{j,t}\Vert^2_{2,0}+ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ & +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2\Vert\nabla u_{j,t}\Vert^2_{2,0} + C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \nonumber\\ &+C\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2}+C\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}\Vert u_{j,t}\Vert_{2,k+1}^{2} +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}\Vert u_{j,tt}\Vert_{2,0}^{2}.\nonumber \end{align} (B24)The next step is the application of the discrete Gronwall inequality (see Girault & Raviart, 1979, p. 176):   \begin{align} \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}&+ \overline{\nu}_{\textrm{min}}\varDelta t\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \|\nabla\xi_{j,h}^{N}\|^{2} +\sum_{n=0}^{N-1}\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \nonumber\\ &+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2}\\ \leqslant& \ e^{\frac{CT}{\overline{\nu}_{\textrm{min}}^{3}}} \bigg\{\frac{1}{2}\|\xi_{j,h}^{0}\|^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} \nonumber\\ &\qquad\quad+C\Big[ \varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}{\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} +\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} +\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ &\qquad\quad+h^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t {\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} + \overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ & \qquad\quad+\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2{\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} + \overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \nonumber \\ &\qquad\quad+\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2} +\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}{\left\vert\!\left\vert\!\left\vert u_{j,t} \right\vert\!\right\vert\!\right\vert}_{2,k+1}^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}{\left\vert\!\left\vert\!\left\vert u_{j,tt} \right\vert\!\right\vert\!\right\vert}_{2,0}^{2}\Big] \bigg\}.\nonumber \end{align} (B25) Recall that $$e_{j}^{n}=\eta _{j}^{n}+\xi _{j,h}^{n}$$. Using the triangle inequality on the error equation to split the error terms into terms of $$\eta _{j}^{n}$$ and $$\xi _{j,h}^{n}$$ gives   \begin{align*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} &+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \\ \leqslant&\ \frac{1}{2}\Vert\xi_{j,h}^{N}\Vert^{2} +\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2} \\ &+\frac{1}{2}\Vert\eta_{j}^{N}\Vert^{2} +\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla \eta_{j}^{n+1}\Vert ^{2} \end{align*}and   \begin{equation*} \begin{aligned} \frac{1}{2}\Vert\xi_{j,h}^{0}\Vert^{2} &+\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla \xi_{j,h} ^{0}\Vert^{2}\\ \leqslant&\ \frac{1}{2}\Vert e_{j}^{0}\Vert^{2} +\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla e_{j}^{0}\Vert^{2}\\ &+\frac{1}{2}\Vert\eta_{j}^{0}\Vert^{2} +\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla \eta_{j}^{0}\Vert^{2}\, . \end{aligned} \end{equation*}Applying inequality (B25), using the previous bounds for $$\eta _{j}^{n}$$ terms, and absorbing constants into a new constant C, completes the proof of Theorem 3.1. © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# An efficient algorithm for simulating ensembles of parameterized flow problems

, Volume Advance Article – May 24, 2018
26 pages

Loading next page...

/lp/ou_press/an-efficient-algorithm-for-simulating-ensembles-of-parameterized-flow-8HOsh8F03q
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/dry029
Publisher site
See Article on Publisher Site

### Abstract

Abstract Many applications of computational fluid dynamics require multiple simulations of a flow under different input conditions. In this paper, a numerical algorithm is developed to efficiently determine a set of such simulations in which the individually independent members of the set are subject to different viscosity coefficients, initial conditions and/or body forces. The proposed scheme, when applied to the flow ensemble, needs to solve a single linear system with multiple right-hand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and long-term stable under certain conditions on the time-step size and a parameter deviation ratio. A rigorous numerical error estimate shows the scheme is of first-order accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results. 1. Introduction Numerical simulations of incompressible viscous flows have important applications in engineering and science. In this paper we consider settings in which one wishes to obtain solutions for several different values of the physical parameters and several different choices for the forcing functions appearing in the partial differential equation (PDE) model. For example, in building low-dimensional surrogates for the PDE solution such as sparse-grid interpolants or proper orthogonal decomposition approximations, one has to first determine expensive approximation of solutions corresponding to several values of the parameters. Sensitivity analyses of solutions often need to determine approximate solutions for several parameter values and/or forcing functions. An important third example is quantifying the uncertainties of outputs from the model equations. Mathematical models should take into account the uncertainties invariably present in the specification of physical parameters and/or forcing functions appearing in the model equations. For flow problems, because the viscosity of the liquid or gas often depends on the temperature, an inaccurate measurement of the temperature would introduce some uncertainty into the viscosity of the flow. Direct measurements of the viscosity using flow meters and measurements of the state of the system are also prone to uncertainties. Of course, forcing functions, e.g., initial condition data, can and usually are also subject to uncertainty. In such cases, due to the lack of exact information, stochastic modeling is used to describe flows subject to a random viscosity coefficient and/or random forcing. Subsequently, numerical methods are employed to quantify the uncertainties in system output. It is known that uncertainty quantification, when a random sampling method such as a Monte Carlo method is used, could be computationally expensive for large-scale problems because each individual realization requires a large-scale computation but on the other hand, many realizations may be needed in order to obtain accurate statistical information about the outputs of interest. Therefore, for all the examples discussed and for many others, how to design efficient algorithms for performing multiple numerical simulations becomes a matter of great interest. The ensemble method that forms the basis for our approach was proposed in Jiang & Layton (2014); there, a set of J solutions of the Navier–Stokes equations (NSE) with distinct initial conditions and forcing terms is considered. All solutions are found, at each time step, by solving a linear system with one shared coefficient matrix and J right-hand sides (RHS), reducing both the storage requirements and computational costs of the solution process. The algorithm of Jiang & Layton (2014) is first-order accurate in time; it is extended to higher-order accurate schemes in Jiang (2015, 2017). Ensemble regularization methods are developed in the studies by Jiang (2015), Jiang & Layton (2015) and Takhirov et al. (2016) for high Reynolds number flows, and a turbulence model based on ensemble averaging is developed in Jiang et al. (2015). The ensemble algorithm has also been extended to magnetohydrodynamics flows in Mohebujjaman & Rebholz (2017), to natural convection problems in Fiordilino (2017a) and to parametrized flow problems in Gunzburger et al. (2017b). Ensemble algorithms incorporating reduced-order modeling techniques are studied in Gunzburger et al. (2017a, 2018). Recently, the ensemble method has been introduced in Fiordilino (2017b); Luo & Wang (2018a,b) for uncertainty quantification problems on random linear parabolic equations. In this paper, we develop a numerical scheme for simulating ensembles of the NSE flow problems in which not only the initial data and body force function, but also the viscosity coefficient, may vary from one ensemble member to another. Specifically, we consider a set of J NSE simulations on a bounded domain subject to no-slip boundary conditions in which the jth individual member solves the system   $$\left\{ \begin{array}{rcll} u_{j, t}+u_{j}\cdot\nabla u_{j}-\nabla \cdot (\nu_j \nabla u_{j})+\nabla p_{j} &=& f_{j}(x,t) \quad &\textrm{in}\ \varOmega\times [0, \infty), \\ \nabla\cdot u_{j} &=& 0 \quad &\textrm{in}\ \varOmega \times [0, \infty), \\ u_{j} &=& 0 \quad &\textrm{on}\ \partial\varOmega \times [0, \infty),\\ u_{j}(x,0) &=& u_{j}^{0}(x) \quad & \textrm{in} \ \varOmega, \end{array}\right.$$ (1.1)which, for each j, corresponds to a different variable kinematic viscosity $$\nu _j= \nu _j(x)$$ and/or distinct initial data $$u_{j}^{0}$$ and/or body forces $$f_{j}$$. In the sequel, it is assumed that $$\nu _j(x) \in L^\infty (\varOmega )$$ and $$\nu _j(x)\geqslant \nu _{j, \textrm{min}}>0$$. Due to the nonlinear convection term, implicit and semiimplicit schemes are invariably used for time integration. For a semiimplicit scheme the associated discrete linear systems would be different for each individually independent simulation, i.e., for each j. As a result, at each time step, J linear systems need to be solved to determine the ensemble, resulting in a huge computational effort. For a fully implicit scheme, the situation is even worse because one would have to solve many more linear systems due to the nonlinear solver iteration. To tackle this issue we propose a novel discretization scheme that results, at each time step, in a common coefficient matrix for all the ensemble members. 1.1 The ensemble-based semiimplicit scheme To focus on the main idea we temporarily ignore the spatial discretization and consider only the ensemble-based implicit–explicit temporal integration scheme   \left\{ \begin{aligned} \frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}+\overline{u}^{n}\cdot\nabla u_{j}^{n+1} +(u_{j}^{n}-\overline{u}^{n})\cdot\nabla u_{j}^{n} -\nabla \cdot (\overline{\nu} \nabla u_{j}^{n+1}) -\nabla \cdot \left((\nu_j-\overline{\nu})\nabla u_{j}^{n}\right) +\nabla p_{j}^{n+1} &=f_{j}^{n+1},\\ \nabla\cdot u_{j}^{n+1}&=0, \end{aligned}\right. (1.2)where $$\overline{u}^{n}$$ and $$\overline{\nu }$$ are the ensemble means of the velocity and viscosity coefficients, respectively, defined as   $$\overline{u}^{n}:=\frac{1}{J}\sum_{j=1}^{J}u_{j}^{n} \qquad \textrm{and}\qquad \overline{\nu}:=\frac{1}{J}\sum_{j=1}^{J}\nu_{j}.$$We also define $$\overline{\nu }_{\textrm{min}}:= \frac{1}{J}\sum _{j=1}^{J}\nu _{j, \textrm{min}}$$. After rearranging the system, we have, at time $$t_{n+1}$$,   \left\{ \begin{aligned} \frac{1}{\varDelta t}u_{j}^{n+1}\!+\overline{u}^{n}\cdot\nabla u_{j}^{n+1} \!-\! \nabla\!\cdot\! (\overline{\nu} \nabla u_{j}^{n+1}) \!+\!\nabla p_{j}^{n+1} &= f_{j}^{n+1}\!+\! \frac{1}{\varDelta t}u_{j}^{n} \!-(u_{j}^{n}\!-\overline{u}^{n})\cdot\nabla u_{j}^{n} \!+\!\nabla \!\cdot\! \left( (\nu_j\!-\!\overline{\nu}) \nabla u_{j}^{n} \right), \\ \nabla\cdot u_{j}^{n+1}&=0. \end{aligned}\right. (1.3)It is clear that the coefficient matrix of the resulting linear system will be independent of j. Thus, for the flow ensemble, to advance all members of the ensemble by one time step, we need to solve only a single linear system with J RHS. Compared with solving J individually independent simulations, this approach used with either a single LU factorization for small-scale problems or a block Krylov subspace method (Gutknecht, 2007; Parks et al., 2016) for large-scale problems is computationally more efficient and significantly reduces the required storage. When the size of the ensemble becomes huge, it can be subdivided into p sub-ensembles so as to balance memory, communication and computational costs and then (1.2) can be applied to each sub-ensemble. The rest of this section is devoted to establishing notation and to providing other preliminary information. Then in Section 2 we prove a conditional stability result for a fully discrete finite element discretization of (1.2). In Section 3 we derive an error estimate for the fully discrete approximation. Results of the preliminary numerical simulations that illustrate the theoretical results are given in Section 4, and Section 5 provides some concluding remarks. 1.2 Notation and preliminaries Let $$\varOmega$$ denote an open, regular domain in $$\mathbb{R}^{d}$$ for d = 2 or 3 having boundary denoted by $$\partial \varOmega$$. The $$L^{2}(\varOmega )$$ norm and inner product are denoted by ∥ ⋅ ∥ and (⋅, ⋅), respectively. The $$L^{p}(\varOmega )$$ norms and the Sobolev $$W^{k}_{p}(\varOmega )$$ norms are denoted by $$\|\cdot \|_{L^{p}}$$ and $$\|\cdot \|_{W_{p}^{k}}$$, respectively. The Sobolev space $$W_{2}^{k}(\varOmega )$$ is simply denoted by $$H^{k}(\varOmega )$$ and its norm by $$\|\cdot \|_{k}$$. For functions v(x, t) defined on (0, T), we define, for $$1\leqslant m <\infty$$,   $$\| v \|_{\infty,k} \textrm{ }:=\textrm{EssSup}_{[0,T]}\| v(\cdot, t)\|_{k}\qquad \textrm{and}\qquad \|v\|_{m,k} \textrm{ }:= \left( \int_{0}^{T}\|v(\cdot, t)\|_{k}^{m}\, \textrm{d}t \right) ^{1/m}.$$Given a time step $$\varDelta t$$, associated discrete norms are defined as   $${\left\vert \! \left\vert \! \left\vert v \right\vert\! \right\vert \! \right\vert}_{\infty,k}=\max\limits_{0\leqslant n\leqslant N}\Vert v^{n}\Vert_{k} \qquad \textrm{and}\qquad{\left\vert \! \left\vert \! \left\vert v \right\vert\! \right\vert \! \right\vert}_{m,k}:= \left(\sum_{n=0}^{N}\|v^{n}\|_{k}^{m}\varDelta t \right)^{1/m},$$where $$v^{n}=v(t_{n})$$ and $$t_{n}=n\varDelta t$$. Denote by $$H^{-1}(\varOmega )$$ the dual space of bounded linear functions on $$H^{1}_{0}(\varOmega )=\{v\in H^{1}\,:\, v=0 \ \textrm{on}\ \partial \varOmega \}$$; a norm on $$H^{-1}(\varOmega )$$ is given by   $$\|\,f\|_{-1}=\sup_{0\neq v\in H^{1}_{0}(\varOmega)}\frac{(\,f,v)}{\Vert\nabla v\Vert}.$$ The velocity space X and pressure space Q are given by   $$X:=[H_{0}^{1}(\varOmega)]^{d} \qquad \textrm{and}\qquad Q:=L_{0}^{2}(\varOmega)=\{q\in L^2(\varOmega)\,\,:\,\, \textstyle\int_\varOmega q\, \textrm{d}x =0\},$$respectively. The space of weakly divergence-free functions is   $$V\textrm{ }:=\{v\in X\,\,:\,\, (\nabla\cdot v,q)=0,\quad \forall\, q\in Q\} .$$ A weak formulation of (1.1) reads as follows: for j = 1, … , J, find $$u_j: \,[0,T]\rightarrow X$$ and $$p_j:\, [0,T]\rightarrow Q$$ for a.e. t ∈ (0, T] satisfying   \begin{equation*} \left\{ \begin{aligned} (u_{j,t},v)+(u_{j}\cdot\nabla u_{j},v)+(\nu_j\nabla u_{j},\nabla v) -(p_{j},\nabla\cdot v) &= (f_{j},v) &\forall\, v\in X, \\ (\nabla\cdot u_{j},q) &= 0 &\forall\, q\in Q, \end{aligned}\right. \end{equation*}with $$u_{j}(x,0)=u_{j}^{0}(x)$$. Our analysis is based on a finite element method (FEM) for spatial discretization. However, the results also extend, without much difficulty, to other variational discretization methods. Let $$X_{h}\subset X$$ and $$Q_{h}\subset Q$$ denote families of conforming velocity and pressure finite element spaces on a regular subdivision of $$\varOmega$$ into simplices; the family is parameterized by the maximum diameter h of any of the simplices. Assume that the pair of spaces $$(X_h,Q_h)$$ satisfy the discrete inf-sup (or $$\textrm{LBB}_h$$) condition required for the stability of the finite element approximation and that the finite element spaces satisfy the approximation properties   $$\inf_{v_h\in X_h}\| v- v_h \| \leqslant C h^{k+1}\Vert u \Vert_{k+1} \quad \quad\quad\quad \forall\, v\in [H^{k+1}(\varOmega)]^d,$$ (1.4)  $$\inf_{v_h\in X_h}\| \nabla ( v- v_h )\| \leqslant C h^k \Vert v\Vert_{k+1}\quad \;\;\,\quad\quad\quad \forall\, v\in [H^{k+1}(\varOmega)]^d,$$ (1.5)  $$\inf_{q_h\in Q_h}\| q- q_h \| \leqslant C h^{s+1}\Vert p\Vert_{s+1} \quad \quad\quad\quad\quad \forall\, q\in H^{s+1}(\varOmega),$$ (1.6)where the generic constant C > 0 is independent of mesh size h. An example for which the $$\textrm{LBB}_h$$ stability condition and the approximation properties are satisfied is the family of Taylor–Hood $$P^{s+1}$$–$$P^{s}$$, $$s\geqslant 1$$ element pairs. For details concerning FEMs see, e.g., Ciarlet (2002) and Girault & Raviart (1979, 1986); Gunzburger (1989); Layton (2008) for FEMs for the NSE. The discretely divergence-free subspace of $$X_{h}$$ is defined as   $$V_{h}\ :=\{v_{h}\in X_{h}\,\,:\,\,(\nabla\cdot v_{h},q_{h})=0\quad \forall\, q_{h}\in Q_{h}\}.$$Note that, in general, $$V_h\not \subset V$$. We assume the mesh and finite element spaces satisfy the standard inverse inequality   $$h\Vert\nabla v_{h}\Vert \leqslant C_{(\textrm{inv})}\Vert v_{h}\Vert \quad\forall\, v_{h}\in X_{h},$$ (1.7)which is known to hold for standard finite element spaces with locally quasi-uniform meshes (Brenner & Scott, 2008). We also define the standard explicitly skew-symmetric trilinear form   $$b^{\ast}(u,v,w):=\tfrac{1}{2}(u\cdot\nabla v,w)-\tfrac{1}{2}(u\cdot\nabla w,v),$$which satisfies the bounds (Layton, 2008)   $$b^{\ast}(u,v,w)\leqslant C \left(\Vert \nabla u\Vert\Vert u\Vert\right)^{1/2}\Vert\nabla v\Vert\Vert\nabla w \Vert \quad \forall\, u, v, w \in X,$$ (1.8)  $$b^{\ast}(u,v,w)\leqslant C \Vert \nabla u\Vert\Vert\nabla v\Vert\left(\Vert\nabla w \Vert\Vert w\Vert\right)^{1/2} \quad \forall\, u, v, w \in X .$$ (1.9)We also denote the exact and approximate solutions at $$t=t^n$$ as $$u_{j}^{n}$$ and $$u_{j, h}^{n}$$, respectively. 2. Stability analysis The fully discrete finite element discretization of (1.2) is given as follows. Given $$u_{j,h}^{0}\in X_{h}$$, for n = 0, 1, … , N − 1, find $$u_{j,h}^{n+1}\in X_{h}$$ and $$p_{j,h}^{n+1}\in Q_{h}$$ satisfying   $$\begin{cases} \left(\frac{u_{j,h}^{n+1}-u_{j,h}^{n}}{\varDelta t},v_{h}\right)+b^{\ast}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h})+b^{\ast}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n}, v_{h})-(p_{j,h}^{n+1},\nabla\cdot v_{h})\\ \qquad\qquad\quad\;\ \ \ \ + \, (\overline{\nu} \nabla u_{j,h}^{n+1},\nabla v_{h}) +\left( (\nu_j-\overline{\nu}) \nabla u_{j,h}^{n},\nabla v_{h} \right) =(f_{j}^{n+1},v_{h})\quad \forall\, v_{h}\in X_{h},\\ \big(\nabla\cdot u_{j,h}^{n+1},q_{h}\big)=0 \quad\forall\, q_{h}\in Q_{h}. \end{cases}$$ (2.1)We begin by proving the conditional, nonlinear, long-time stability of scheme (2.1) under a time-step condition and a parameter deviation condition. Theorem 2.1 (Stability). For all j = 1, … , J, if for some $$\mu$$, $$0\leqslant \mu <1$$, and some $$\varepsilon$$, $$0< \varepsilon \leqslant 2-2\sqrt{\mu }$$, the following time-step condition and parameter deviation condition both hold,   $$C\frac{\varDelta t}{ h \overline{\nu}_{\textrm{min}}}\left\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\right\Vert{}^{2} \leqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)},$$ (2.2)  $$\frac{\|\nu_j -\overline{\nu}\|_{\infty} }{\overline{\nu}_{\textrm{min}}} \leqslant \sqrt{\mu},\quad$$ (2.3)then scheme (2.1) is nonlinearly, long-time stable. In particular, for j = 1, … , J and for any $$N\geqslant 1$$, we have   $$\begin{split} &\quad\frac{1}{2}\|u_{j,h}^{N}\|^{2}+\frac{1}{4}\sum_{n=0}^{N-1}\|u_{j,h} ^{n+1}-u_{j,h}^{n}\|^{2}+\sum_{n=0}^{N-1}\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\\ &\qquad+\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_\textrm{{min}}}\right) \|\nabla u_{j,h}^{N}\|^{2}\\ &\qquad\qquad\leqslant\sum_{n=0}^{N-1}\frac{2\varDelta t}{\overline{\nu}_{\textrm{min}}}\|\,f_{j}^{n+1}\|_{-1}^{2}+ \frac{1} {2}\|u_{j,h}^{0}\|^{2} +\overline{\nu}_\textrm{{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\|\nabla u_{j,h}^{0}\|^{2}. \end{split}$$ (2.4) Proof. The proof is given in Appendix A. Remark 2.2 It is seen from (2.2) that the upper bound in the time-step condition increases as $$\varepsilon$$ decreases. As $$\varepsilon \rightarrow 0$$ the bound approaches $$1-\sqrt{\mu }$$. Because the relative deviation of the viscosity coefficient in (2.3) is bounded by $$\sqrt{\mu }$$, the two stability conditions are oppositional to each other. Remark 2.3 Noting that condition (2.2) depends only on known quantities such as the solution at $$t_n$$ and that scheme (2.1) is a one-step method, (2.2) can be used to adapt $$\varDelta t$$ in order to guarantee stability for the ensemble simulations. Remark 2.4 When the viscosity $$\nu _j$$ is a constant instead of being variable, the relative viscosity coefficient deviation ratio required for stability is still bounded by $$\sqrt{\mu }$$, and condition (2.3) becomes $$\max_{j}\frac{\left |\nu _j-\overline{\nu }\right |}{\overline{\nu }}\leqslant \sqrt{\mu }$$. 3. Error analysis In this section we give a detailed error analysis of the proposed method under the same type of time-step condition (with possibly a different constant C on the left-hand side of the inequality) and the same parameter deviation condition. Assuming that $$X_{h}$$ and $$Q_{h}$$ satisfy the $$\textrm{LBB}_{h}$$ condition, scheme (2.1) is equivalent to the following: given $$u_{j,h}^{0}\in V_{h}$$, for n = 0, 1, … , N − 1 find $$u_{j,h} ^{n+1}\in V_{h}$$ such that   \begin{align} \left(\frac{u_{j,h}^{n+1}-u_{j,h}^{n}}{\varDelta t},v_{h}\right)&+b^{\ast}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h})+b^{\ast}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h} ^{n},v_{h})\nonumber\\ &+\overline{\nu}(\nabla u_{j,h}^{n+1},\nabla v_{h})+((\nu_j-\overline{\nu})\nabla u_{j,h}^{n},\nabla v_{h})=(\,f_{j}^{n+1},v_{h})\quad\forall\, v_{h}\in V_{h}. \end{align} (3.1)To analyse the rate of convergence of the approximation we assume the following regularity for the exact solutions:   \begin{align*} u_{j} \in L^{\infty}(0,T;H^{k+1}(\varOmega))\cap H^{1}(0,T;H^{k+1}(\varOmega))\cap H^{2}(0,T;L^{2}(\varOmega)),\\ p_{j} \in L^{2}(0,T;H^{s+1}(\varOmega))\quad \textrm{and}\quad f_{j} \in L^{2} (0,T;L^{2}(\varOmega)). \end{align*}Let $$e_{j}^{n}=u_{j}^{n}-u_{j,h}^{n}$$ denote the approximation error of the jth simulation at the time instance $$t_n$$. We then have the following error estimates. Theorem 3.1 (Convergence of scheme (2.1)). For all j = 1, … , J, if for some $$\mu$$, $$0\leqslant \mu <1$$ and some $$\varepsilon$$, $$0< \varepsilon \leqslant 2-2\sqrt{\mu }$$ the following time-step condition and parameter deviation condition both hold,   $$C\frac{\varDelta t}{\overline{\nu}_{\textrm{min}} h}\left\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\right\Vert{}^{2} \leqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)},$$ (3.2)  $$\frac{\|\nu_j -\overline{\nu}\|_{\infty} }{\overline{\nu}_{\textrm{min}}} \leqslant \sqrt{\mu},\quad$$ (3.3)then there exists a positive constant C independent of the time step such that   \begin{align*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} &+\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left( 1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \\ \leqslant&\, e^{\frac{CT}{\overline{\nu}_{\textrm{min}}^{3}}} \left\{\frac{1}{2}\Vert e_{j}^{0}\Vert^{2} +\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla e_{j}^{0}\Vert^{2}\right. \\ &\qquad\qquad +\frac{1}{2}h^{2k+2}\Vert u_j^0 \Vert_{ k+1}^2+\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\overline{\nu}_{\textrm{min}}\varDelta t h^{2k} \Vert u_{j}^{0}\Vert_{k+1}^{2}\\ &\qquad\qquad+C\varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla u_{j,t}\Vert^2_{2,0}+C\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert \! \left\vert u_j \right\vert\!\right\vert \! \right\vert}^2_{2,k+1} +C\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \\ &\qquad\qquad+Ch^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t \Vert \nabla u_{j,t}\Vert^2_{2,0} + C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \\ &\qquad\qquad+C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2\Vert \nabla u_{j,t}\Vert^2_{2,0}+ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \\ &\qquad\qquad+C\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2}+C\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}\Vert u_{j,t}\Vert_{2,k+1}^{2} +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}\Vert u_{j,tt}\Vert_{2,0}^{2}\Bigg\} \\ &+\frac{1}{2}h^{2k+2}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}_{\infty, k+1}^2+\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left( 1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}. \end{align*} Proof. The proof is given in Appendix B. In particular, when Taylor–Hood elements (k = 2, s = 1) are used, i.e., the $$C^{0}$$ piecewise-quadratic velocity space $$X_{h}$$ and the $$C^{0}$$ piecewise-linear pressure space $$Q_{h}$$, we have the following estimate. Corollary 3.2 Assume that $$\Vert e_{j}^{0}\Vert$$ and $$\Vert \nabla e_j^0\Vert$$ are both $$\mathcal{O}$$(h) accurate or better. Then if $$(X_{h},Q_{h})$$ is chosen as the $$(P_2, P_1)$$ Taylor–Hood element pair, we have   \begin{equation*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} +\frac{1}{15} \frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\Big( 1-\frac{\sqrt{\mu}}{2}\Big)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \leqslant C (h^{2} + \varDelta t^2 +h \varDelta t). \end{equation*} 4. Numerical experiments In this section we illustrate the effectiveness of our proposed method (1.2) and the associated theoretical analyses in Sections 2 and 3 by considering two examples: a Green–Taylor vortex problem and a flow between two offset cylinders. The first problem has a known exact solution that is used to illustrate the error analysis. The second example does not have an analytic solution but has complex flow structures; it is used to check the stability analysis and demonstrate the necessity and efficiency of the proposed method. 4.1 The Green–Taylor vortex problem The Green–Taylor vortex flow is commonly used for testing convergence rates, e.g., see Chorin (1968), Tafti (1996), John & Layton (2002), Berselli (2005), Barbato et al. (2007) and Jiang & Tran (2015). The Green–Taylor vortex solution given by   \begin{align} & u(x,y,t)=-\cos(\omega\pi x)\sin(\omega\pi y)e^{-2\omega^{2}\pi^{2}t/\tau }, \nonumber\\ & v(x,y,t)=\sin(\omega\pi x)\cos(\omega\pi y)e^{-2\omega^{2}\pi^{2}t/\tau }, \\ & p(x,y,t)=-\tfrac{1}{4}(\cos(2\omega\pi x)+\cos(2\omega\pi y))e^{-4\omega ^{2}\pi^{2}t/\tau} \nonumber \end{align} (4.1)satisfies the NSE in $$\varOmega =(0,1)^{2}$$ for $$\tau =\textrm{Re}$$ and initial condition   $$u^{0}=\big(\!\!-\cos(\omega\pi x)\sin(\omega\pi y), \sin(\omega\pi x)\cos(\omega\pi y)\big)^{\textrm{T}}.$$The solution consists of an $$\omega \times \omega$$ array of oppositely signed vortices that decay as $$t\rightarrow \infty$$. In the following numerical tests we take $$\omega =1$$, $$\nu =1/\textrm{Re}$$, T = 1, h = 1/m and $$\varDelta t/h=2/5$$. The boundary condition is assumed to be inhomogeneous Dirichlet, that is, the boundary values match that of the exact solution. We consider an ensemble of two members, $$u_{1}$$ and $$u_{2}$$, corresponding to two incompressible NSE simulations with different viscosity coefficients $$\nu _j$$ and initial conditions $$u_{j, 0}$$. We investigate the ensemble simulations and compare them with independent simulations. For j = 1, 2 we define by $$\mathcal{E}^{\, E}_ {\, j} = u_{j}-u_{j, h}$$ the approximation error of the jth member of the ensemble simulation and by $$\mathcal{E}^{\, S}_{\, j} = u_{j}-u_{j, h}$$ the approximation error of the jth independently determined simulation. Here, the superscript ‘E’ stands for ‘ensemble’ whereas ‘S’ stands for ‘independent’. Case 1. We set the viscosity coefficient $$\nu _1=0.2$$ and initial condition $$u_{1, 0}= (1+\varepsilon )u^0$$ for the first member $$u_1$$ and $$\nu _2=0.3$$ and $$u_{2, 0}= (1-\varepsilon )u^0$$ for the second member $$u_2$$, where $$\varepsilon = 10^{-3}$$. For this choice of parameters, we have $$\vert \nu _j -\overline{\nu }\vert /\overline{\nu } =\frac{1}{5}$$ for both j = 1 and j = 2 so that condition (2.3) is satisfied. We first apply the ensemble algorithm; results are shown in Table 1. It is seen that the convergence rate for $$u_{1}$$ and $$u_{2}$$ is first order. Table 1 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors for ensemble simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  View Large Table 1 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors for ensemble simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$1.05\cdot 10^{-2}$$  –  $$4.17\cdot 10^{-2}$$  –  $$7.36\cdot 10^{-3}$$  –  $$2.53\cdot 10^{-2}$$  –  40  $$5.86\cdot 10^{-3}$$  0.85  $$2.21\cdot 10^{-2}$$  0.91  $$3.87\cdot 10^{-3}$$  0.93  $$1.31\cdot 10^{-2}$$  0.95  80  $$3.10\cdot 10^{-3}$$  0.92  $$1.14\cdot 10^{-2}$$  0.95  $$2.02\cdot 10^{-3}$$  0.94  $$6.70\cdot 10^{-3}$$  0.97  160  $$1.59\cdot 10^{-3}$$  0.96  $$5.81\cdot 10^{-3}$$  0.97  $$1.03\cdot 10^{-3}$$  0.97  $$3.39\cdot 10^{-3}$$  0.98  View Large We next compare the ensemble simulations with independent simulations. To this end, we perform each NSE simulation independently using the same discretization setup. The associated approximation errors are listed in Table 2. Comparing with Table 1, we observe that the ensemble simulation is able to achieve accuracies close to those of the independent, more costly simulations. Table 2 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors in independent simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  View Large Table 2 For the Green–Taylor vortex problem (Case 1) and for a sequence of uniform grid sizes h, errors in independent simulations of two members with inputs $$\nu _1=0.2$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.3$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$1.01\cdot 10^{-2}$$  –  $$3.88\cdot 10^{-2}$$  –  $$7.88\cdot 10^{-3}$$  –  $$2.76\cdot 10^{-2}$$  –  40  $$5.47\cdot 10^{-3}$$  0.89  $$2.04\cdot 10^{-2}$$  0.93  $$4.24\cdot 10^{-3}$$  0.90  $$1.44\cdot 10^{-2}$$  0.93  80  $$2.85\cdot 10^{-3}$$  0.94  $$1.05\cdot 10^{-2}$$  0.96  $$2.22\cdot 10^{-3}$$  0.93  $$7.41\cdot 10^{-3}$$  0.96  160  $$1.46\cdot 10^{-3}$$  0.97  $$5.30\cdot 10^{-3}$$  0.98  $$1.13\cdot 10^{-3}$$  0.97  $$3.76\cdot 10^{-3}$$  0.98  View Large Case 2. We now set $$\nu _1= 0.01$$ and $$\nu _2= 0.49$$ while keeping the same initial conditions as for Case 1. With this choice of parameters, $$\vert \nu _j -\overline{\nu }\vert /\overline{\nu } =\frac{24}{25}$$ for both j = 1 and j = 2, which still satisfies (2.3) but is closer to the upper limit. The ensemble simulation errors are listed in Table 3, which shows that the rate of convergence for the second member is nearly 1 and for the first member is approaching 1. Table 3 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in ensemble simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  View Large Table 3 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in ensemble simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ E}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ E}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ E}_{\ 2}\|_{2,0}$$  Rate  20  $$2.91\cdot 10^{-2}$$  –  $$2.96\cdot 10^{-1}$$  –  $$3.50\cdot 10^{-3}$$  –  $$9.94\cdot 10^{-3}$$  –  40  $$1.86\cdot 10^{-2}$$  0.65  $$1.80\cdot 10^{-1}$$  0.71  $$1.65\cdot 10^{-3}$$  1.08  $$4.97\cdot 10^{-3}$$  1  80  $$1.08\cdot 10^{-2}$$  0.78  $$1.02\cdot 10^{-1}$$  0.83  $$8.53\cdot 10^{-4}$$  0.95  $$2.52\cdot 10^{-3}$$  0.98  160  $$5.89\cdot 10^{-3}$$  0.87  $$5.46\cdot 10^{-2}$$  0.90  $$4.32\cdot 10^{-4}$$  0.98  $$1.27\cdot 10^{-3}$$  0.98  View Large The approximation errors for two independent simulations under the same discretization setup are listed in Table 4. Comparing the ensemble simulation results in Table 3 with the independent simulations, we find that the accuracy of the first member in the ensemble simulation degrades slightly whereas that of the second member in the ensemble simulation improves a bit. Overall, the ensemble simulation is able to achieve the same order of accuracy as the independent simulations. Table 4 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in independent simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  View Large Table 4 For the Green–Taylor vortex problem (Case 2) and for a sequence of uniform grid sizes h, errors in independent simulations of two members: $$\nu _1=0.01$$, $$u_{1, 0}= (1+10^{-3})u^0$$ and $$\nu _2=0.49$$, $$u_{2, 0}= (1-10^{-3})u^0$$ 1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  1/h  $$\|\mathcal{E}^{\ S}_{\ 1}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 1}\|_{2,0}$$  Rate  $$\|\mathcal{E}^{\ S}_{\ 2}\|_{\infty , 0}$$  Rate  $$\|\nabla \mathcal{E}^{\ S}_{\ 2}\|_{2,0}$$  Rate  20  $$3.19\cdot 10^{-2}$$  –  $$2.95\cdot 10^{-1}$$  –  $$5.49\cdot 10^{-3}$$  –  $$1.79\cdot 10^{-2}$$  –  40  $$1.67\cdot 10^{-2}$$  0.93  $$1.54\cdot 10^{-1}$$  0.93  $$3.03\cdot 10^{-3}$$  0.86  $$9.38\cdot 10^{-3}$$  0.94  80  $$8.56\cdot 10^{-3}$$  0.97  $$7.90\cdot 10^{-2}$$  0.97  $$1.59\cdot 10^{-3}$$  0.93  $$4.81\cdot 10^{-3}$$  0.96  160  $$4.33\cdot 10^{-3}$$  0.98  $$3.99\cdot 10^{-2}$$  0.98  $$8.18\cdot 10^{-4}$$  0.96  $$2.44\cdot 10^{-3}$$  0.98  View Large 4.2 Flow between two offset cylinders Next we check the stability of our algorithm by considering the problem of a flow between two offset circles (see, e.g., Jiang & Layton, 2014, 2015; Jiang, 2015; Jiang et al., 2015). The domain is a disk with a smaller off-center obstacle inside. Letting $$r_{1}=1$$, $$r_{2}=0.1$$ and $$c=(c_{1},c_{2})=\big(\frac{1}{2} ,0\big)$$, the domain is given by   $$\varOmega=\big\{(x,y)\,\,:\,\,x^{2}+y^{2}\leqslant r_{1}^{2} \,\,\textrm{and}\,\, (x-c_{1})^{2} +(y-c_{2})^{2}\geqslant r_{2}^{2}\big\}.$$ The flow is driven by a counterclockwise rotational body force   $$f(x,y,t)=\big(-6y(1-x^{2}-y^{2}), 6x(1-x^{2}-y^{2})\big)^{\textrm{T}}$$with no-slip boundary conditions imposed on both circles. The flow between the two circles shows interesting structures interacting with the inner circle. A von K$$\acute{\textrm{a}}$$rm$$\acute{\textrm{a}}$$n vortex street is formed behind the inner circle and then re-interacts with that circle and with itself, generating complex flow patterns. We consider multiple numerical simulations of the flow with different viscosity coefficients using the ensemble-based algorithm (2.1). For spatial discretization, we apply the Taylor–Hood element pair on a triangular mesh that is generated by Delaunay triangulation with 80 mesh points on the outer circle and 60 mesh points on the inner circle and with refinement near the inner circle, resulting in 18,638 degrees of freedom; see Fig. 1. In order to illustrate the stability analysis, we select two different sets of viscosity coefficients for Case 1:$$\nu _1=0.005$$, $$\nu _2=0.039$$, $$\nu _3=0.016$$; Case 2:$$\nu _1=0.005$$, $$\nu _2=0.041$$, $$\nu _3=0.014$$. Note that the viscosity coefficients in Case 2 are obtained by slightly perturbing those in Case 1. The average of the viscosity coefficients is $$\overline{\nu }= 0.02$$ for both cases. However, the stability condition (2.3) is satisfied in the first case but is not satisfied in the second one, i.e., we have Case 1:$$\frac{\vert \nu _1-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{4}<1$$, $$\frac{\vert \nu _2-\overline{\nu }\vert }{\overline{\nu }}=\frac{19}{20}<1$$, $$\frac{\vert \nu _3-\overline{\nu }\vert }{\overline{\nu }}=\frac{1}{5}<1$$; Case 2:$$\frac{\vert \nu _1-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{4}<1$$, $$\frac{\vert \nu _2-\overline{\nu }\vert }{\overline{\nu }}=\frac{21}{20}>1$$, $$\frac{\vert \nu _3-\overline{\nu }\vert }{\overline{\nu }}=\frac{3}{10}<1$$. The second member of Case 2 has a perturbation ratio greater than 1. Simulations of both cases are subject to the same initial condition and body forces for all ensemble members. In particular, the initial condition is generated by solving the steady Stokes problem with viscosity $$\nu =0.02$$ and the same body force f(x, y, t). All the simulations are run over the time interval [0, 5] with a time-step size $$\varDelta t= 0.01$$. This time-step size is smaller than the one ensuring a stable ensemble simulation in Case 1, thus condition (2.2) holds. Furthermore, the same $$\varDelta t$$ is used in Case 2. Because the viscosity coefficients in Case 2 are slightly perturbed from those in Case 1, and $$\varDelta t$$ is chosen small, we believe condition (2.2) still holds. But condition (2.3) no longer holds. Therefore, we expect the ensemble simulation to be unstable even when using the small time-step size $$\varDelta t=0.01$$. For testing the stability, we use the kinetic energy as a criterion and compare the ensemble simulation results with independent simulations using the same mesh and time-step size. Fig. 1. View largeDownload slide Mesh for the flow between two offset cylinders example. Fig. 1. View largeDownload slide Mesh for the flow between two offset cylinders example. Fig. 2. View largeDownload slide For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 2. View largeDownload slide For the flow between two offset cylinders, Case 1, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 3. View largeDownload slide For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). Fig. 3. View largeDownload slide For the flow between two offset cylinders, Case 2, the energy evolution of the ensemble (Ens.) and independent simulations (Ind.). The comparison of the energy evolution of ensemble-based simulations with the corresponding independent simulations is shown in Figs 2 and 3. It is seen that for Case 1 the ensemble simulation is stable, but for Case 2 it becomes unstable. This phenomenon coincides with our stability analysis since condition (2.3) holds for all members of Case 1, but does not hold for the second member of Case 2. Indeed, it is observed from Fig. 3 that the energy of the second member in Case 2 blows up after t = 3.7, then affecting the other two members, and this results in their energy dramatically increasing after t = 4.7. Next we use this test example to illustrate the necessity of ensemble simulations. Indeed, the ensemble calculation is much needed when flow problems are under parameter uncertainty. This is because even though one could obtain a ‘best estimate’ of the parameter that is close enough to the true parameter value, the corresponding individual solution may vary significantly from the true solution due to the nonlinearity of the problems. However, the ensemble mean of the model problems at several selected parameter samples, which are drawn from a probability distribution (usually experimentally determined), tends to smooth out possible errors and gives a better forecast than the single run using the ‘best estimate’ parameter value. Here we consider a simple computational setting in which the true viscosity $$\nu = 0.01$$. Due to the errors in measurements, the ‘best estimate’ of the viscosity is 0.00995. Suppose the viscosity value empirically follows a uniform distribution on [0.009, 0.011]. We first run the simulation with viscosity $$\nu =0.01$$ on the domain partitioned by a Delaunay triangular mesh with 80 mesh points on the outer circle and 40 mesh points on the inner circle, take the uniform time-step size to be $$\varDelta t=0.002$$ and label the output solution ‘true’ as it is the numerical solution associated with the true datum. We then run the simulation associated to the ‘best estimate’ $$\nu =0.00995$$ and label its solution ‘best estimate’, which is to be compared with the ensemble forecast. For the latter, we draw J (= 16, 32, 64) independent, uniformly distributed samples of $$\nu$$ from [0.009, 0.011], and seek the mean of these J simulation results by using the proposed ensemble algorithm. Time evolution of the kinetic energy for the single ‘best estimate’ ensemble means of J simulations are shown in Fig. 4, together with that of the ‘true’ solution. It is observed that although the ‘best estimate’ of $$\nu$$ is very close to the true value, the difference between these two solutions is relatively large. However, the ensemble mean is able to yield a more accurate forecast. The greater the number of realizations, the better forecast that ensemble mean could achieve. Fig. 4. View largeDownload slide For the flow between two offset cylinders, comparison of ‘true’ solution, ‘best estimate’ solution and ensemble forecast. Fig. 4. View largeDownload slide For the flow between two offset cylinders, comparison of ‘true’ solution, ‘best estimate’ solution and ensemble forecast. Finally, we illustrate the efficiency of the proposed ensemble algorithm. To this end, we solve a group of J flow problems using a sequence of individual simulations and one ensemble simulation, respectively. These problems are subject to fixed initial and boundary conditions and a forcing function, but different viscosity parameter values. We then compare the efficiency of these two approaches by considering the CPU time. Assume the viscosity parameter $$\nu _j$$ to be a random variable that distributes uniformly on the interval [0.4, 0.5]. We generate a set of J (= 2, 4, 8, 16) random samples using MATLAB, partition the domain by a Delaunay triangular mesh with 80 mesh points on the outer circle and 40 mesh points on the inner circle and take the uniform time-step size to be $$\varDelta t=0.002$$. Both individual and ensemble simulations are implemented by using FreeFem++ with UMFPACK, which solves the linear system with multifrontal LU factorization. Numerical results together with the CPU time for simulations (measured in seconds) are listed in Table 5, where $$\Vert \overline{u}^{E}(T)\Vert$$ is the $$L_2$$ norm of the mean of the velocity obtained by the ensemble algorithm (2.1) at the final simulation time T = 0.1 and $$\Vert \overline{u}^{S}(T)\Vert$$ is that of individual simulations. Table 5 Efficiency comparison between the ensemble simulations and individual simulations on a set of J parametrized NSE flow between two offset cylinders problems at T = 0.1. It is observed that the ensemble algorithm outperforms the individual simulations as the former spends less CPU times than the latter but keeps the same accuracy J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  View Large Table 5 Efficiency comparison between the ensemble simulations and individual simulations on a set of J parametrized NSE flow between two offset cylinders problems at T = 0.1. It is observed that the ensemble algorithm outperforms the individual simulations as the former spends less CPU times than the latter but keeps the same accuracy J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  J  $$\Vert \bar{u}^{E}(T)\Vert$$  CPU (s)  $$\Vert \bar{u}^{S}(T)\Vert$$  CPU (s)  2  1.9870  44.2404  1.9870  101.244  4  1.8674  60.9567  1.8674  199.142  8  1.8893  95.1164  1.8893  397.261  16  1.9227  164.637  1.9227  788.407  View Large Define the speedup factor of the ensemble algorithm to be the ratio of execution time for the sequence of individual simulations to that of the ensemble simulation. It is observed that for this simple test, the computational time for completing the sequence of individual simulations increases almost linearly as the size of the ensemble increases, which is expected as all the work including assembly of matrices and LU decomposition needs to be done repeatedly for each ensemble member. However, the ensemble algorithm takes significantly less CPU time while maintaining the same accuracy. It is seen that for J = 2 the speedup factor is 2.30, while for J = 16 the speedup factor increases to 4.81. This saving rate will keep increasing as the size of the ensemble increases. 5. Conclusions In this paper, we consider a set of Navier–Stokes simulations in which each member may be subject to a distinct viscosity coefficient, initial conditions and/or body forces. An ensemble algorithm is developed for the group by which all the flow ensemble members, after discretization, share a common linear system with different right-hand-side vectors. This leads to a great saving in both storage requirements and computational costs. The stability and accuracy of the ensemble method are analysed. Two numerical experiments are presented. The first is for Green–Taylor flow and serves to illustrate the first-order accuracy in time of the ensemble-based scheme. The second is for a flow between two offset cylinders and serves to show that our stability analysis is sharp and the proposed method is effective. As a next step, we will investigate higher-order accurate schemes for the flow ensemble simulations. Funding US Air Force Office of Scientific Research (FA9550-15-1-0001); US Department of Energy Office of Science (DE-SC0009324, DE-SC0016591, DE-SC0016540); US National Science Foundation (DMS-1522672, DMS-1720001); University of Missouri Research Board. References Barbato, D., Berselli, L. & Grisanti, C. ( 2007) Analytical and numerical results for the rational large eddy simulation model. J. Math. Fluid Mech. , 9, 44– 74. Google Scholar CrossRef Search ADS   Berselli, L. ( 2005) On the large eddy simulation of the Taylor-Green vortex. J. Math. Fluid Mech. , 7, S164– S191. Google Scholar CrossRef Search ADS   Brenner, S. & Scott, R. ( 2008) The Mathematical Theory of Finite Element Methods , 3rd edn. New York: Springer-Verlag Google Scholar CrossRef Search ADS   Chorin, A. ( 1968) Numerical solution for the Navier–Stokes equations. Math. Comp. , 22, 745– 762. Google Scholar CrossRef Search ADS   Ciarlet, P. ( 2002) The Finite Element Method for Elliptic Problems . Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Fiordilino, J. ( 2017a) A second order ensemble timestepping algorithm for natural convection. Preprint, arXiv:1708.00488. Fiordilino, J. ( 2017b) Ensemble timestepping algorithms for the heat equation with uncertain conductivity. Preprint, arXiv:1708.00893. Girault, V. & Raviart, P.-A. ( 1979) Finite Element Approximation of the Navier-Stokes Equations . Lecture Notes in Mathematics , vol. 749. Berlin: Springer. Google Scholar CrossRef Search ADS   Girault, V. & Raviart, P.-A. ( 1986) Finite Element Methods for Navier-Stokes Equations - Theory and Algorithms . Berlin: Springer. Google Scholar CrossRef Search ADS   Gunzburger, M. ( 1989) Finite Element Methods for Viscous Incompressible Flows - A Guide to Theory, Practice, and Algorithms.  London: Academic Press. Gunzburger, M., Jiang, N. & Schneier, M. ( 2017a) An ensemble-proper orthogonal decomposition method for the nonstationary Navier-Stokes equations. SIAM J. Numer. Anal. , 55, 286– 304. Google Scholar CrossRef Search ADS   Gunzburger, M., Jiang, N. & Schneier, M. ( 2018) A higher-order ensemble/proper orthogonal decomposition method for the nonstationary Navier-Stokes equations. Int. J. Numer. Anal. Model. , 15, 608– 627. Gunzburger, M., Jiang, N. & Wang, Z. ( 2017b) A second-order time-stepping scheme for simulating ensembles of parameterized flow problems. Comput. Meth. Appl. Math.  (in press). https://doi.org/10.1515/cmam-2017-0051. Gutknecht, M. ( 2007) Block Krylov space methods for linear systems with multiple right-hand sides: an introduction. Modern Mathematical Models, Methods and Algorithms for Real World Systems (A. H. Siddiqi, I. S. Duff and O. Christensen eds). New Delhi: Anamaya Publishers, pp. 420– 447. Jiang, N. ( 2015) A higher order ensemble simulation algorithm for fluid flows. J. Sci. Comput. , 64, 264– 288. Google Scholar CrossRef Search ADS   Jiang, N. ( 2017) A second-order ensemble method based on a blended backward differentiation formula timestepping scheme for time-dependent Navier-Stokes equations. Numer. Meth. Partial. Diff. Eqs. , 33, 34– 61. Google Scholar CrossRef Search ADS   Jiang, N. & Layton, W. ( 2014) An algorithm for fast calculation of flow ensembles. Int. J. Uncertain. Quantif. , 4, 273– 301. Google Scholar CrossRef Search ADS   Jiang, N. & Layton, W. ( 2015) Numerical analysis of two ensemble eddy viscosity numerical regularizations of fluid motion. Numer. Meth. Partial Diff. Eqs. , 31, 630– 651. Google Scholar CrossRef Search ADS   Jiang, N. & Tran, H. ( 2015) Analysis of a stabilized CNLF method with fast slow wave splittings for flow problems. Comput. Meth. Appl. Math. , 15, 307– 330. Jiang, N., Kaya, S. & Layton, W. ( 2015) Analysis of model variance for ensemble based turbulence modeling. Comput. Meth. Appl. Math. , 15, 173– 188. John, V. & Layton, W. ( 2002) Analysis of numerical errors in large eddy simulation. SIAM J. Numer. Anal. , 40, 995– 1020. Google Scholar CrossRef Search ADS   Layton, W. ( 2008) Introduction to the Numerical Analysis of Incompressible Viscous Flows . Philadelphia: Society for Industrial and Applied Mathematics (SIAM) . Google Scholar CrossRef Search ADS   Luo, Y. & Wang, Z. ( 2018a) An ensemble algorithm for numerical solutions to deterministic and random parabolic PDEs. SIAM J. Numer. Anal. , 56, 859– 876. Google Scholar CrossRef Search ADS   Luo, Y. & Wang, Z. ( 2018b) A multilevel Monte Carlo ensemble scheme for solving random parabolic PDEs. Preprint, arXiv:1802.05743. Mohebujjaman, M. & Rebholz, L. ( 2017) An efficient algorithm for computation of MHD flow ensembles. Comput. Meth. Appl. Math. , 17, 121– 138. Google Scholar CrossRef Search ADS   Parks, M., Soodhalter, K. M. & Szyld, D. B. ( 2016) A block recycled GMRES method with investigations into aspects of solver performance. Preprint, arXiv:1604.01713. Tafti, D. ( 1996) Comparison of some upwind-biased high-order formulations with a second order central-difference scheme for time integration of the incompressible Navier-Stokes equations. Comput. & Fluids , 25, 647– 665. Google Scholar CrossRef Search ADS   Takhirov, A., Neda, M. & Waters, J. ( 2016) Time relaxation algorithm for flow ensembles. Numer. Meth. Partial Diff. Eqs. , 32, 757– 777. Google Scholar CrossRef Search ADS   Thomée, V. ( 1997) Galerkin Finite Element Methods for Parabolic Problems . Berlin, Heidelberg: Springer. Google Scholar CrossRef Search ADS   Appendix A. Proof of Theorem 2.1 Setting $$v_{h}=u_{j,h}^{n+1}$$ and $$q_h = p_{j, h}^{n+1}$$ in (2.1) and then adding two equations we obtain   \begin{align*} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert ^{2}+\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2}+ \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n}, u_{j,h}^{n},u_{j,h}^{n+1}) \\ &+ \varDelta t \Vert \overline{\nu}^{\frac{1}{2}} \nabla u_{j,h}^{n+1}\Vert^{2}= \varDelta t (f_{j}^{n+1}, u_{j,h}^{n+1}) - \varDelta t \left( \left( \nu_j - \overline{\nu}\right)\nabla u_{j,h}^n, \nabla u_{j,h}^{n+1}\right). \end{align*}Note that $$\overline{\nu }\geqslant \overline{\nu }_{\textrm{min}}>0$$. We apply Young’s inequality to the terms on the RHS and have, for all $$\alpha , \beta>0$$,   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} + \overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla u_{j,h}^{n+1}\Vert^{2} + \varDelta t b^{*} \left(u_{j,h}^{n} -\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1} -u_{j,h}^{n} \right)\nonumber\\ \leqslant &\frac{\alpha\overline{\nu}_{\textrm{min}}\varDelta t}{8} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} + \frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2} +\frac{\beta\overline{\nu}_{\textrm{min}}\varDelta t}{4} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} +\frac{\|\nu_j-\overline{\nu}\|^2_{\infty}\varDelta t}{\beta\overline{\nu}_{\textrm{min}}}\Vert \nabla u_{j,h}^n\Vert^2. \end{align} (A1)Both $$\frac{\beta \overline{\nu }_{\textrm{min}}\varDelta t}{4} \Vert \nabla u_{j,h}^{n+1} \Vert ^{2}$$ and $$\frac{\|\nu _j-\overline{\nu }\|_{\infty }^2\varDelta t}{\beta \overline{\nu }_{\textrm{min}}}\Vert \nabla u_{j,h}^n\Vert ^2$$ on the RHS of (A1) need to be absorbed into $$\overline{\nu }_{\textrm{min}}\varDelta t \Vert \nabla u_{j,h}^{n+1}\Vert ^{2}$$ on the LHS. To this end, we minimize $$\frac{\beta \overline{\nu }_{\textrm{min}}\varDelta t}{4}+\frac{\|\nu _j-\overline{\nu }\|_{\infty }^2\varDelta t}{\beta \overline{\nu }_{\textrm{min}}}$$ by selecting $$\beta = \frac{2 \| \nu _j-\overline{\nu }\|_{\infty } }{\overline{\nu }_{\textrm{min}}}$$ so that (A1) becomes   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{2}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} + \overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla u_{j,h}^{n+1}\Vert^{2} + \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})\nonumber\\ \leqslant&\, \frac{\alpha\overline{\nu}_{\textrm{min}}\varDelta t}{8} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} + \frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert f_{j}^{n+1}\Vert_{-1}^{2} +\frac{\| \nu_j-\overline{\nu}\|_{\infty}\varDelta t}{2} \Vert\nabla u_{j,h}^{n+1} \Vert^{2} +\frac{\| \nu_j-\overline{\nu}\|_{\infty} \varDelta t}{2}\Vert \nabla u_{j,h}^n\Vert^2. \end{align} (A2)Next we bound the trilinear term using the inequality (1.9) and the inverse inequality (1.7), obtaining   \begin{align*} - \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})&\leqslant C\varDelta t\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert\Vert\nabla u_{j,h} ^{n}\Vert\left[\Vert \nabla ( u_{j,h}^{n+1}-u_{j,h}^{n})\Vert\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert\right]^{1/2}\\ &\leqslant C\varDelta t\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert\Vert\nabla u_{j,h} ^{n}\Vert Ch^{-\frac{1}{2}} \Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert. \end{align*}Using Young’s inequality again gives   $$- \varDelta t b^{*}(u_{j,h}^{n}-\overline{u}_{h}^{n},u_{j,h}^{n},u_{j,h}^{n+1}-u_{j,h} ^{n})\leqslant C \frac{\varDelta t^{2}}{h} \Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2} \Vert\nabla u_{j,h}^{n}\Vert^{2}+ \frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h} ^{n}\Vert^{2}.$$ (A3)Substituting (A3) into (A2) and combining like terms, we have   \begin{align*} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\\&+\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2} \leqslant\frac{2\varDelta t}{\alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2} + C \frac{\varDelta t^{2}}{h} \Vert\nabla( u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2} \Vert\nabla u_{j,h}^{n}\Vert^{2} \\&+\frac{\| \nu_j-\overline{\nu}\|_{\infty} \varDelta t}{2}\Vert \nabla u_{j,h}^n\Vert^2. \end{align*}For any $$0<\sigma <1$$,   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+ \overline{\nu}_{\textrm{min}}\varDelta t \Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)\left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert^{2}\right)\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \left[(1-\sigma)\Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla u_{j,h}^{n}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}} \varDelta t \left[\sigma\Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla u_{j,h}^n\Vert^2 \leqslant\frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert f_{j}^{n+1}\Vert_{-1}^{2}\text{ .} \end{align} (A4)Select $$\alpha =4-\frac{2(\sigma +1)}{\sigma }\sqrt{\mu }$$. Since $$\alpha$$ is supposed to be greater than 0, we have   $$\sigma>\frac{\sqrt{\mu}}{2-\sqrt{\mu}} \in (0,1).$$ (A5)Now taking $$\sigma = \frac{\sqrt{\mu }+\varepsilon }{2-\sqrt{\mu }}$$, where $$\varepsilon \in (0, 2-2\sqrt{\mu })$$, (A4) becomes   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\alpha}{8}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+ \overline{\nu}_{\textrm{min}}\varDelta t \Big(1-\frac{\alpha}{4}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big)\left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \bigg[\Big( (1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu}-(1-\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\Big) -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\bigg] \Vert\nabla u_{j,h}^{n}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}} \varDelta t \Big[ \frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\Big]\Vert \nabla u_{j,h}^n\Vert^2\leqslant\frac{2\varDelta t}{ \alpha\overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert _{-1}^{2}\text{ .} \end{align} (A6)Stability follows if the following conditions hold:   $$(1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu} -(1-\sigma)\frac{1}{2}\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{\overline{\nu}_{\textrm{min}}} -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \geqslant0,$$ (A7)  $$\frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\geqslant 0.$$ (A8)Using assumption (2.3), we have   \begin{equation*} \frac{\sigma+1}{2}\sqrt{\mu}-(1+\sigma)\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}} = \frac{2+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}-\frac{\|\nu_j-\overline{\nu}\|_{\infty}}{2\overline{\nu}_{\textrm{min}}}\right) \geqslant 0 \end{equation*}so that (A7) holds. Together with assumption (2.2), we then have   \begin{align*} &(1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu} -(1-\sigma)\frac{1}{2}\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{\overline{\nu}_{\textrm{min}}} -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \\ &\qquad\geqslant (1-\sigma)\frac{\sigma+1}{2\sigma}\sqrt{\mu}-(1-\sigma)\frac{1}{2}\sqrt{\mu}-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n} -\overline{u}_{h}^{n})\Vert^{2} \\ &\qquad=\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}=0, \end{align*}so (A8) holds. Therefore, assuming that (2.2) and (2.3) hold, (A6) reduces to   \begin{align} \frac{1}{2}\Vert u_{j,h}^{n+1}\Vert^{2}&-\frac{1}{2}\Vert u_{j,h}^{n}\Vert^{2} +\frac{1}{4}\Vert u_{j,h}^{n+1}-u_{j,h}^{n}\Vert^{2} +\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ &+\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \left(\Vert\nabla u_{j,h}^{n+1}\Vert^{2}-\Vert\nabla u_{j,h}^{n}\Vert^{2}\right) \leqslant \frac{2\varDelta t}{ \overline{\nu}_{\textrm{min}}} \Vert \,f_{j}^{n+1}\Vert_{-1}^{2}\text{ .} \end{align} (A9)Summing (A9) from n = 0 to n = N − 1 results in   \begin{align} \frac{1}{2}\|u_{j,h}^{N}\|^{2} &+\frac{1}{4}\sum_{n=0}^{N-1}\|u_{j,h}^{n+1}-u_{j,h}^{n}\|^{2} +\sum_{n=0}^{N-1}\frac{\varepsilon(2-\sqrt{\mu})}{4(\sqrt{\mu}+\varepsilon)}\overline{\nu}_{\textrm{min}}\varDelta t\Vert\nabla u_{j,h}^{n+1}\Vert^{2}\nonumber\\ & +\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon} -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \|\nabla u_{j,h}^{N}\|^{2}\nonumber\\ \leqslant& \sum_{n=0}^{N-1}\frac{2\varDelta t}{\overline{\nu}_{\textrm{min}}}\|f_{j}^{n+1}\|_{-1}^{2}+ \frac{1} {2}\|u_{j,h}^{0}\|^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon} -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\|\nabla u_{j,h}^{0}\|^{2}. \end{align} (A10)This completes the proof of stability. Appendix B. Proof of Theorem 3.1 The weak solution of the NSE $$u_{j}$$ satisfies   \begin{align} \left(\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}, v_{h}\right)&+ b^{*}(u_{j}^{n+1}, u_{j}^{n+1}, v_{h})+ (\nu_j\nabla u_{j}^{n+1}, \nabla v_{h})\nonumber\\ &- (p_{j} ^{n+1},\nabla\cdot v_{h})=(\,f_{j}^{n+1}, v_{h}) + \textrm{Intp}(u_{j}^{n+1};v_{h})\qquad\forall\, v_{h}\in V_{h}, \end{align} (B1)where $$\textrm{Intp}(u_{j}^{n+1};v_{h})= \left (\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t} (t^{n+1}),v_{h}\right)$$. Let   $$e_{j}^{n}=u_{j}^{n}-u_{j,h}^{n}=(u_{j}^{n}-R_{h} u_{j}^{n})+(R_{h} u_{j} ^{n}-u_{j,h}^{n})=\eta_{j}^{n}+\xi_{j,h}^{n},$$where $$R_{h} u_{j}^{n} \in V_{h}$$ is the Ritz projection of $$u_{j}^{n}$$ onto $$V_{h}$$ defined by $$\big (\overline{\nu }~\nabla (u_{j}^{n}-R_{h} u_{j}^{n}), \nabla v_h\big ) = 0$$. The associated Ritz projection error satisfies (see, e.g., Thomée, 1997)   $$\|R_h u_j^n - u_j^n\| + h \|\nabla (R_h u_j^n - u_j^n)\| \leqslant C h^{k+1} \|u_j^n\|_{k+1}.$$ (B2)Subtracting (3.1) from (B1) gives   \begin{align*} \left(\frac{\xi_{j,h}^{n+1}-\xi_{j,h}^{n}}{ \varDelta t},v_{h}\right) &+(\overline{\nu}\nabla\xi_{j,h}^{n+1},\nabla v_{h}) +((\nu_j-\overline{\nu})\nabla (u_{j}^{n+1}-u_{j}^n),\nabla v_{h}) +((\nu_j-\overline{\nu})\nabla \xi_{j,h}^n,\nabla v_h) \nonumber\\ & +b^{*}(u_{j}^{n+1},u_{j}^{n+1},v_{h})-b^{*}(\overline{u}_{h}^{n},u_{j,h}^{n+1},v_{h}) -b^{*}(u_{j,h}^{n}-\overline{u}_{h}^n,u_{j,h} ^{n},v_{h})-(p_{j}^{n+1},\nabla\cdot v_{h}) \nonumber\\ =&-\left(\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{\varDelta t},v_{h}\right) -(\overline{\nu}\nabla\eta_{j}^{n+1},\nabla v_{h})-((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla v_{h}) +\textrm{Intp}(u_j^{n+1};v_{h}).\nonumber \end{align*}Setting $$v_{h}=\xi _{j,h}^{n+1}\in V_{h}$$ and rearranging the nonlinear terms, we have   $$\begin{split} &\frac{1}{\varDelta t}\left(\frac{1}{2}\|\xi_{j,h}^{n+1}\|^{2}-\frac{1}{2}\|\xi _{j,h}^{n}\|^{2}+\frac{1}{2}\|\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\|^{2}\right)+ \|\overline{\nu}^{\frac{1}{2}}\nabla\xi_{j,h}^{n+1}\|^{2}\\ &\qquad=-\,((\nu_j-\overline{\nu})\nabla (u_{j}^{n+1}-u_{j}^n),\nabla \xi_{j,h}^{n+1}) -((\nu_j-\overline{\nu})\nabla \xi_{j,h}^n,\nabla \xi_{j,h}^{n+1}) -((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla \xi_{j,h}^{n+1}) \\ &\qquad\quad-b^{*}(u_{j,h}^n-\overline{u}_h^n,u_{j,h}^{n+1}-u_{j,h}^{n},\xi_{j,h}^{n+1}) -b^{*}(u_{j}^{n+1},u_{j}^{n+1},\xi_{j,h}^{n+1})+b^{*}(u_{j,h}^{n}, u_{j,h}^{n+1},\xi_{j,h}^{n+1})\\ &\qquad\quad+(p_{j}^{n+1},\nabla\cdot\xi_{j,h}^{n+1})-\left(\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t},\xi_{j,h}^{n+1}\right) +\textrm{Intp}(u_j^{n+1};\xi_{j,h}^{n+1}). \end{split}$$ (B3)We first bound the viscous terms on the RHS of (B3):   \begin{align} -((\nu_j-\overline{\nu}) \nabla (u_{j}^{n+1}-u_j^n), \nabla\xi_{j,h}^{n+1}) &\leqslant \| \nu_j-\overline{\nu}\|_{\infty} \|\nabla (u_j^{n+1}-u_j^n)\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\|\nabla(u_{j}^{n+1}-u_j^n)\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber \\ &\leqslant \frac{\varDelta t}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \left(\int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2}\, \textrm{d}t\right)+ C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B4)  \begin{align} -((\nu_j-\overline{\nu})\nabla\eta_{j}^{n},\nabla\xi_{j,h}^{n+1}) &\leqslant\| \nu_j-\overline{\nu}\|_{\infty} \|\nabla\eta_{j} ^{n}\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_0}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B5)  \begin{align} -(\nu_j-\overline{\nu})(\nabla\xi_{j,h}^{n},\nabla\xi_{j,h}^{n+1}) &\leqslant\| \nu_j-\overline{\nu}\|_{\infty} \|\nabla\xi_{j,h}^{n}\| \|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{1}{4C_1}\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \|\nabla\xi_{j,h}^{n}\|^{2} + C_1\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber \\ &\leqslant \frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2} \|\nabla\xi_{j,h}^{n}\|^{2}+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2} \|\nabla\xi_{j,h}^{n+1}\|^{2}, \end{align} (B6)in which we note that both terms on the RHS of (B6) need to be hidden in the LHS of the error equation, thus $$C_1= \frac{\|\nu _j-\overline{\nu }\|_{\infty }}{2\overline{\nu }_{\textrm{min}}}$$ is selected to minimize the summation. Next we analyse the nonlinear terms on the RHS of (B3) one by one. For the first nonlinear term we have   \begin{align}\!\! -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j,h}^{n+1}-u_{j,h}^n, \xi_{j,h}^{n+1}) =&-\,b^{\ast}(u_{j,h}^n-\overline{u}_h^n, e_{j}^{n+1}-e_{j}^n, \xi_{j,h}^{n+1}) +b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j}^{n+1}-u_{j}^n, \xi_{j,h}^{n+1})\nonumber\\ =&-\,b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n+1}, \xi_{j,h}^{n+1})+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n}, \xi_{j,h}^{n+1})\nonumber\\ &+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \xi_{j}^{n}, \xi_{j,h}^{n+1})+b^{\ast}(u_{j,h}^n-\overline{u}_h^n, u_{j}^{n+1}-u_{j}^n, \xi_{j,h}^{n+1})\,. \end{align} (B7)Using inequality (1.8) and Young’s inequality, we have the estimates   \begin{align} -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla\eta_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla\eta_{j}^{n+1}\|^{2} +C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B8)and   \begin{align} -b^{\ast}(u_{j,h}^n-\overline{u}_h^n, \eta_{j}^{n}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla\eta_{j}^{n}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1} \|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B9)Because $$b^{\ast }(\cdot ,\cdot ,\cdot )$$ is skew-symmetric we have   \begin{align*} b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n},\xi_{j,h}^{n+1})&=b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi _{j,h}^{n}-\xi_{j,h}^{n+1},\xi_{j,h}^{n+1}) \\& =b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n+1},\xi_{j,h}^{n+1}-\xi_{j,h}^n)\,. \end{align*}Then by inequality (1.8) we obtain   \begin{align} b^{\ast}(u_{j,h}^n-\overline{u}_h^n,\xi_{j,h}^{n},\xi_{j,h}^{n+1}) &\leqslant \Vert \nabla (u_{j,h}^n-\overline{u}_h^n)\Vert\Vert \nabla \xi_{j,h}^{n}\Vert \Vert \nabla (\xi_{j,h}^{n+1}-\xi_{j,h}^n)\Vert^{1/2}\Vert \xi_{j,h}^{n+1}-\xi_{j,h}^n\Vert^{1/2}\nonumber\\ &\leqslant C\Vert \nabla (u_{j,h}^n-\overline{u}_h^n)\Vert\Vert \nabla \xi_{j,h}^{n}\Vert h^{-1/2}\Vert \xi_{j,h}^{n+1}-\xi_{j,h}^n\Vert\nonumber\\ &\leqslant \frac{1}{4\varDelta t}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2}+\left[ C\frac{\varDelta t}{h}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}. \end{align} (B10)For the last term of (B7) we have   \begin{align} b^{*}(u_{j,h}^n-\overline{u}_h^n,u_{j}^{n+1}-u_{j}^{n},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|\|\nabla(u_{j}^{n+1}-u_{j}^{n})\|\|\nabla\xi_{j,h}^{n+1} \|\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\|\nabla (u_{j}^{n+1}-u_{j}^{n})\|^{2} +C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &\leqslant \frac{C^2}{4C_0}\varDelta t \overline{\nu}_{\textrm{min}}^{-1}\|\nabla (u_{j,h}^n-\overline{u}_h^n)\|^{2}\left(\int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \, \mathrm{d}t\right) +C_0\overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} .\nonumber \end{align} (B11)Next we bound the last two nonlinear terms on the RHS of (B3) as   \begin{align} &-b^{*}(u_{j}^{n+1}, u_{j}^{n+1},\xi_{j,h}^{n+1}) +b^{*}(u_{j,h}^{n},u_{j,h}^{n+1},\xi_{j,h}^{n+1})\nonumber\\ &\qquad= -b^{*}(e_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u_{j,h}^{n},e_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u^{n+1}_{j}-u_j^{n}, u_{j}^{n+1}, \xi_{j,h}^{n+1}) \nonumber\\ &\qquad= -b^{*}(\eta^{n}_j,u_{j}^{n+1},\xi_{j,h}^{n+1})-b^{*}(\xi_{j,h}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) -b^{*}(u^{n}_{j,h}, \eta_{j}^{n+1}, \xi_{j,h}^{n+1})-b^{*}(u^{n+1}_{j}-u^{n}_j, u_{j}^{n+1}, \xi_{j,h}^{n+1}), \end{align} (B12)where, with the assumption $$u_j^{n+1}\in L^{\infty }(0,T; H^1(\varOmega ))$$, we have   \begin{align} -b^{*}(\eta_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla \eta_{j}^{n}\|\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\| \nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla\eta_{j}^{n}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B13)Using inequality (1.8), Young’s inequality and $$u_j^{n+1}\in L^{\infty }(0,T; H^1(\varOmega ))$$ we get   \begin{align} -b^{*}(\xi^{n}_{j,h}, u_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\| \nabla\xi ^{n}_{j,h} \|^{1/2} \Vert\xi^{n}_{j,h}\Vert^{1/2}\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\leqslant C\| \nabla\xi^{n}_{j,h} \|^{1/2} \Vert\xi^{n}_{j,h}\Vert ^{1/2}\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant C\Big(\frac{1}{4\alpha}\|\nabla\xi^{n}_{j,h} \|\|\xi^{n}_{j,h} \| + \alpha\|\nabla\xi_{j,h}^{n+1}\|^{2}\Big)\\ &\leqslant C\Big[\frac{1}{4\alpha} \big(\frac{\varDelta}{2}\|\nabla\xi^{n}_{j,h} \|^{2}+\frac{1}{2\varDelta}\|\xi^{n}_{j,h} \|^2 \big) + \alpha\|\nabla\xi_{j,h}^{n+1}\|^{2} \Big]\nonumber\\ &\leqslant C_0\overline{\nu}_{\textrm{min}} \|\nabla \xi^{n}_{j,h} \|^{2}+\frac{C^4}{64C_0^3\overline{\nu}_{\textrm{min}}^{3}}\|\xi^{n}_{j,h} \|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2}\,,\nonumber \end{align} (B14)where we set $$\alpha = \frac{C_0\overline{\nu }_{\textrm{min}}}{C}$$ and $$\varDelta = \frac{8C_0^2\overline{\nu }_{\textrm{min}}^2}{C^2}$$. By Young’s inequality, (1.8) and the result (A10) from the stability analysis, i.e., $$\Vert u_{j,h}^n \Vert ^2 \leqslant C$$, we also have   \begin{align} b^{*}(u^{n}_{j,h}, \eta_{j}^{n+1}, \xi_{j,h}^{n+1}) &\leqslant C\|\nabla u^{n}_{j,h}\|^{1/2}\Vert u_{j,h}^n\Vert^{1/2}\|\nabla \eta_{j}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla u_{j,h}^{n}\|\|\nabla\eta^{n+1}_{j}\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B15)and   \begin{align} b^{*}(u_{j}^{n+1}-u_{j}^{n},u_{j}^{n+1},\xi_{j,h}^{n+1}) &\leqslant C\|\nabla (u_{j}^{n+1}-u_{j}^{n})\|\|\nabla u_{j}^{n+1}\|\|\nabla\xi_{j,h} ^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\|\nabla(u_{j}^{n+1}-u_{j}^{n})\|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &= \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1} \varDelta t^2 \left\|\frac{\nabla u_{j}^{n+1}- \nabla u_{j}^{n}}{\varDelta t}\right\|^{2} + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2} \\ &= \frac{C^2}{4C_0} \overline{\nu}_{\min}^{-1} \varDelta t^2 \int_{\varOmega}\left(\frac{1}{\varDelta t}\int_{t^{n}}^{t^{n+1}} \nabla u_{j,t} \, \mathrm{d}t\right)^{2} \textrm{d}\varOmega+ C_0 \overline{\nu} _{\min}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t} \|^{2}\, \mathrm{d}t + C_0 \overline{\nu}_{\textrm{min}} \|\nabla\xi_{j,h}^{n+1}\|^{2}. \nonumber \end{align} (B16)For the pressure term in (B3), because $$\xi _{j,h}^{n+1}\in V_{h}$$ we have, for any $$q_{j, h}^{n+1}\in Q_h$$,   \begin{align} (p_{j}^{n+1},\nabla\cdot\xi_{j,h}^{n+1})&=(p_{j}^{n+1}-q_{j,h}^{n+1}, \nabla\cdot\xi_{j,h}^{n+1})\leqslant\sqrt{d}\,\|p_{j}^{n+1}-q_{j,h}^{n+1}\|\|\nabla\xi_{j,h}^{n+1}\|\\ &\leqslant \frac{1}{4d\, C_0}\overline{\nu}_{\textrm{min}}^{-1} \|p_{j}^{n+1}-q_{j,h}^{n+1}\|^{2} + C_0\,\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \text{ .}\nonumber \end{align} (B17) The other terms are bounded as   \begin{align} \Big( \frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t},\xi_{j,h}^{n+1} \Big) &\leqslant C \Big\|\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t} \Big\| \|\nabla\xi_{j,h} ^{n+1}\|\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{\eta_{j}^{n+1}-\eta_{j}^{n}}{ \varDelta t} \right\|^{2} +C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{1}{\varDelta t}\int_{t^{n}}^{t^{n+1}} \eta_{j,t} \textrm{ } \mathrm{d}t \right \|^2+C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\| \eta_{j,t}\|^{2}\textrm{ } \mathrm{d}t+C_0 \overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \end{align} (B18)and   \begin{align} \textrm{Intp}(u_{j}^{n+1};\xi_{j,h}^{n+1}) &=\left(\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1}),\xi_{j,h}^{n+1}\right)\leqslant C\left\|\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1})\right\| \|\nabla \xi_{j,h}^{n+1}\|\nonumber\\ &\leqslant \frac{C^2}{4C_0} \overline{\nu}_{\textrm{min}}^{-1}\left \|\frac{u_{j}^{n+1}-u_{j}^{n}}{\varDelta t}-u_{j,t}(t^{n+1}) \right\|^{2} + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2} \nonumber \\ &\leqslant \frac{C^2}{4C_0}\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\|u_{j,tt}\|^{2}\, \mathrm{d}t + C_0\overline{\nu}_{\textrm{min}}\|\nabla\xi_{j,h}^{n+1}\|^{2}. \end{align} (B19)Combining (B4)–(B19), we have   \begin{align}&\frac{1}{\varDelta t}\left(\frac{1}{2}||\xi_{j,h}^{n+1}||^{2}-\frac{1}{2}||\xi _{j,h}^{n}||^{2}+\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \right)\nonumber\\ &\quad+ C_0 \overline{\nu}_{\textrm{min}} \Vert \nabla \xi_{j,h}^{n+1}\Vert^2+ \overline{\nu}_{\textrm{min}} \left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)\left(\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2}-\Vert\nabla \xi_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[(1-\sigma)\left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) -\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\sigma\left(1-14C_0-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla \xi_{j,h}^n\Vert^2\nonumber\\ &\leqslant C \bigg[ \frac{1}{\overline{\nu}_{\textrm{min}}^{3}}\Vert\xi_{j,h}^{n}\Vert^{2} +\varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \mathrm{d}t \\ &\quad+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2} \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1} }\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla u_{j,h}^{n}\Vert\Vert\nabla\eta_{j}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ &\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\mathrm{d}t \nonumber \bigg]. \end{align} (B20)Note that the generic constant C independent of $$\varDelta t$$ is used on the RHS. It, however, depends on the geometry and mesh due to the use of an inverse inequality in (B10). Similarly to the stability analysis we take $$C_0= \frac{1}{14}(1-\frac{\sigma +1}{2\sigma }\sqrt{\mu }) = \frac{1}{14} \frac{\varepsilon }{\sqrt{\mu }+\varepsilon }( 1-\frac{\sqrt{\mu }}{2})$$ with $$\sigma =\frac{\sqrt{\mu }+\varepsilon }{2-\sqrt{\mu }}$$. Then (B20) becomes   \begin{align}& \frac{1}{\varDelta t}\left(\frac{1}{2}||\xi_{j,h}^{n+1}||^{2}-\frac{1}{2}||\xi _{j,h}^{n}||^{2}+\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \right)\nonumber\\&\quad+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right) \overline{\nu}_{\textrm{min}} \Vert \nabla \xi_{j,h}^{n+1}\Vert^2 \nonumber\\&\quad+ \overline{\nu}_{\textrm{min}} \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\left(\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2}-\Vert\nabla \xi_{j,h}^{n}\Vert ^{2}\right)\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}-\overline{u}_{h}^{n})\Vert^{2}\right] \Vert\nabla \xi_{j,h}^{n}\Vert^{2}\nonumber\\ &\quad+\overline{\nu}_{\textrm{min}} \left[\frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\Vert \nabla \xi_{j,h}^n\Vert^2\nonumber\\ \leqslant&\ C\bigg[ \frac{1}{\overline{\nu}_{\textrm{min}}^{3}}\Vert\xi_{j,h}^{n}\Vert^{2} +\varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \textrm{ }\mathrm{d}t +\overline{\nu}_{\textrm{min}}\Vert\nabla\eta_{j}^{n+1}\Vert^{2} \\ \quad&\quad+ \frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2} \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1} }\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla u_{j,h}^{n}\Vert\Vert\nabla\eta_{j}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \nonumber\\ \quad&\quad+ \overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\mathrm{d}t\nonumber \bigg]. \end{align} (B21)By the convergence condition (3.3) we have   \begin{equation*} \frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\geqslant \frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\sqrt{\mu}}{2}> 0 \end{equation*}and   \begin{align*}& \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) -\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\\ &\qquad\geqslant \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\sqrt{\mu}}{2}\right)-\frac{\sqrt{\mu}}{2}\\ &\qquad\geqslant \frac{\sqrt{\mu}+\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2-\sqrt{\mu}}{\sqrt{\mu}+\varepsilon}\right)-\frac{\sqrt{\mu}}{2} \geqslant \frac{\sqrt{\mu}}{2}-\frac{\sqrt{\mu}}{2}=0 \end{align*}and by the convergence conditions (3.2) and (3.3) we have   \begin{align*} &\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2+\varepsilon}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right)-\frac{C \varDelta t}{\overline{\nu}_{\textrm{min}} h}\Vert\nabla(u_{j,h}^{n}- \overline{u}_{h}^{n})\Vert^{2}\\ &\qquad\geqslant\frac{2-2\sqrt{\mu}-\varepsilon}{2-\sqrt{\mu}}\left(\frac{\sqrt{\mu}}{2}\frac{2-\sqrt{\mu}}{\sqrt{\mu}+\varepsilon}\right)-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}\\ &\qquad\geqslant \frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}-\frac{(2-2\sqrt{\mu}-\varepsilon)\sqrt{\mu}}{2(\sqrt{\mu}+\varepsilon)}\geq 0. \end{align*}Summing (B20) from n = 1 to N − 1 and multiplying both sides by $$\varDelta t$$ gives   \begin{align*}& \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}+\frac{1}{4}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n+1} -\xi_{j,h}^{n}\Vert^{2}\nonumber\\&\quad+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2} \\ &\quad+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{N}\|^{2}\\ \leqslant&\ \frac{1}{2}\|\xi_{j,h}^{0}\|^{2}+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} +\frac{C\varDelta t}{\overline{\nu}_{\textrm{min}}^{3}}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n}\Vert^{2}\\ &\quad+C\varDelta t\sum_{n=0}^{N-1} \bigg[ \varDelta t\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}} \int_{t^{n}}^{t^{n+1}}\| \nabla u_{j,t}\|^{2} \,\mathrm{d}t +\overline{\nu}_{\textrm{min}}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\\ &\qquad\qquad\qquad+\frac{\| \nu_j-\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_{h}^n)\Vert^{2} \Vert\nabla\eta_{j}^{n+1}\Vert^{2}\\ &\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \Vert\nabla (u_{j,h}^n-\overline{u}_h^n)\Vert^{2} \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\, \mathrm{d}t \\ \quad&\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n}\Vert^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t \int_{t^{n}}^{t^{n+1}}\Vert\nabla u_{j,t}\Vert^{2}\,\mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert \\ \quad&\qquad\qquad\qquad+\overline{\nu}_{\textrm{min}}^{-1}\Vert p_{j}^{n+1}-q_{j,h}^{n+1}\Vert^{2} + \overline{\nu}_{\textrm{min}}^{-1} \varDelta t^{-1}\int_{t^{n}}^{t^{n+1}}\Vert\eta_{j,t}\Vert^{2}\, \mathrm{d}t +\overline{\nu}_{\textrm{min}}^{-1} \varDelta t \int_{t^{n}}^{t^{n+1}}\Vert u_{j,tt}\Vert^{2}\, \mathrm{d}t \bigg] \text{ .} \end{align*}Using the Ritz projection error (B2) and the result (2.4) from the stability analysis, i.e., $$\varDelta t \sum _{n=0}^{N-1}\Vert \nabla u_{j,h}^{n} \Vert ^2 \leqslant C$$, we have   $$C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t\sum_{n=0}^{N-1}\Vert\nabla\eta_{j}^{n+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert \leqslant C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}\varDelta t \sum_{n=0}^{N-1}\Vert u_j^{n+1}\Vert^2_{k+1}\Vert^{2}\Vert\nabla u_{j,h}^{n}\Vert$$ (B22)  \begin{align} \leqslant &\ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}\left( \varDelta t \sum_{n=0}^{N-1}\Vert u_{j}^{n+1}\Vert_{k+1}^4+\varDelta t \sum_{n=0}^{N-1}\Vert\nabla u_{j,h}^{n}\Vert^2 \right)\\ \leqslant &\ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}.\nonumber \end{align} (B23)By applying again the Ritz projection error (B2) and the interpolation error (1.6) we have   \begin{align} \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}&+ \overline{\nu}_{\textrm{min}}\varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}_{\textrm{min}}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{N}\|^{2} +\sum_{n=0}^{N-1}\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \nonumber\\ & +\frac{1}{15}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2}\nonumber\\ \leqslant\ &\frac{1}{2}\|\xi_{j,h}^{0}\|^{2}+\overline{\nu}_{\textrm{min}} \varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} +\frac{C\varDelta t}{\overline{\nu}_{\textrm{min}}^{3}}\sum_{n=0}^{N-1}\Vert\xi_{j,h}^{n}\Vert^{2} \\ &+ C\varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}\Vert\nabla u_{j,t}\Vert^2_{2,0} +C\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber \\ &+Ch^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t \Vert\nabla u_{j,t}\Vert^2_{2,0}+ C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ & +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2\Vert\nabla u_{j,t}\Vert^2_{2,0} + C\overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+C\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \nonumber\\ &+C\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2}+C\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}\Vert u_{j,t}\Vert_{2,k+1}^{2} +C\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}\Vert u_{j,tt}\Vert_{2,0}^{2}.\nonumber \end{align} (B24)The next step is the application of the discrete Gronwall inequality (see Girault & Raviart, 1979, p. 176):   \begin{align} \frac{1}{2}\|\xi_{j,h}^{N}\|^{2}&+ \overline{\nu}_{\textrm{min}}\varDelta t\left(\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right) \|\nabla\xi_{j,h}^{N}\|^{2} +\sum_{n=0}^{N-1}\frac{1}{4}\Vert\xi_{j,h}^{n+1}-\xi_{j,h}^{n}\Vert^{2} \nonumber\\ &+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}}\varDelta t\sum_{n=0}^{N-1}\|\nabla\xi_{j,h}^{n+1}\|^{2}\\ \leqslant& \ e^{\frac{CT}{\overline{\nu}_{\textrm{min}}^{3}}} \bigg\{\frac{1}{2}\|\xi_{j,h}^{0}\|^{2} +\overline{\nu}_{\textrm{min}}\varDelta t \left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\|\nabla\xi_{j,h}^{0}\|^{2} \nonumber\\ &\qquad\quad+C\Big[ \varDelta t^2\frac{\| \nu_j -\overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}{\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} +\overline{\nu}_{\textrm{min}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} +\frac{\| \nu_j - \overline{\nu}\|_{\infty}^2}{\overline{\nu}_{\textrm{min}}}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ &\qquad\quad+h^{2k+1}\varDelta t^{-1}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1}+C h \varDelta t {\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} + \overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^2_{2,k+1} \nonumber\\ & \qquad\quad+\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^2{\left\vert\!\left\vert\!\left\vert \nabla u_{j,t} \right\vert\!\right\vert\!\right\vert}^2_{2,0} + \overline{\nu}_{\textrm{min}}^{-1}h^{2k}{\left\vert\!\left\vert\!\left\vert u_j \right\vert\!\right\vert\!\right\vert}^4_{4, k+1}+\overline{\nu}_{\textrm{min}}^{-1}h^{2k} \nonumber \\ &\qquad\quad+\overline{\nu}_{\textrm{min}}^{-1} h^{2s+2}{\left\vert\!\left\vert\!\left\vert p_{j} \right\vert\!\right\vert\!\right\vert}_{2,s+1}^{2} +\overline{\nu}_{\textrm{min}}^{-1} h^{2k+2}{\left\vert\!\left\vert\!\left\vert u_{j,t} \right\vert\!\right\vert\!\right\vert}_{2,k+1}^{2} +\overline{\nu}_{\textrm{min}}^{-1}\varDelta t^{2}{\left\vert\!\left\vert\!\left\vert u_{j,tt} \right\vert\!\right\vert\!\right\vert}_{2,0}^{2}\Big] \bigg\}.\nonumber \end{align} (B25) Recall that $$e_{j}^{n}=\eta _{j}^{n}+\xi _{j,h}^{n}$$. Using the triangle inequality on the error equation to split the error terms into terms of $$\eta _{j}^{n}$$ and $$\xi _{j,h}^{n}$$ gives   \begin{align*} \frac{1}{2}\Vert e_{j}^{N}\Vert^{2} &+\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla e_{j}^{n+1}\Vert^{2} \\ \leqslant&\ \frac{1}{2}\Vert\xi_{j,h}^{N}\Vert^{2} +\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla \xi_{j,h}^{n+1}\Vert^{2} \\ &+\frac{1}{2}\Vert\eta_{j}^{N}\Vert^{2} +\frac{1}{14}\frac{\varepsilon}{\sqrt{\mu}+\varepsilon}\left(1-\frac{\sqrt{\mu}}{2}\right)\overline{\nu}_{\textrm{min}} \varDelta t\sum_{n=0}^{N-1}\Vert\nabla \eta_{j}^{n+1}\Vert ^{2} \end{align*}and   \begin{equation*} \begin{aligned} \frac{1}{2}\Vert\xi_{j,h}^{0}\Vert^{2} &+\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\| \nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla \xi_{j,h} ^{0}\Vert^{2}\\ \leqslant&\ \frac{1}{2}\Vert e_{j}^{0}\Vert^{2} +\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla e_{j}^{0}\Vert^{2}\\ &+\frac{1}{2}\Vert\eta_{j}^{0}\Vert^{2} +\left[\frac{\sqrt{\mu}}{2}\frac{(2+\varepsilon)}{\sqrt{\mu}+\varepsilon}-\frac{\nu_j-\overline{\nu}\|_{\infty}}{2 \overline{\nu}_{\textrm{min}}}\right]\overline{\nu}_{\textrm{min}}\varDelta t \Vert\nabla \eta_{j}^{0}\Vert^{2}\, . \end{aligned} \end{equation*}Applying inequality (B25), using the previous bounds for $$\eta _{j}^{n}$$ terms, and absorbing constants into a new constant C, completes the proof of Theorem 3.1. © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: May 24, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off