\section{Introduction}\label{sec-1}
	We consider the  following time-dependent stochastic Navier-Stokes equations:
	\begin{subequations}\label{eq1.1}
		\begin{alignat}{2} \label{eq1.1a}
			d{\bf u} &=\bigl[\nu\Delta {\bf u} - {\bf u}\cdot \nabla{\bf u} -\nabla p + { {\bf f}}\bigr] dt +  {\bf g}(t)d W(t)  &&\qquad\mbox{a.s. in}\, D_T,\\
			\div {\bf u} &=0 &&\qquad\mbox{a.s. in}\, D_T,\label{eq1.1b}\\

			{\bf u}(0)&= {\bf u}_0 &&\qquad\mbox{a.s. in}\, D,\label{eq1.1d}
		\end{alignat}
	\end{subequations}
	where $D = (0,L)^2 \subset \mathbb{R}^2 \,$ represents a period of the periodic domain in $\mathbb{R}^2$, ${\bf u}$ and $p$ stand for respectively the velocity field and the pressure of the fluid, $\{W(t); t\geq 0\}$ denotes real-valued
	Wiener process, and ${\bf f}$ is a body force function. In addition, ${\bf g}$ is the diffusion coefficient (see section \ref{sub5.2.2} for its precise definition).  Here we seek periodic-in-space solutions $({\bf u},p)$ with period $L$, that is,
	${\bf u}(t,{\bf x} + L{\bf e}_i) = {\bf u}(t,{\bf x})$ and $p(t,{\bf x}+L{\bf e}_i)=p(t,{\bf x})$
	almost surely
	and for any $(t, {\bf x})\in (0,T)\times \mathbb{R}^d$  and $1\leq i\leq d$, where
	$\{\bf e_i\}_{i=1}^d$ denotes the canonical basis of $\mathbb{R}^d$.

	  Numerical analysis of \eqref{eq1.1} has been studied by several researchers.  In  \cite{BCP12} the authors established stability and convergence of the standard mixed finite element method of \eqref{eq1.1}. Later in \cite{CP2012} the authors proved the rates of convergence in probability for the velocity approximation in the case of multiplicative noise. The main difficulty for establishing a strong convergence for any numerical approximation of \eqref{eq1.1} is the interplay between nonlinearity and  stochasticity of the equations. To compute or estimate quantities of stochastic interests such as the expectation and moments, all the norms must have another layer of the integration which is the main reason why the classical Gronwall inequality argument fails. To overcome this difficulty,  in \cite{CP2012} the authors introduced a sequence of sub-sample spaces which converges to the sample space under the probability measure. The error estimates on these sub-sample spaces are computed with partial expectations and are considered as weak convergence, and some improved error estimates of the same type were recently obtained in \cite{Qui_2022}.
	  %
	  In \cite{BM18, BM_space} the authors showed the strong $L^2$-convergence of the mixed finite element method for \eqref{eq1.1} by estimating the error estimates on the complements of these sub-sample spaces, which then leads to  the strong convergence with a logarithmic rate.  Moreover, the authors  were able to establish in \cite{BM2021}  strong convergence with a polynomial rate  in the case of a divergence-free additive noise.
	  We note that all the above mentioned error estimates  for \eqref{eq1.1}  are second moment estimates and most  are only for the velocity approximation.
	  No high moment and pathwise error estimates have been reported in the literature so far.  These missing error estimates are important to know because they provide different
	  quantities of stochastic interests in practice.

	  The primary goals of this paper is to fill such a void in the case of general additive noise and to develop the analysis techniques for deriving high moment and pathwise error estimates for numerical nonlinear stochastic PDEs in general. It should be noted that
	  the desired high moment and pathwise error estimates will be obtained for both the velocity and pressure approximations of the fully discrete mixed finite element method for \eqref{eq1.1}.  Our main ideas are to obtain the former based on an exponential
	  stability estimate, which is inspired by a similar idea first introduced in  \cite{BM2021}, and a bootstrap technique,  and to obtain the latter by using the
	  Kolmogorov Theorem (see Theorem \ref{kolmogorov}).

	 The remainder of this paper is organized as follows. In section \ref{sec-2}, we present some preliminaries including the definition
	 of variational solutions to \eqref{eq1.1} and the assumptions on the diffusion function ${\bf g}$. In section \ref{section5.3},  we introduce  the time discretization for \eqref{eq1.1} in Algorithm 1 and establish some stability estimates for its solution, including an exponential stability estimate which plays a crucial role in our error analysis for the velocity approximation in Theorems \ref{theorem_semi_chapter5}, and \ref{thm_pathwise_semi}, and the error estimates for the pressure approximation in Theorems \ref{semi_pressure_error} and \ref{thm_semi_pathwise_pressure}.  In section \ref{section5.4}, we formulate the fully discrete mixed finite element method in Algorithm 2. The highlight of this section is to establish the desired  high moment and pathwise error estimates for both velocity and pressure approximations. Finally,
	 we present three numerical experiment results in section \ref{sec-5} to validate the error estimates and their sharpness.

	\section{Preliminaries}\label{sec-2}
	\subsection{Notations}\label{sec-2.1}
	Standard function and space notation will be adopted in this paper.
	Let ${\bf H}^1_0(D)$ denote the subspace of ${\bf H}^1(D)$ whose ${\mathbb R}^d$-valued functions have zero trace on $\partial D$, and $(\cdot,\cdot):=(\cdot,\cdot)_D$ denote the standard $L^2$-inner product, with induced norm $\Vert \cdot \Vert$. We also denote ${\bf L}^p_{per}(D)$ and ${\bf H}^{k}_{per}(D)$ as the Lebesgue and Sobolev spaces of the functions that are periodic.

	 $C$ denotes a generic constant which is independent of the mesh parameter $h$ and $k$.

	Let $(\Omega,\mathcal{F}, \{\mathcal{F}_t\},\mathbb{P})$ be a filtered probability space with the probability measure $\mathbb{P}$, the
	$\sigma$-algebra $\mathcal{F}$ and the continuous  filtration $\{\mathcal{F}_t\} \subset \mathcal{F}$. For a random variable $v$
	defined on $(\Omega,\mathcal{F}, \{\mathcal{F}_t\},\mathbb{P})$,
	${\mathbb E}[v]$ denotes the expected value of $v$.
	For a vector space $X$ with norm $\|\cdot\|_{X}$,  and $1 \leq p < \infty$, we define the Bochner space
	$\bigl(L^p(\Omega;X); \|v\|_{L^p(\Omega;X)} \bigr)$, where
	$\|v\|_{L^p(\Omega;X)}:=\bigl({\mathbb E} [ \Vert v \Vert_X^p]\bigr)^{\frac1p}$.
	We also define
	\begin{align*}
		{\mathbb H} := \bigl\{{\bf v}\in  {\bf L}^2_{per}(D) ;\,\div {\bf v}=0 \mbox{ in }D\, \bigr\}\, , \quad
		{\mathbb V} :=\bigl\{{\bf v}\in  {\bf H}^1_{per}(D) ;\,\div {\bf v}=0 \mbox{ in }D \bigr\}\, .
	\end{align*}

	We recall from \cite{Girault_Raviart86} the (orthogonal) Helmholtz projection
	${\bf P}_{{\mathbb H}}: {\bf L}^2_{per}(D) \rightarrow {\mathbb H}$  and define the Stokes operator ${\bf A} := -{\bf P}_{\mathbb H} \Delta: {\mathbb V} \cap {\bf H}^2_{per}(D) \rightarrow {\mathbb H}$.

	\subsection{Some useful facts and inequalities} In this subsection, we cite some   useful facts and inequalities which will  be used in later sections.

	First of all, we recall the Kolmogorov Criteria/Theorem for a pathwise continuity of stochastic processes (cf. \cite{PZ1992}).

	\begin{theorem}[Kolmogorov's criteria]\label{kolmogorov}
		Let ${\bf X}(t), t \in [0,T]$, be a stochastic process with values in a separable Banach space $E$ such that, for some positive constant $C>0$, $\alpha > 0, \beta > 0$ and all $t,s \in [0,T]$,
		\begin{align}
			\mathbb{E}\bigl[\|{\bf X}(t) - {\bf X}(s)\|^{\beta}\bigr] \leq C|t-s|^{1+\alpha}.
		\end{align}
		Then for each $T>0$, almost every $\omega$ and each $0< \gamma < \frac{\alpha}{\beta}$ there exists a constant $K = K(\omega,\gamma,T)$ such that
		\begin{align}
			\|{\bf X}(t,\omega) - {\bf X}(s,\omega)\| \leq K|t-s|^{\gamma}\qquad\mbox{ for all } t,s \in [0,T].
		\end{align}
		Moreover, $\mathbb{E}\bigl[|K|^\beta\bigr] < \infty$ for all $\beta >0$.
	\end{theorem}

	\begin{lemma}[Burkholder-Davis-Gundy inequality, \cite{PZ1992}]
		Let $\pmb{\phi}(t)$ be a stochastic process for all $t \in [0,T]$. For any $p >0$, there exists a positive constant $C_b$ such that:
		\begin{align}
			\mathbb{E}\biggl[\max_{0\leq t \leq T}\biggl\|\int_0^t \pmb{\phi}(\xi)\, dW(\xi)\biggr\|_2^p\biggr] \leq C_b\, \mathbb{E}\biggl[\biggl(\int_0^T \|\pmb{\phi}(\xi)\|^2_{L^2}\, d\xi\biggr)^{\frac{p}{2}}\biggr],
		\end{align}
		where $C_b = C_b(T,p)$.
\end{lemma}

	The next lemma recalls the well-known It\^o isometry and also states a helpful inequality for stochastic processes 	(cf. \cite{PZ1992}).

	\begin{lemma}\label{lemma2.3}
		Let $\pmb{\phi}(t)$ be a stochastic process for all $t \in [0,T]$. Define $\displaystyle {\bf X}_t = \int_0^t \pmb{\phi}(\xi)\, dW(\xi)$. We have:
		\begin{enumerate}[{\rm (i)}]
			\item If $\pmb{\phi} \in L^2(\Omega;L^2(0,T; {\bf L}^2(D)))$, then
			\begin{align}\label{ito2}
				\mathbb{E}\bigl[\|{\bf X}_t\|^2_{{\bf L}^2}\bigr] = \mathbb{E}\biggl[\int_0^t\|\pmb{\phi}(\xi)\|^2_{{\bf L}^2}\, d\xi\biggr].
			\end{align}
			\item If $\pmb{\phi} \in L^p(\Omega;L^p(0,T;{\bf L}^2(D)))$ for $p > 2$, then
			\begin{align}\label{ito4}
				\mathbb{E}\bigl[\|{\bf X}_t\|^p_{{\bf L}^2}\bigr] \leq C(t,p)\,\mathbb{E}\biggl[\int_0^t \|\pmb{\phi}(\xi)\|^p_{{\bf L}^2}\, d\xi\biggr],
			\end{align}
			where $\displaystyle C(t,p) = \frac{C_b}{2}(p-1)(p-2) t^{\frac{p}{2}} + (p-1)C_b$.
		\end{enumerate}
	\end{lemma}

The next lemma states  both continuous and discrete Gronwall inequalities.

	\begin{lemma}[Gronwall inequalities]\label{gronwall}
		\begin{enumerate}[{\rm (a)}]
			\item Let $y$ and $g$ be nonnegative and integrable and $c$ a nonnegative constant. If
			\begin{align*}
				y(t) \leq c + \int_0^t g(s)y(s)\, ds\qquad\forall t \geq 0,
			\end{align*}
			then $\displaystyle y(t) \leq c\exp\biggl(\int_0^t g(s)\, ds\biggr)$ for all $t \geq 0$.

			\item Let $\{y_n\}$ and $\{g_n\}$ be nonnegative sequences and $c$ a nonnegative constant. If
			\begin{align*}
				y_n \leq c + \sum_{0\leq j < n} g_jy_j\qquad\forall n \geq 0,
			\end{align*}
			then $\displaystyle y_n \leq c \prod_{0\leq j <n}(1+g_j) \leq c\exp\biggl(\sum_{0\leq j < n}g_j\biggr)$ for all $n \geq 0$.
		\end{enumerate}
	\end{lemma}

	\smallskip

	\begin{lemma}[Hoeffding's lemma]
		Let $X$ be a real-valued random variable with $\mathbb{E}[X] = \mu$, such that $a\leq X \leq b$ a.s., for some $a,b \in \mathbb{R}$. Then, for all $\lambda \in \mathbb{R}$,
		\begin{align}
			\mathbb{E}\bigl[e^{\lambda X}\bigr] \leq \exp\Bigl(\lambda\mu + \frac{\lambda^2}{8}(b-a)^2\Bigr).
		\end{align}
	\end{lemma}

	Finally, we recall the following property of  the $\mathbb{R}$-valued Wiener process.

	\begin{align}\label{mean_wiener}
		\mathbb{E}\Bigl[|W(t) - W(s)|^{2m}\Bigr] \leq C_m|t-s|^{m}\qquad\forall m \in \mathbb{N},
	\end{align}
	where, for $m=1$, we obtain the equality and $C_m =1$. See \cite{Ichikawa} for higher dimensional Wiener process.

 In addition, there exists a deterministic constant $C>0$ such that:
	\begin{align}\label{a.s_wiener}
		|W(t) - W(s)| \leq C|t-s|^{\alpha},\qquad\mathbb{P}-a.s.,
	\end{align}
	where,  when $m=1$,  the equality holds and $C_1 =1$. See \cite{Ichikawa} for high dimensional generalizations.

	\subsection{Variational solutions}\label{sub5.2.2}
	Let ${\bf g}(\cdot,\cdot): D\times [0,T] \rightarrow \mathbb{R}$ and ${\bf g} \in
	L^{\infty}\bigl(0,T;{\bf H}^1(D)\bigr)$. We make the following additional assumptions on ${\bf g}$:
	\begin{enumerate}[(G1)]
		\item ${\bf g}$ is Lipschitz continuous in time, that is, there is a constants $C_{{\bf g}} > 0$ such that
		\begin{align*}
			\|{\bf g}(t) - {\bf g}(s)\|_{{\bf L}^2} \leq C_{{\bf g}}|t-s|\qquad\forall t,s \in [0,T].
		\end{align*}
		\item Let $\sigma >0$. Assume that there exists a constant $ 0< K_0 \leq \frac{\nu}{8 C_L\sigma}$ such that $\|{\bf g}\|_{L^{\infty}(0,T;{\bf H}^1)} \leq K_0$, where ${C}_L$ is the constant from the inequality \eqref{interpolation}.
	\end{enumerate}

\smallskip
	\begin{remark}
		The assumption (G2) will be crucially used to obtain the exponential stability estimates in Lemma \ref{lemma_exp_discrete}. On the other hand, both (G1) and (G2) will be needed to derive the order of convergence in Theorem \ref{theorem_semi_chapter5} and Theorem \ref{theorem_fully_chapter5}.
	\end{remark}

	Next, we introduce the variational solution concept for \eqref{eq1.1} and refer the reader to \cite{Chow07,PZ1992} for a proof of its existence and uniqueness.

	\smallskip
	\begin{definition}\label{def2.1}
		Given $(\Omega,\mathcal{F}, \{\mathcal{F}_t\},\mathbb{P})$, let $W$ be an ${\mathbb R}$-valued Wiener process on it.
		Suppose ${\bf u}_0\in L^2(\Omega, {\mathbb V})$.
		An $\{\mathcal{F}_t\}$-adapted stochastic process  $\{{\bf u}(t) ; 0\leq t\leq T\}$ is called
		a variational solution of \eqref{eq1.1} if ${\bf u} \in  L^2\bigl(\Omega; C([0,T]; {\mathbb V}))
		\cap L^2\bigl(\Omega;0,T;{\bf H}^2_{per}(D)\bigr)$,
		and satisfies $\mathbb{P}$-a.s.~for all $t\in (0,T]$
		\begin{align}\label{equu2.8a}
			\bigl({\bf u}(t),  {\bf v} \bigr)& + \int_0^t  \nu\bigl(\nabla {\bf u}(s), \nabla {\bf v} \bigr)
			\,  ds + \int_0^t \big({\bf u}(s)\cdot\nabla{\bf u}(s),{\bf v}\big)\, ds
			\\\nonumber
			&=({\bf u}_0, {\bf v}) + {  \Bigl(\int_0^t {\bf g}(s)\, dW(s), {\bf v} \Bigr)}  \qquad\forall  \, {\bf v}\in {\mathbb V}\, .
		\end{align}
	\end{definition}

	We also recall some properties of the convection term. Define
	\begin{align*}
		b({\bf u}, {\bf v}, {\bf w}) := \bigl({\bf u}\cdot\nabla{\bf v},{\bf w}\bigr)\qquad \forall {\bf u}, {\bf v},{\bf w} \in {\bf H}^1_{per}(D).
	\end{align*}
	It is easy to check that $b({\bf u},\cdot,\cdot)$ is skew-symmetric for ${\bf u} \in \mathbb{V}$, that is,
	\begin{align}
		b({\bf u},{\bf v},{\bf w}) = - b({\bf u},{\bf w},{\bf v})\qquad\forall {\bf v}, {\bf w} \in {\bf H}^1_{per}
	\end{align}
	and
	\begin{align*}
		b({\bf u},{\bf v},{\bf v}) = 0 \qquad\forall {\bf v} \in {\bf H}^1_{per}(D).
	\end{align*}

	In addition, we also recall the following interpolation inequality for ${\bf H}^1$-functions:
	\begin{align}\label{interpolation}
		b({\bf u},{\bf v},{\bf w})
		\leq C_L\|{\bf u}\|^{1/2}_{{\bf L}^2}\|\nabla{\bf u}\|^{1/2}_{{\bf L}^2}\|\nabla{\bf v}\|_{{\bf L}^2} \|{\bf w}\|^{1/2}_{{\bf L}^2}\|\nabla{\bf w}\|^{1/2}_{{\bf L}^2}.
	\end{align}

	Definition \ref{def2.1} only defines the velocity $\mathbf{u}$ for \eqref{eq1.1},
	its associated pressure $p$ is subtle to define. We do so in the following theorem.

	\begin{theorem}\label{thm 2.2}
		Let $\{{\bf u}(t) ; 0\leq t\leq T\}$ be a variational solution of \eqref{eq1.1}. There exists a unique adapted process
	 $P\in {L^2\bigl(\Omega; L^2(0,T; H^1_{per}(D)/\mathbb{R})\bigr)}$ such that $(\mathbf{u}, P)$ satisfies
		$\mathbb{P}$-a.s.~for all $t\in (0,T]$
		\begin{subequations}\label{equu2.100}
			\begin{align}\label{equu2.10a}
				&\bigl({\bf u}(t),  {\bf v} \bigr) + \nu\int_0^t  \bigl(\nabla {\bf u}(s), \nabla {\bf v} \bigr) \, ds + \int_0^t \big({\bf u}(s)\cdot\nabla{\bf u}(s),{\bf v}\big)\, ds
				- \bigl(  \div \mathbf{v}, P(t) \bigr) \\
				&=({\bf u}_0, {\bf v}) + \int_0^t \big({\bf f}(s), {\bf v}\big) \, ds
				+  {\int_0^t  \bigl( {\bf B}\bigl({\bf u}(s)\bigr), {\bf v} \bigr)\, dW(s)}  \,\,\, \forall  \, {\bf v}\in {\bf H}^1_{per}(D)\, , \nonumber \\

				&\bigl(\div {\bf u}, q \bigr) =0 \qquad\forall \, q\in  L^2_{per}(D)/\mathbb{R} .  \label{equu2.10b}
			\end{align}
		\end{subequations}
	\end{theorem}
Since the proof is similar \cite[Theorem 1.3]{FPL2021}, so we omit it.

	\begin{lemma}[\cite{Breit,CP2012}]\label{stability_pdes}
		Let ${\bf u}$ be a solution defined in Definition \ref{def2.1}. Then we have
		\begin{enumerate}[{\rm (a)}]
			\item Assume that ${\bf u}_0 \in L^r\bigl(\Omega; \mathbb{V}\bigr)$ for some $r \geq 2$. Then there hold
			\begin{align*}
				\mathbb{E}\biggl[\Bigl(\sup_{0\leq t \leq T} \|\nabla{\bf u}(t)\|^2_{{\bf L}^2} + \int_0^T \nu\|\nabla^2{\bf u}(t)\|^2_{{\bf L}^2}\, dt\Bigr)^{\frac{r}{2}} \biggr] \leq C_r \mathbb{E}\Bigl[\|\nabla {\bf u}_0\|^r_{{\bf L}^2}\Bigr].
			\end{align*}
			\item Assume that ${\bf u}_0 \in L^r\bigl(\Omega; \mathbb{V}\cap{\bf H}^2(D)\bigr)\cap L^{5r}\bigl(\Omega,\mathbb{V}\bigr)$ for some $r \geq 2$. Then we have
			\begin{align*}
				\mathbb{E}\biggl[\Bigl(\sup_{0\leq t \leq T} \|\nabla^2{\bf u}(t)\|^2_{{\bf L}^2} + \int_0^T \nu\|\nabla^3{\bf u}(t)\|^2_{{\bf L}^2}\, dt\Bigr)^{\frac{r}{2}} \biggr] \leq C_r \mathbb{E}\biggl[\Bigl(\|{\bf u}_0\|^2_{{\bf H}^2} + \|\nabla{\bf u}_0\|^{10}_{{\bf L}^2}\Bigr)^{\frac{r}{2}}\biggr].
			\end{align*}
		\end{enumerate}
	\end{lemma}

	\subsection{H\"older continuity of the variational solution}
	We cite the following high moment H\"older continuity estimate for the variational solution whose proof
	can be found in \cite{Breit,CHP2012}.
	\begin{lemma}[\cite{Breit,CHP2012}]\label{lemma2.2}
		Suppose ${\bf u}_0$ satisfies the assumptions in Lemma \ref{stability_pdes} (b) for some $r >2$ and ${\bf g} \in L^2(\Omega;L^{\infty}(0,T;{\bf H}_{per}^2(D)))$. Then there exists a constant $C \equiv C(D_T, {\bf u}_0, {\bf g},r)>0$, such that the variational solution to problem \eqref{eq1.1} satisfies
		for $s,t \in [0,T]$
		\begin{align}\label{equu2.20a}
			\mathbb{E}\Bigl[\|{\bf u}(t) - {\bf u}(s)\|^{r}_{\mathbb{V}}\Bigr] \leq C|t-s|^{r\gamma}\qquad\forall \gamma \in \Bigl(0,\frac12\Bigr).
		\end{align}
	\end{lemma}

	\section{Semi-discretization in time}\label{section5.3}
	In this section, we consider the time discretization of \eqref{equu2.8a} that is based on the fully implicit Euler-Maruyama method.

	\subsection{Formulation and stability of the time discretization}\,
	The Euler-Maruyama time discretization of \eqref{equu2.8a} is given by the following algorithm.

	\textbf{Algorithm 1.} Let ${\bf u}^0 = {\bf u}_0$ be a given $\mathbb{V}$-valued random variable. Find $\bigl({\bf u}^{n+1}, p^{n+1}\bigr) \in L^2(\Omega;\mathbb{V} \times L^2_{per}(D)/\mathbb{R})$ such that $\mathbb{P}$-a.s.
	\begin{align}\label{equu3.1}
		\bigl({\bf u}^{n+1} - {\bf u}^{n},\pmb{\phi}\bigr) + \nu k\,\bigl(\nabla{\bf u}^{n+1},\nabla\pmb{\phi}\bigr) + k\,\bigl({\bf u}^{n+1}\cdot\nabla{\bf u}^{n+1},\pmb{\phi}\bigr)
		&- k\, \bigl(p^{n+1},\div \pmb{\phi}\bigr) \\\nonumber
		&= \bigl({\bf g}^n\Delta W_{n+1},\pmb{\phi}\bigr),\\
		\bigl(\div {\bf u}^{n+1},\psi\bigr) &=0
	\end{align}
	$\forall \pmb{\phi} \in {\bf H}^1_{per}(D)$ and $\psi \in L^2_{per}(D)$. Where ${\bf g}^n = {\bf g}(t_n)$ and $\Delta W_{n+1} = W(t_{n+1})- W(t_n)$.

	If $\pmb{\phi} \in \mathbb{V}$, then \eqref{equu3.1} reduces to
	\begin{align}\label{eq_reforms}
		\bigl({\bf u}^{n+1} - {\bf u}^{n},\pmb{\phi}\bigr) + \nu k\,\bigl(\nabla{\bf u}^{n+1},\nabla\pmb{\phi}\bigr) + k\,\bigl({\bf u}^{n+1}\cdot\nabla{\bf u}^{n+1},\pmb{\phi}\bigr)
		= \bigl({\bf g}^n\Delta W_{n+1},\pmb{\phi}\bigr).
	\end{align}

	We quote the following stability estimates for $\{{\bf u}^n\}$ given in \cite{BCP12}.

	\begin{lemma}\label{stability_means}
		Let ${\bf u}_0 \in L^{2^q}(\Omega;\mathbb{V})$ for an integer $1 \leq q <\infty$ be given, such that $\mathbb{E}\bigl[\|{\bf u}_0\|^{2^q}_{\mathbb{V}}\bigr] \leq C$. Then there exists a constant $C_{T,q} = C(T, q, {\bf u}_0)$ such that the following estimations hold:
		\begin{enumerate}[{\rm (i)}]
			\item $\displaystyle\mathbb{E}\biggl[\max_{1\leq n \leq M}\|{\bf u}^n\|^{2^q}_{\mathbb{V}} + \nu k\sum_{n=1}^M \|{\bf u}^n\|^{2^q-2}_{\mathbb{V}}\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\biggr] \leq C_{T,q}$.
			\item $\displaystyle\mathbb{E}\biggl[\sum_{n=1}^M \|{\bf u}^n - {\bf u}^{n-1}\|^2_{\mathbb{V}}\|{\bf u}^n\|^2_{\mathbb{V}}\biggr] \leq C_{T,2}$.
			\item $\displaystyle\mathbb{E}\Biggl[\biggl(\sum_{n=1}^M \|{\bf u}^{n} - {\bf u}^{n-1}\|^2_{\mathbb{V}}\biggr)^q + \biggl(\nu k\sum_{n=1}^M \|{\bf A}{\bf u}^{n}\|^2_{\mathbb{V}}\biggr)^q\Biggr] \leq C_{T,q}$.

		\end{enumerate}
	\end{lemma}

	Stability estimates for the pressure approximation $\{p^n\}$ can also be obtained accordingly using  the stochastic {\em inf-sup}
	estimate.

	\begin{lemma}\label{stability_pressures}
		Let $\{p^n\}_{n=1}^M$ be the pressure approximation from Algorithm 1. Under the assumptions of Lemma \ref{stability_means} and for $1\leq q < \infty$, we have
		\begin{align}\label{stability_pressure2}
			\mathbb{E}\biggl[\biggl(k\sum_{n=1}^M \|\nabla p^{n}\|^2_{{\bf L}^2}\biggr)^q\biggr] \leq \frac{C_{T,q}}{k^q}.
		\end{align}
	\end{lemma}

	Next, we present an exponential stability estimate for the velocity approximation from Algorithm 1. Such an estimate was first
	derived  in \cite[Theorem 8.3]{BM_space}. However, here we provide a new proof which eliminates the restriction that $M$ must  be sufficiently large required in \cite{BM_space}. This estimate will be used in the proofs of Theorem \ref{theorem_semi_chapter5} and Theorem \ref{theorem_fully_chapter5}.

	\begin{lemma}\label{lemma_exp_discrete}
		Assume that $\mathbb{E}[\exp(4\sigma \|\nabla{\bf u}_0\|^2_{{\bf L}^2})] \leq C$ for some $\sigma > 0$. Let $\{{\bf u}^n\}$ be generated by Algorithm 1 and assume that ${\bf g}$ satisfies $(G2)$. Then, there holds
		\begin{align}
			\mathbb{E}\Bigl[\exp\Bigl(&\sigma\max_{1\leq \ell \leq M}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} + \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\Bigr] \leq C_2,
		\end{align}
		where $C_2 = C_2(\sigma,\nu,T,\mathbb{E}[\exp(4\sigma \|\nabla{\bf u}_0\|^2_{{\bf L}^2})])$.
	\end{lemma}

 \begin{proof}
		Choosing $\pmb{\phi} = {\bf A}{\bf u}^{n+1}$ in \eqref{eq_reforms} and using integration by parts, the binomial formula $2(a,a-b) = \|a\|^2 - \|b\|^2 + \|a-b\|^2$, we obtain
		\begin{align}\label{equu3.13}
			\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2} &- \|\nabla{\bf u}^n\|^2_{{\bf L}^2} + \|\nabla({\bf u}^{n+1} - {\bf u}^n)\|^2_{{\bf L}^2} + 2\nu k \|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} \\\nonumber
			&=2\bigl(\nabla{\bf g}^{n}\Delta W_{n+1}, \nabla({\bf u}^{n+1}-{\bf u}^n)\bigr) +2\bigl(\nabla{\bf g}^{n}\Delta W_{n+1}, \nabla{\bf u}^n\bigr)\\\nonumber
			&\leq \|\nabla{\bf g}^n\Delta W_{n+1}\|^2_{{\bf L}^2} + \|\nabla({\bf u}^{n+1} - {\bf u}^n)\|^2_{{\bf L}^2} \\\nonumber
			&\qquad\qquad\qquad\qquad+ 2\bigl(\nabla{\bf g}^{n}\Delta W_{n+1}, \nabla{\bf u}^n\bigr).
		\end{align}
		Therefore,
		\begin{align}
			\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2} &- \|\nabla{\bf u}^n\|^2_{{\bf L}^2} + 2\nu k \|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} \\\nonumber
			&\leq \|\nabla{\bf g}^n\Delta W_{n+1}\|^2_{{\bf L}^2} + 2\bigl(\nabla{\bf g}^{n}\Delta W_{n+1}, \nabla{\bf u}^n\bigr).
		\end{align}
		Next, we lower the index by one and apply the summation $ \sum_{n=1}^{\ell}$ for any $1 \leq \ell \leq M$ to get
		\begin{align}\label{equu3.18}
			\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} +2\nu k\sum_{n=1}^{\ell} \|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} &\leq \|\nabla{\bf u}^0\|^2_{{\bf L}^2} + \sum_{n=1}^{\ell}\|\nabla{\bf g}^{n-1}\Delta W_{n}\|^2_{{\bf L}^2}\\\nonumber
			&+ 2\sum_{n=1}^{\ell}\bigl(\nabla{\bf g}^{n}\Delta W_{n+1}, \nabla{\bf u}^n\bigr).
		\end{align}

		Taking $\displaystyle \max_{1\leq \ell \leq M}$ on both sides of \eqref{equu3.18} leads to
		\begin{align}\label{equu3.19}
			\max_{1\leq \ell \leq M}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} +2\nu k\sum_{n=1}^{M} \|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} &\leq \|\nabla{\bf u}^0\|^2_{{\bf L}^2} + \sum_{n=1}^{M}\|\nabla{\bf g}^{n-1}\Delta W_n\|^2_{{\bf L}^2}\\\nonumber
			&+2\max_{1\leq \ell \leq M}\sum_{n=1}^{\ell}\bigl(\nabla{\bf g}^{n-1}\Delta W_n,\nabla{\bf u}^{n-1}\bigr)\\\nonumber
			&\leq 2\|\nabla{\bf u}^0\|^2_{{\bf L}^2} + 2K_0 \sum_{n=1}^{M}|\Delta W_n|^2\\\nonumber
			&+2\max_{2\leq \ell \leq M}\sum_{n=2}^{\ell}\bigl(\nabla{\bf g}^{n-1}\Delta W_n,\nabla{\bf u}^{n-1}\bigr).
		\end{align}
		For any $\sigma > 0$, denote $Z_{\ell} =2\sigma\sum_{n=2}^{\ell}\bigl(\nabla{\bf g}^{n-1}\Delta W_n,\nabla{\bf u}^{n-1}\bigr)$, for $2\leq \ell\leq M$. Define the following martingale by using a piece-wise linear function that defines as follow: for $s \in [t_n,t_{n+1})$, with $n =1,2,\cdots, M-1$, set $\underline{s} = t_n$, ${\bf u}^{\underline{s}} = {\bf u}^n$. Then $\displaystyle Z_{\ell} = \tilde{Z}_{t_{\ell}}$, where
		\begin{align}
			\tilde{Z}_{t} = 2\sigma\int_{t_1}^{t}\bigl(\nabla{\bf g}^{\underline{s}}\,d W(s),\nabla {\bf u}^{\underline{s}}\bigr)\qquad\forall t\in [t_1,T].
		\end{align}
		Then, $\tilde{Z}_t$ is a $\mathcal{F}_t$-martingale for $t \in [t_1,T]$ and its quadratic variation satisfies
		\begin{align}\label{equu3.21}
			\bigl<\tilde{Z}\bigr>_{t_{\ell}}  &\leq 4\sigma^2\int_{t_1}^{t_{\ell}} \|\nabla{\bf g}^{\underline{s}}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{\underline{s}}\|^2_{{\bf L}^2}\, ds\\\nonumber
			&\leq 4C_L\sigma^2K_0\, k\sum_{n=1}^{\ell}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2},
		\end{align}
		where $C_L$ is the constant from the inequality $\|\nabla{\bf u}^n\|^2_{{\bf L}^2} \leq C_L\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}$.

		We multiply \eqref{equu3.19} by $\sigma >0$ and then take the exponential on both sides to get
		\begin{align}\label{equu323}
			\exp\Bigl(&\sigma\max_{1\leq \ell \leq M}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} + \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\Bigr) \\\nonumber
			&\leq \exp\bigl(2\sigma\|\nabla{\bf u}^0\|^2_{{\bf L}^2}\bigr)\times \exp\Bigl(2K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\\\nonumber
			&\times\exp\Bigl(\max_{2\leq \ell \leq M}\bigl[Z_{\ell} -2\bigl<\tilde{Z}\bigr>_{t_{\ell}}\bigr]\Bigr)\\\nonumber
			&\times \exp\biggl(2\max_{2\leq \ell \leq M}\bigl<\tilde{Z}\bigr>_{t_{\ell}} - \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\biggr).
		\end{align}
		Moreover, from \eqref{equu3.21} we have
		\begin{align}\label{equu3.24}
			&2\max_{2\leq \ell \leq M}\bigl<\tilde{Z}\bigr>_{t_{\ell}} - \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2} \\\nonumber
			&\leq 8 C_L\sigma^2K_0\, k\sum_{n=1}^{\ell}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} - \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\\\nonumber
			&\leq \bigl(8 K_0C_L\sigma^2 - \nu\sigma\bigr) k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}.
		\end{align}
		Now, by using (G2), we have $8 K_0C_L\sigma^2 - \nu\sigma \leq 0$,
		which also implies that
		\begin{align}\label{equu3.26}
			\exp\biggl(2\max_{2\leq \ell \leq M}\bigl<\tilde{Z}\bigr>_{t_{\ell}} - \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\biggr)  \leq 1\qquad\mathbb{P}-a.s.
		\end{align}
		Next, taking the expectation on \eqref{equu323} and using \eqref{equu3.26} we obtain
		\begin{align}\label{equu3.27}
			\mathbb{E}\Bigl[\exp\Bigl(&\sigma\max_{1\leq \ell \leq M}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} + \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\Bigr] \\\nonumber
			&\leq \mathbb{E}\Bigl[\exp\bigl(2\sigma\|\nabla{\bf u}^0\|^2_{{\bf L}^2}\bigr)\times \exp\Bigl(2K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\\\nonumber
			&\times\exp\Bigl(\max_{2\leq \ell \leq M}\bigl[\tilde{Z}_{t_{\ell}} -2\bigl<\tilde{Z}\bigr>_{t_{\ell}}\bigr]\Bigr)\Bigr].
		\end{align}
		Using the H\"older inequality to separate the product on the right-hand side of \eqref{equu3.27} yields
		\begin{align}\label{equu3.28}
			\mathbb{E}\Bigl[\exp\Bigl(&\sigma\max_{1\leq \ell \leq M}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2} + \sigma\nu k\sum_{n=1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\Bigr]\\\nonumber
			&\leq \Bigl(\mathbb{E}\bigl[\exp\bigl(4\sigma\|\nabla{\bf u}^0\|^2_{{\bf L}^2}\bigr)\bigr]\Bigr)^{1/2}\\\nonumber
			&\qquad\times \biggl(\mathbb{E}\Bigl[\exp\Bigl(8K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\Bigr]\biggr)^{1/4}\\\nonumber
			&\qquad\times \Bigl(\mathbb{E}\Bigl[\exp\Bigl(\max_{2\leq \ell \leq M}\bigl[4\tilde{Z}_{t_{\ell}} -\frac{1}{2}\bigl<4\tilde{Z}\bigr>_{t_{\ell}}\bigr]\Bigr)\Bigr]\Bigr)^{1/4}.
		\end{align}
		Since $\Bigl\{\exp\Bigl(4\tilde{Z}_{t} -\frac{1}{2}\bigl<4\tilde{Z}\bigr>_{t}\Bigr)\Bigr\}$ is an exponential martingale on $[t_1,T]$, so it suffices to show that $\mathbb{E}\Bigl[\exp\Bigl(8K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\Bigr]$ is bounded. To the end, let $\beta = 8K_0\sigma$ and $X_n = |\Delta W_n|^2$,\,$ 1 \leq n \leq M$.
		Then, by using the independence of the increments of $W(t)$, we can rewrite
		\begin{align}\label{equation3.23}
			\mathbb{E}\Bigl[\exp\Bigl(8K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\Bigr] &= \mathbb{E}\Bigl[\exp\Bigl(\beta \sum_{n=1}^M X_n\Bigr)\Bigr] = \prod_{n=1}^M \mathbb{E}\bigl[\exp(\beta  X_n)\bigr].
		\end{align}
		In addition, by using \eqref{mean_wiener}, we have that for each $1\leq n\leq M$, $\mathbb{E}[X_n] = \mathbb{E}[|\Delta W_n|^2] = k$ and then by \eqref{a.s_wiener}, we also have $0 \leq X_n \leq C k^{\alpha}$\,\, a.s. \, $\forall \alpha \in (0,1)$.

		Next, we use Hoeffding's lemma to conclude that
		\begin{align}\label{equation3.24}
			\mathbb{E}\bigl[\exp(\beta X_n)\bigr] \leq \exp\bigl(\beta  \mathbb{E}[X_n] + \frac{C^2\beta^2}{8} k ^{2\alpha}\bigr) = \exp\bigl(\beta k + \frac{C^2\beta^2}{8} k ^{2\alpha}\bigr).
		\end{align}
		Substituting \eqref{equation3.24} into the right-hand side of \eqref{equation3.23}, we arrive at
		\begin{align*}
			\mathbb{E}\Bigl[\exp\Bigl(8K_0\sigma\sum_{n=1}^M|\Delta W_n|^2\Bigr)\Bigr] \leq  \exp\bigl(\beta M k + \frac{C^2\beta^2}{8} M k ^{2\alpha}\bigr)  =  \exp\bigl(\beta T + \frac{C^2\beta^2}{8} T k ^{2\alpha-1}\bigr).
		\end{align*}

		The proof is complete by choosing $\frac12 \leq \alpha < 1$.
	\end{proof}

	\subsection{High moment and pathwise error estimates for the velocity approximation}
	In this subsection, we present the first main result of this paper which establishes the optimal order
	high moment and sub-optimal order pathwise error estimates for the velocity approximation generated by Algorithm 1.

	\begin{theorem}\label{theorem_semi_chapter5}
		Let ${\bf u}$ be the variational solution to \eqref{equu2.100} and $\{{\bf u}^{n}\}_{n=1}^M$ be generated by  Algorithm 1. Assume that ${\bf u}_0 \in L^{q}(\Omega; \mathbb{V})$ and . Then there exists $C_q = C(T,q,{\bf u}_0,{\bf f})>0$ for any {integer} $2 \leq q < \infty$ and {real number} $0 < \gamma < \frac12$ such that
		\begin{align}\label{equu310}
			\bigl(\mathbb{E}\bigl[\max_{1\leq n\leq M}\|{\bf u}(t_n) - {\bf u}^n\|^{q}_{{\bf L}^2}\bigr]\bigr)^{\frac1q} + \biggl(\mathbb{E}\biggl[\Bigl(\nu k \sum_{n = 1}^M \|\nabla({\bf u}(t_n) &- {\bf u}^n)\|^2_{{\bf L}^2}\Bigr)^{q/2}\biggr]\biggr)^{\frac1q} \\\nonumber
			&\leq C_q\, k^{\frac{1}{2}-\gamma}.
		\end{align}
	\end{theorem}

	\smallskip

  The proof of Theorem \ref{theorem_semi_chapter5} follows immediately  from Lemma \ref{lemma3.4} and \ref{lemma3.5} below.

	\smallskip

	\begin{lemma}\label{lemma3.4}
		Let ${\bf e}^n := {\bf u}(t_n) - {\bf u}^n$. 	Under the assumptions of Theorem \ref{theorem_semi_chapter5}, there holds $\mathbb{P}$-a.s.
		\begin{align}
			\max_{1\leq \ell \leq M}\|{\bf e}^{\ell}\|^2_{{\bf L}^2} &+ {\nu k}\sum_{n=1}^{M} \|\nabla{\bf e}^{n}\|^2_{{\bf L}^2} \leq \bigl(A + B\bigr) \exp\biggl(k\sum_{n=1}^{M}D^n\biggr),
		\end{align}
		where
		\begin{align*}\label{equu3.288}
			A &:= 2\nu\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+ \frac{2C_L^2}{\nu} k \max_{1\leq \ell \leq M}\|{\bf e}^{\ell}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2},
		\end{align*}
		\begin{align*}
			B &:= 2\sum_{n=1}^{M}\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2\\\nonumber
			&\qquad+ 2\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\biggl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\biggr)\Bigr|,\\\nonumber
			D^n &:= \frac{2C_L^2}{\nu} \|\nabla{\bf u}^n\|^2_{{\bf L}^2}\qquad\forall 1\leq n\leq M.
		\end{align*}
	\end{lemma}

	\begin{proof}
		Subtracting \eqref{eq_reforms} from \eqref{equu2.8a} we obtain
		\begin{align}
			\nonumber	\bigl({\bf e}^{n+1}-{\bf e}^n,{\bf v}\bigr) + \nu k\bigl(\nabla{\bf e}^{n+1},\nabla{\bf v}\bigr) &= \nu \int_{t_n}^{t_{n+1}}\bigl(\nabla({\bf u}(s)-{\bf u}(t_{n+1})),\nabla{\bf v}\bigr)\, ds \\
			&+ \int_{t_n}^{t_{n+1}}\bigl({\bf u}^{n+1}\cdot\nabla{\bf u}^{n+1} - {\bf u}(s)\cdot\nabla{\bf u}(s),{\bf v}\bigr)\, ds\\\nonumber
			&+\biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s)-{\bf g}(t_n)\bigr)\,dW(s),{\bf v}\biggr).
		\end{align}
		Choosing ${\bf v} = {\bf e}^{n+1}$ yields
		\begin{align}\label{error_eq}
			\nonumber	\bigl({\bf e}^{n+1}-{\bf e}^n,{\bf e}^{n+1}\bigr) + \nu k\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2} &= \nu \int_{t_n}^{t_{n+1}}\bigl(\nabla({\bf u}(s)-{\bf u}(t_{n+1})),\nabla{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&+ \int_{t_n}^{t_{n+1}}\bigl({\bf u}^{n+1}\cdot\nabla{\bf u}^{n+1} - {\bf u}(s)\cdot\nabla{\bf u}(s),{\bf e}^{n+1}\bigr)\, ds\\
			&+\biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s)-{\bf g}(t_n)\bigr)\,dW(s),{\bf e}^{n+1}\biggr)\\\nonumber
			&:= {\tt I + II + III}.
		\end{align}
		The left-hand side of \eqref{error_eq} can be easily controlled by using the formula $2(a,a-b) = \|a\|^2 - \|b\|^2 + \|a-b\|^2$. The right-hand side of \eqref{error_eq} can be bounded as follows.

		Using Cauchy-Schwarz inequality, we obtain
		\begin{align}\label{equu311}
			{\tt I} &\leq \nu \int_{t_n}^{t_{n+1}} \|\nabla({\bf u}(t_{n+1}) - {\bf u}(s))\|_{{\bf L}^2}\|\nabla{\bf e}^{n+1}\|_{{\bf L}^2}\, ds\\\nonumber
			&\leq \nu\int_{t_n}^{t_{n+1}}\|\nabla({\bf u}(t_{n+1}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds + \frac{\nu k}{4}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2}.

		\end{align}
		We rewrite the second term as
		\begin{align}\label{equu312}
			{\tt II} &= \int_{t_n}^{t_{n+1}}\bigl(({\bf u}(s) - {\bf u}(t_{n+1}))\cdot\nabla {\bf u}(s),{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&\qquad+\int_{t_n}^{t_{n+1}} \bigl({\bf u}(t_{n+1})\cdot\nabla({\bf u}(s) - {\bf u}(t_{n+1})),{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&\qquad+ \int_{t_n}^{t_{n+1}}\bigl({\bf e}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&\qquad+ \int_{t_n}^{t_{n+1}}\bigl({\bf u}^{n+1}\cdot\nabla{\bf e}^{n+1},{\bf e}^{n+1}\bigr)\, ds\\\nonumber
			&=\int_{t_n}^{t_{n+1}}\bigl(({\bf u}(s) - {\bf u}(t_{n+1}))\cdot\nabla {\bf u}(s),{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&\qquad+\int_{t_n}^{t_{n+1}} \bigl({\bf u}(t_{n+1})\cdot\nabla({\bf u}(s) - {\bf u}(t_{n+1})),{\bf e}^{n+1}\bigr)\, ds \\\nonumber
			&\qquad+ \int_{t_n}^{t_{n+1}}\bigl({\bf e}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf e}^{n+1}\bigr)\, ds\\\nonumber
			&= {\tt II_1 + II_2 +II_3}.
		\end{align}
		Here we have used the facts that $\bigl({\bf u}^{n+1}\cdot\nabla{\bf e}^{n+1},{\bf e}^{n+1}\bigr) = 0$ and $\bigl({\bf e}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf e}^{n+1}\bigr) = \bigl({\bf e}^{n+1}\cdot\nabla{\bf u}(t_{n+1}),{\bf e}^{n+1}\bigr)$.
		In addition,
		\begin{align*}
			{\tt II_1} &\leq \int_{t_n}^{t_{n+1}} \|{\bf u}(s) - {\bf u}(t_{n+1})\|_{{\bf L}^4}\|{\bf u}(s)\|_{{\bf L}^4}\|\nabla{\bf e}^{n+1}\|_{{\bf L}^2}\, ds\\\nonumber
			&\leq \frac{2}{\nu}\int_{t_n}^{t_{n+1}} \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf L}^4}\|{\bf u}(s)\|^2_{{\bf L}^4} \, ds + \frac{\nu  k}{8}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\leq \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds + \frac{\nu k}{8}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2};

		\end{align*}
		\begin{align*}
			{\tt II_2} &\leq  \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n+1})\|^2_{{\bf L}^2} \, ds + \frac{\nu k}{8}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2};\\\nonumber
			{\tt II_3} &\leq k \|{\bf e}^{n+1}\|^2_{{\bf L}^4}\|\nabla{\bf u}(t_{n+1})\|_{{\bf L}^2} \\\nonumber
			&\leq C_L k\|\nabla{\bf e}^{n+1}\|_{{\bf L}^2}\|{\bf e}^{n+1}\|_{{\bf L}^2}\|\nabla{\bf u}(t_{n+1})\|_{{\bf L}^2}\\\nonumber
			&\leq \frac{C_L^2}{\nu}\,k \, \|{\bf e}^{n+1}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2} + \frac{\nu k}{4}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2}.
		\end{align*}
		Thus, \begin{align}\label{equu313}
			{\tt II} &\leq \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n+1})\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n+1})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+ \frac{\nu k}{2}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2} + \frac{C_L^2}{\nu}\, k\, \|{\bf e}^{n+1}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}.
		\end{align}
		Now, we estimate the noise term as below
		\begin{align}\label{equu314}
			{\tt III} &= \biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) -{\bf g}(t_n)\bigr)\, dW(s), {\bf e}^{n+1} - {\bf e}^n\biggr) \\\nonumber
			&\qquad+ \biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) -{\bf g}(t_n)\bigr)\, dW(s), {\bf e}^n\biggr)\\\nonumber
			&\leq \Bigl\|\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) - {\bf g}(t_n)\bigr)\, dW(s)\Bigr\|^2_{{\bf L}^2} + \frac14\|{\bf e}^{n+1} - {\bf e}^n\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+\biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) -{\bf g}(t_n)\bigr)\, dW(s), {\bf e}^n\biggr).
		\end{align}
		Substituting \eqref{equu311},\eqref{equu313} and \eqref{equu314} into \eqref{error_eq}, we obtain
		\begin{align}\label{equu315}
			\frac12 \bigl[\|{\bf e}^{n+1}\|^2_{{\bf L}^2} &- \|{\bf e}^n\|^2_{{\bf L}^2}\bigr] + \frac14 \|{\bf e}^{n+1} - {\bf e}^n\|^2_{{\bf L}^2} + \frac{\nu k}{4}\|\nabla{\bf e}^{n+1}\|^2_{{\bf L}^2} \\\nonumber
			&\leq \nu\int_{t_n}^{t_{n+1}}\|\nabla({\bf u}(t_{n+1}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{2}{\nu}\int_{t_n}^{t_{n+1}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n+1})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n+1})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+\frac{C_L^2}{\nu}\, k\, \|{\bf e}^{n+1}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \Bigl\|\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) - {\bf g}(t_n)\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2 \\\nonumber
			&\qquad+ \biggl(\int_{t_n}^{t_{n+1}}\bigl({\bf g}(s) -{\bf g}(t_n)\bigr)\, dW(s), {\bf e}^n\biggr).
		\end{align}
		Next, lowering the index $n$ in \eqref{equu315} by $1$ and applying the summation operator $ \sum_{n=1}^{\ell}$ to both sides of \eqref{equu315} for any $1\leq \ell \leq M$, we obtain
		\begin{align}\label{equu316}
			&\|{\bf e}^{\ell}\|^2_{{\bf L}^2} + \frac{1}{2}\sum_{n=1}^{\ell} \|{\bf e}^{n} - {\bf e}^{n-1}\|^2_{{\bf L}^2} + \frac{\nu k}{2}\sum_{n=1}^{\ell} \|\nabla{\bf e}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\leq 2\nu\sum_{n=1}^{\ell}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{\ell}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{\ell}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+\frac{2C_L^2}{\nu}\, k\,\sum_{n=1}^{\ell} \|{\bf e}^{n}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ 2\sum_{n=1}^{\ell}\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2\\\nonumber
			&\qquad+ 2\sum_{n=1}^{\ell}\biggl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\biggr)\\\nonumber
			&\leq 2\nu\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+\frac{2C_L^2}{\nu}\, k\,\sum_{n=1}^{\ell-1} \|{\bf e}^{n}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + \frac{2C_L^2k}{\nu}\max_{1\leq \ell \leq M}\|{\bf e}^{\ell}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ 2\sum_{n=1}^{M}\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2\\\nonumber
			&\qquad+ 2\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\biggl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\biggr)\Bigr|\\\nonumber
			&= A + B + k\sum_{n = 1}^{\ell-1} D^n\|{\bf e}^n\|^2_{{\bf L}^2},
		\end{align}
		where $A,B$ are constants in $\ell$ and are defined by
		\begin{align}\label{equu317}
			A &:= 2\nu\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
			&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n})\|^2_{{\bf L}^2} \, ds\\\nonumber
			&\qquad+ \frac{2C_L^2}{\nu} k \max_{1\leq \ell \leq M}\|{\bf e}^{\ell}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2}, \\\nonumber
			B &:= 2\sum_{n=1}^{M}\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2\\\nonumber
			&\qquad+ 2\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\biggl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\biggr)\Bigr|,\\\nonumber
			D^n &:= \frac{2C_L^2}{\nu} \|\nabla{\bf u}^n\|^2_{{\bf L}^2}\qquad\forall 1\leq n\leq M.
		\end{align}
		Finally, applying the discrete Gronwall inequality	to \eqref{equu316}, we obtain
		\begin{align}\label{equu318}
			\|{\bf e}^{\ell}\|^2_{{\bf L}^2} &+  \sum_{n=1}^{M} \|{\bf e}^{n} - {\bf e}^{n-1}\|^2_{{\bf L}^2} + {\nu k}\sum_{n=1}^{M} \|\nabla{\bf e}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\leq \bigl(A + B\bigr)\exp\biggl(k\sum_{n=1}^{\ell-1} D^n\biggr) \qquad\mathbb{P}-\mbox{a.s.}
		\end{align}
		The proof is complete by taking $\displaystyle\max_{1\leq \ell \leq M}$ on both sides.
	\end{proof}

	Next, we estimate $ \mathbb{E}\biggl[\bigl(A + B\bigr)^{\frac{q}{2}}\exp\biggl(\frac{kq}{2}\sum_{n=1}^{M} D^n\biggr)\biggr]$ in the following lemma to extract the convergent rates as stated in Theorem \ref{theorem_semi_chapter5}.

	\smallskip
	\begin{lemma}\label{lemma3.5}
		Assume that (G2) satisfies with $\sigma = \frac{2C_L^2}{\nu}$. Then, under the assumptions of Theorem \ref{theorem_semi_chapter5} and Lemma \ref{lemma3.4}, there holds
		\begin{align}
			\biggl(\mathbb{E}\biggl[\bigl(A + B\bigr)^{\frac{q}{2}}\exp\biggl(\frac{kq}{2}\sum_{n=1}^{M} D^n\biggr)\biggr]\biggr)^{\frac1q} \leq C_q\, k^{\frac12-\gamma}.
		\end{align}
		Here, $A,B,D^n$ are defined as in \eqref{equu317} and $0<\gamma<\frac12$. $C_q$ is the same constant as in Theorem \ref{theorem_semi_chapter5} which depends on $\nu, q, {\bf g}, D_T, {\bf u}_0$.
	\end{lemma}

	\begin{proof}
		By using the simple inequality $(a+b)^p \leq 2^{p-1}(a^p + b^p)$ for $1\leq p <\infty$ we have
		\begin{align}\label{equu3.40}
			\mathbb{E}\biggl[\bigl(A + B\bigr)^{\frac{q}{2}}\exp\biggl(\frac{kq}{2}\sum_{n=1}^{M} D^n\biggr)\biggr] &\leq C_q\mathbb{E}\biggl[A^{\frac{q}{2}}\exp\Bigl(\frac{kq}{2}\sum_{n = 1}^M D^n\Bigr)\biggr] \\\nonumber
			&\qquad+ C_q\mathbb{E}\biggl[B^{\frac{q}{2}}\exp\Bigl(\frac{kq}{2}\sum_{n = 1}^M D^n\Bigr)\biggr]\\\nonumber
			&\leq C_q\bigl(\mathbb{E}[A^q]\bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\exp\Bigl(qk\sum_{n = 1}^{M}D^n\Bigr)\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\quad+C_q\bigl(\mathbb{E}[B^q]\bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\exp\Bigl(qk\sum_{n = 1}^{M}D^n\Bigr)\Bigr]\Bigr)^{\frac12}.
		\end{align}
		Therefore, there are three factors in \eqref{equu3.40} which we need to estimate, namely, $ \mathbb{E}[A^q],\, \mathbb{E}[B^q]$ and $ \mathbb{E}\biggl[\exp\biggl(qk\sum_{n=1}^{M}D^n\biggr)\biggr]$.

		\begin{enumerate}[(a)]
			\item To estimate $\mathbb{E}[A^q]$, by \eqref{equu317} and the discrete Cauchy-Schwarz inequality, we have
			\begin{align*}
				\mathbb{E}[A^q] &= \mathbb{E}\biggl[\biggl(2\nu\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds \\\nonumber
				&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2} \, ds \\\nonumber
				&\qquad+ \frac{4}{\nu}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} C_L^2\, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1}\|\nabla{\bf u}(t_{n})\|^2_{{\bf L}^2} \, ds\\\nonumber
				&\qquad+ \frac{2C_L^2}{\nu} k \max_{1\leq \ell \leq M}\|{\bf e}^{\ell}\|^2_{{\bf L}^2}\|\nabla{\bf u}^{\ell}\|^2_{{\bf L}^2}	\biggr)^q\biggr]\\\nonumber
				&\leq C_q\mathbb{E}\biggl[\biggl(\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^2_{{\bf L}^2}\, ds\biggr)^q\biggr]\\\nonumber
				&\qquad+ C_q\mathbb{E}\biggl[\biggl(\sup_{s\in [0,T]}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} \, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1} \, ds\biggr)^q\biggr]\\\nonumber
				&\qquad+ C_qk^q\mathbb{E}\bigl[\max_{1\leq n \leq M}\|{\bf e}^n\|^q_{{\bf L}^2}\|\nabla{\bf u}^n\|^q_{{\bf L}^2}\bigr]\\\nonumber
				& = A_1 +A_2 +A_3.
			\end{align*}

			By using Cauchy-Schwarz inequality and \eqref{equu2.20a}, we have
			\begin{align*}
				A_1 &= C_q\mathbb{E}\biggl[\biggl(\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^{2}_{{\bf L}^2}\, ds\biggr)^q\biggr]\\\nonumber
				&\leq C_q\mathbb{E}\biggl[\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}}\|\nabla({\bf u}(t_{n}) - {\bf u}(s))\|^{2q}_{{\bf L}^2}\, ds\biggr]\\\nonumber
				&\leq C_qk^{2q(\frac12-\gamma)}.
			\end{align*}
			By using Cauchy-Schwarz inequality, \eqref{stability_means} and \eqref{equu2.20a}, we also obtain
			\begin{align*}
				A_2 &= C_q\mathbb{E}\biggl[\biggl(\sup_{s\in [0,T]}\|\nabla{\bf u}(s)\|^2_{{\bf L}^2}\sum_{n=1}^{M}\int_{t_{n-1}}^{t_{n}} \, \|{\bf u}(s) - {\bf u}(t_{n})\|^2_{{\bf H}^1} \, ds\biggr)^q\biggr]\\\nonumber
				&\leq C_q\mathbb{E}\biggl[\sup_{s\in [0,T]}\|\nabla{\bf u}(s)\|^{2q}_{{\bf L}^2}\sum_{n=1}^M\int_{t_{n-1}}^{t_n}\|{\bf u}(s) - {\bf u}(t_n)|^{2q}_{{\bf H}^1}\, ds\biggr]\\\nonumber
				&\leq C_q\bigl(\mathbb{E}\bigl[\sup_{s\in [0,T]}\|\nabla{\bf u}(s)\|^{4q}_{{\bf L}^2}\bigr]\bigr)^{\frac12}\biggl(\mathbb{E}\biggl[\sum_{n=1}^M\int_{t_{n-1}}^{t_n}\|{\bf u}(s) - {\bf u}(t_n)|^{4q}_{{\bf H}^1}\, ds\biggr]\biggr)^{\frac12}\\\nonumber
				&\leq C_qk^{2q(\frac12-\gamma)}.
			\end{align*}
			In addition, by using \eqref{stability_means}, $A_3 \leq C_q k^q$. In summary, $\mathbb{E}\bigl[A^q\bigr] \leq C_q k^{2q(\frac12-\gamma)}$.

			\smallskip
			\item To estimate $\mathbb{E}[B^q]$, by \eqref{equu317}, we have
			\begin{align*}
				\mathbb{E}[B^q] &= \mathbb{E}\biggl[\biggl(2\sum_{n=1}^{M}\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|_{{\bf L}^2}^2\\\nonumber
				&\qquad+ 2\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\Bigl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\Bigr)\Bigr|\biggr)^q\biggr]\\\nonumber
				&\leq C_q\mathbb{E}\biggl[\biggl(\sum_{n = 1}^M\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|^2_{{\bf L}^2}\biggr)^q\biggr]\\\nonumber
				&\qquad + C_q\mathbb{E}\biggl[\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\Bigl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\Bigr)\Bigr|^q\biggr]\\\nonumber
				&= B_1 + B_2.
			\end{align*}
			Next, by using the discrete H\"older inequality, \eqref{ito4} and (G1) we have

			\begin{align*}
				B_1 &= C_q\mathbb{E}\biggl[\biggl(\sum_{n = 1}^M\Bigl\|\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|^2_{{\bf L}^2}\biggr)^q\biggr]\\\nonumber
				&\leq C_q M^{q-1}\mathbb{E}\biggl[\sum_{n = 1}^M\Bigl\|\int_{t_{n-1}}^{t_n}\bigl({\bf g}(s) - {\bf g}(t_{n-1})\bigr)\, dW(s)\Bigr\|^{2q}_{{\bf L}^2}\biggr]\\\nonumber
				&\leq C_qM^{q-1}\mathbb{E}\biggl[\sum_{n = 1}^M\int_{t_{n-1}}^{t_n} \|{\bf g}(s) - {\bf g}(t_{n-1})\|^{2q}_{{\bf L}^2}\, ds\biggr]\\\nonumber&\leq C_q k^{q+1}.
			\end{align*}

			To estimate $B_2$, we use the Burkholder--Davis--Gundy inequality and (G1) to get
			\begin{align*}
				B_2 &= C_q\mathbb{E}\biggl[\max_{1\leq \ell \leq M}\Bigl|\sum_{n=1}^{\ell}\Bigl(\int_{t_{n-1}}^{t_{n}}\bigl({\bf g}(s) -{\bf g}(t_{n-1})\bigr)\, dW(s), {\bf e}^{n-1}\Bigr)\Bigr|^q\biggr]\\\nonumber
				&\leq C_q\mathbb{E}\biggl[\Bigl(\sum_{n = 1}^M\int_{t_{n-1}}^{t_n}\|{\bf g}(s) - {\bf g}(t_{n-1})\|^2_{{\bf L}^2}\|{\bf e}^{n-1}\|^{2}_{{\bf L}^2}\, ds\Bigr)^{q/2}\biggr]\\\nonumber
				&\leq C_qk^q \mathbb{E}\biggl[\Bigl(k\sum_{n = 1}^M\|{\bf e}^{n-1}\|^2_{{\bf L}^2}\Bigr)^{q/2}\biggr]\leq C_qk^q.
			\end{align*}
			In summary, $$\mathbb{E}\bigl[B^q\bigr] \leq C_q k^q. $$

			\smallskip
			\item To estimate $ \mathbb{E}\biggl[\exp\biggl(qk\sum_{n=1}^{M}D^n\biggr)\biggr]$ by \eqref{equu317} and the discrete Jensen inequality, we have
			\begin{align}\label{equu3.23}
			\nonumber 	\mathbb{E}\biggl[\exp\biggl(qk\sum_{n=1}^{M}D^n\biggr)\biggr] &\leq \mathbb{E}\biggl[\exp\biggl(\frac{2qC_L^2}{\nu}k\sum_{n=0}^{M}\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\biggr)\biggr]\\\nonumber
				&\leq \mathbb{E}\biggl[k\sum_{n=0}^M \exp\biggl(\frac{2qC_L^2}{\nu}\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\biggr)\biggr]\\
				&\leq k \sum_{n=0}^M \mathbb{E}\biggl[\max_{0\leq n \leq M}\exp\biggl(\frac{2qC_L^2}{\nu}\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\biggr)\biggr]
				 \leq CT.
			\end{align}
			Here, we have used Lemma \ref{lemma_exp_discrete} with $\sigma = \frac{2qC^2_L}{\nu}$ to get the last inequality of \eqref{equu3.23} .

		\end{enumerate}

		Finally, combining (a), (b), (c) and \eqref{equu3.40}, we get
		\begin{align}\label{equu325}
			\mathbb{E}\biggl[\bigl(A + B\bigr)^{\frac{q}{2}}\exp\biggl(\frac{kq}{2}\sum_{n=1}^{M} D^n\biggr)\biggr]\leq C_q k^{(\frac12-\gamma)q}.
		\end{align}
		The proof is complete.
	\end{proof}

	We are ready to state our first pathwise error estimate for the velocity approximation,
	such an estimate has not been obtained before in the literature.

	\begin{theorem}\label{thm_pathwise_semi}
		Assume that the assumptions of Theorem \ref{theorem_semi_chapter5} hold. Let $2 < q < \infty$ and $0<\gamma < \frac{1}{2}$ such that $\frac12 - \gamma - \frac{1}{q} >0$. Then, for $0 < \gamma_1 < \frac12 -  \gamma - \frac{1}{q}$, there exists a random variable $K = K(\omega)$ with $\mathbb{E}\bigl[|K|^{q}\bigr] <\infty$ such that there holds $\mathbb{P}$-a.s.
		\begin{align}\label{equu3.41}
			\max_{1\leq n \leq M}\|{\bf u}(t_n) - {\bf u}^n\|_{{\bf L}^2}  + \Bigl(\nu k\sum_{n=1}^M\|\nabla\bigl({\bf u}(t_n) - {\bf u}^n\bigr)\|^2_{{\bf L}^2}\Bigr)^{\frac12} \leq K\, k^{\gamma_1}.
		\end{align}
	\end{theorem}

	\begin{proof}
		\eqref{equu3.41} is an immediate consequence of Theorem \ref{theorem_semi_chapter5} and Kolmogorov Theorem \ref{kolmogorov}.
	\end{proof}

	\subsection{High moment and pathwise error estimates for the pressure approximation}
	We state and prove the following pathwise error estimate for the pressure approximation. Such a result has not been obtained before in the literature.
	\begin{theorem}\label{semi_pressure_error}
		Assume that the assumptions of Theorem \ref{theorem_semi_chapter5} hold. Then we have
		\begin{align}
			\biggl(\mathbb{E}\biggl[\Bigl\|P(t_m) - k\sum_{n=1}^m p^n\Bigr\|^q_{{\bf L}^2}\biggr]\biggr)^{\frac1q} \leq C_qk^{\frac12 - \gamma}.
		\end{align}
	\end{theorem}
	\begin{proof} The proof is based on the well-known inf-sup (LBB) condition which is quoted below. There exists a constant $\beta >0$ such that
		\begin{align}
			\beta \|\xi\|_{L^2} \leq \sup_{\pmb{\phi} \in {\bf H}^1_{per}(D)}\frac{\bigl(\xi,\div \pmb{\phi}\bigr)}{\| \nabla\pmb{\phi} \|_{{\bf L}^2}}\,\,\qquad\qquad\forall \xi \in L_{per}^2(D).
		\end{align}
		Set
		\begin{align*}
			{\bf e}^m := {\bf u}(t_m) - {\bf u}^m \mbox{ and }\mathcal{E}^m := P(t_m) - P^m, \mbox{ where } P^m := k\sum_{n = 1}^m p^{n}.
		\end{align*}
		Summing \eqref{equu3.1} from $1$ to $m$ we obtain
		\begin{align}\label{eq_summation_semi}
			\bigl({\bf u}^m,\pmb{\phi}\bigr) &+ \nu k\sum_{n = 1}^m\bigl(\nabla{\bf u}^n,\nabla\pmb{\phi}\bigr) + k\sum_{n = 1}^m\bigl({\bf u}^n\cdot\nabla{\bf u}^n,\pmb{\phi}\bigr) - \bigl(P^m,\div \pmb{\phi}\bigr) \\\nonumber
			&=\bigl({\bf u}^0,\pmb{\phi}\bigr) + \sum_{n = 1}^m\bigl({\bf g}^{n-1}\Delta W_n,\pmb{\phi}\bigr)\,\,\qquad\qquad\forall\pmb{\phi}\in {\bf H}^1_{per}(D).
		\end{align}
		Subtracting \eqref{eq_summation_semi} from \eqref{equu2.10a}, we obtain the following error equations:
		\begin{align}
			\bigl({\bf e}^m,\pmb{\phi}\bigr) &+ \Bigl(\nu k\sum_{n=1}^m \nabla{\bf e}^n,\nabla\pmb{\phi}\Bigr) + \nu k \sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\bigl(\nabla({\bf u}(s)-{\bf u}(t_n)),\nabla\pmb{\phi}\bigr)\, ds \\\nonumber
			&+ k\sum_{n = 1}^m \bigl({\bf u}(t_n)\cdot\nabla{\bf u}(t_n) - {\bf u}^n\cdot\nabla{\bf u}^n, \pmb{\phi}\bigr) \\\nonumber
			&+ \sum_{n=1}^m \int_{t_{n-1}}^{t_n}\bigl({\bf u}(s)\cdot\nabla{\bf u}(s) - {\bf u}(t_n)\cdot\nabla{\bf u}(t_n),\pmb{\phi}\bigr)\, ds - \bigl(\mathcal{E}^m,\div \pmb{\phi}\bigr)\\\nonumber
			&= \Bigl(\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}({\bf g}(s) - {\bf g}^{n-1})\, dW(s),\pmb{\phi}\Bigr).
		\end{align}
		Therefore,
		\begin{align}\label{equu3.50}
			\nonumber	\bigl(\mathcal{E}^m,\div \pmb{\phi}\bigr) &= \bigl({\bf e}^m,\pmb{\phi}\bigr) + \Bigl(\nu k\sum_{n=1}^m \nabla{\bf e}^n,\nabla\pmb{\phi}\Bigr) + \nu  \sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\bigl(\nabla({\bf u}(s)-{\bf u}(t_n)),\nabla\pmb{\phi}\bigr)\, ds \\\nonumber
			&\qquad+  k\sum_{n = 1}^m \bigl({\bf u}(t_n)\cdot\nabla{\bf u}(t_n) - {\bf u}^n\cdot\nabla{\bf u}^n, \pmb{\phi}\bigr) \\
			&\qquad+ \sum_{n=1}^m \int_{t_{n-1}}^{t_n}\bigl({\bf u}(s)\cdot\nabla{\bf u}(s) - {\bf u}(t_n)\cdot\nabla{\bf u}(t_n),\pmb{\phi}\bigr)\, ds\\\nonumber
			&\qquad+ \Bigl(\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}({\bf g}(t_{n-1}) - {\bf g}(s))\, dW(s),\pmb{\phi}\Bigr)\\\nonumber
			&:= {\tt I + II + III + IV + V +VI}.
		\end{align}
		Next, we use Cauchy-Schwarz and Poincar\'e inequalities to estimate each term on the right-hand side of \eqref{equu3.50}.
		\begin{align*}
			{\tt I + II} &\leq C\|{\bf e}^m\|_{{\bf L}^2}\|\nabla\pmb{\phi}\|_{{\bf L}^2} + \nu k \sum_{n = 1}^m\|\nabla{\bf e}^n\|_{{\bf L}^2}\|\nabla\pmb{\phi}\|_{{\bf L}^2}. \\\nonumber
			{\tt III} &\leq \nu \sum_{n=1}^m\int_{t_{n-1}}^{t_n}\|\nabla({\bf u}(s) - {\bf u}(t_n))\|_{{\bf L}^2}\, ds \|\nabla\pmb{\phi}\|_{{\bf L}^2}.\\\nonumber
			{\tt IV} &= k\sum_{n = 1}^m \bigl({\bf u}(t_n)\cdot\nabla{\bf u}(t_n) - {\bf u}^n\cdot\nabla{\bf u}^n, \pmb{\phi}\bigr)\\\nonumber
			&=k\sum_{n = 1}^m \bigl({\bf u}(t_n)\cdot\nabla{\bf u}(t_n) - {\bf u}^n\cdot\nabla{\bf u}(t_n), \pmb{\phi}\bigr) \\\nonumber
			&\qquad+ k\sum_{n = 1}^m \bigl({\bf u}^n\cdot\nabla{\bf u}(t_n) - {\bf u}^n\cdot\nabla{\bf u}^n, \pmb{\phi}\bigr) \\\nonumber
			&= k\sum_{n = 1}^m \bigl({\bf e}^n\cdot\nabla{\bf u}(t_n), \pmb{\phi}\bigr) + k\sum_{n = 1}^m \bigl({\bf u}^n\cdot\nabla{\bf e}^n, \pmb{\phi}\bigr)\\\nonumber
			&\leq k\sum_{n = 1}^m \|{\bf e}^n\|_{{\bf L}^4}\|{\bf u}(t_n)\|_{{\bf L}^4}\|\nabla\pmb{\phi}\|_{{\bf L}^2} + k\sum_{n = 1}^m\|{\bf u}^n\|_{{\bf L}^4}\|{\bf e}^n\|_{{\bf L}^4}\|\nabla\pmb{\phi}\|_{{\bf L}^2}. \\

			{\tt V} &= \sum_{n=1}^m \int_{t_{n-1}}^{t_n}\bigl({\bf u}(s)\cdot\nabla{\bf u}(s) - {\bf u}(t_n)\cdot\nabla{\bf u}(s),\pmb{\phi}\bigr)\, ds \\\nonumber
			&\qquad+ \sum_{n=1}^m \int_{t_{n-1}}^{t_n}\bigl({\bf u}(t_n)\cdot\nabla{\bf u}(s) - {\bf u}(t_n)\cdot\nabla{\bf u}(t_n),\pmb{\phi}\bigr)\, ds\\\nonumber
			&\leq \sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\|{\bf u}(s) - {\bf u}(t_n)\|_{{\bf L}^4}\|{\bf u}(s)\|_{{\bf L}^4}\|\nabla\pmb{\phi}\|_{{\bf L}^2}\, ds \\\nonumber
			&\qquad+ \sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\|{\bf u}(t_n)\|_{{\bf L}^4}\|{\bf u}(s) - {\bf u}(t_n)\|_{{\bf L}^4}\|\nabla\pmb{\phi}\|_{{\bf L}^2}\, ds.\\\nonumber
			{\tt VI} &\leq C\Bigl\|\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}({\bf g}(s) -{\bf g}(t_{n-1}))\, dW(s)\Bigr\|_{{\bf L}^2}\|\nabla\pmb{\phi}\|_{{\bf L}^2}.
		\end{align*}
		Substituting the above estimates into \eqref{equu3.50} and using the inf-sup condition yield
		\begin{align}\label{equu3.57}
			\nonumber	\beta\|\mathcal{E}^m\|_{L^2} &\leq \sup_{\pmb{\phi}\in {\bf H}^1_{per}(D)}\frac{\bigl(\mathcal{E}^m,\div \pmb{\phi}\bigr)}{\|\nabla\pmb{\phi}\|_{{\bf L}^2}}\\\nonumber
			&\leq C\|{\bf e}^m\|_{{\bf L}^2} + \nu k \sum_{n = 1}^m\|\nabla{\bf e}^n\|_{{\bf L}^2} +\nu \sum_{n=1}^m\int_{t_{n-1}}^{t_n}\|\nabla({\bf u}(s) - {\bf u}(t_n))\|_{{\bf L}^2}\, ds \\\nonumber
			&\qquad + k\sum_{n = 1}^m \|{\bf e}^n\|_{{\bf L}^4}(\|{\bf u}(t_n)\|_{{\bf L}^4} + \|{\bf u}^n\|_{{\bf L}^4})\\
			&\qquad + \sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\|{\bf u}(s) - {\bf u}(t_n)\|_{{\bf L}^4}\bigl(\|{\bf u}(s)\|_{{\bf L}^4} + \|{\bf u}(t_n)\|_{{\bf L}^4}\bigr)\, ds\\\nonumber
			&\qquad + C\Bigl\|\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}({\bf g}(s) -{\bf g}(t_{n-1}))\, dW(s)\Bigr\|_{{\bf L}^2}.
		\end{align}
		Taking the expectation on \eqref{equu3.57} and then using Theorem \ref{theorem_semi_chapter5}, Lemmas \ref{lemma2.2}--\ref{stability_means} and It\^o isometry, we obtain
		\begin{align*}
			\nonumber	\beta^q\mathbb{E}\bigl[\|\mathcal{E}^m\|^q_{L^2}\bigr] &\leq C_q\mathbb{E}[\|{\bf e}^m\|^q_{{\bf L}^2}] + C_q\mathbb{E}\Bigl[\Bigl(\nu k\sum_{n=1}^m\|\nabla{\bf e}^n\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigr] \\\nonumber
			&+ C_q \mathbb{E}\Bigl[\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\|\nabla({\bf u}(s) - {\bf u}(t_n))\|^q_{{\bf L}^2}\Bigr]\\\nonumber
			&+ C_q\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^m \|\nabla{\bf e}^n\|_{{\bf L}^2}\bigl(\|\nabla{\bf u}(t_n)\|_{{\bf L}^2} + \|\nabla{\bf u}^n\|_{{\bf L}^2}\bigr)\Bigr)^{q}\Bigr]\\\nonumber
			&+ C_q\mathbb{E}\Bigl[\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\|\nabla({\bf u}(s) - {\bf u}(t_n))\|^q_{{\bf L}^2}\bigl(\|\nabla{\bf u}(s)\|^q_{{\bf L}^2} + \|\nabla{\bf u}(t_n)\|^q_{{\bf L}^2}\bigr)\, ds\Bigr]\\
			&+ C_q\mathbb{E}\Bigl[\Bigl(\sum_{n = 1}^M\int_{t_{n-1}}^{t_n}\|{\bf g}(s) - {\bf g}(t_{n-1})\|^2_{{\bf L}^2}\, ds\Bigr)^{\frac{q}{2}}\Bigr]\\\nonumber
			&\leq C_q k^{(\frac12 - \gamma)q} \\\nonumber
			&+ C_q\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^m\|\nabla{\bf e}^n\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigl(k\sum_{n = 1}^m\bigl(\|\nabla{\bf u}(t_n)\|_{{\bf L}^2} + \|\nabla{\bf u}^n\|_{{\bf L}^2}\bigr)^2\Bigr)^{\frac{q}{2}}\Bigr]\\\nonumber
			&+ C_q\bigl(\mathbb{E}\bigl[\max_{s \in [0,T]}\|\nabla{\bf u}(s)\|^{2q}_{{\bf L}^2}\bigr]\bigr)^{\frac12}\sum_{n = 1}^m\int_{t_{n-1}}^{t_n}\bigl(\mathbb{E}\bigl[\|\nabla({\bf u}(s) - {\bf u}(t_n))\|^{2q}_{{\bf L}^2}\, \bigr]\bigr)^{1/2}ds\\\nonumber
			&\leq C_q k^{(\frac12 - \gamma)q} \\\nonumber
			&+ C_q\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^m\|\nabla{\bf e}^n\|^2_{{\bf L}^2}\Bigr)^{q}\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^m\bigl(\|\nabla{\bf u}(t_n)\|_{{\bf L}^2} + \|\nabla{\bf u}^n\|_{{\bf L}^2}\bigr)^2\Bigr)^{q}\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\leq C_q k^{(\frac12 - \gamma)q}.
		\end{align*}
		Finally, the proof is complete after  dividing  both sides by $\beta^q$.
	\end{proof}

	Next, we state a pathwise error estimate for the pressure approximation. To the best of our knowledge, this is the first pathwise convergence result for  the pressure approximation.
	\begin{theorem}\label{thm_semi_pathwise_pressure}
		Assume the assumptions of Theorem \ref{semi_pressure_error} hold. Let $2 < q < \infty$ and $0 < \gamma < \frac12$ such that $\frac12 - \gamma - \frac{1}{q} > 0$. Then, for $0 < \gamma_1 < \frac12 - \gamma - \frac{1}{q}$, there exists a random variable $K = K(\omega)$ with $\mathbb{E}\bigl[|K|^{q}\bigr] <\infty$ such that for all $1\leq \ell \leq M$, there holds $\mathbb{P}$-a.s.
		\begin{align}\label{equu3.52}
			\biggl\| P(t_{\ell}) - k\sum_{n=1}^{\ell} p^n\biggr\|_{L^2} \leq K\, k^{\gamma_1}.
		\end{align}
	\end{theorem}
	\begin{proof}
		The assertion follows immediately from an application of Theorem \ref{kolmogorov} using the high moment error estimate of Theorem \ref{semi_pressure_error}.
	\end{proof}

	\section{Mixed finite element discretization in space}\label{section5.4}
	 In this section, we consider the spatial approximation for the Algorithm 1 by using finite element method. Let $\mathcal{T}_h$ be a quasi-uniform triangular mesh of the given domain $D \subset \mathbb{R}^2$ with mesh size $h > 0$. We introduce the following Taylor-Hood mixed finite element spaces:
	\begin{align*}
		\mathbb{H}_h &:= \bigl\{{\bf v}_h \in C(\overline{D}) \cap {\bf H}^1_{per}(D);\,{\bf v}_h \in  [\mathcal{P}_2(K)]^2\, \,\forall K \in \mathcal{T}_h\bigr\},\\
		L_h &:= \bigl\{q_h \in L^2_{per}(D)/\mathbb{R};\, q_h \in \mathcal{P}_1(K)\,\,\forall K \in \mathcal{T}_h\bigr\},
	\end{align*}
	where $\mathcal{P}_i(K)$ denotes the space of all polynomials on $K$ of degree at most $i$. It's well-known \cite{Brezzi_Fortin91} that the above pair $(\mathbb{H}_h,L_h)$   satisfies the Ladyzh\v{e}nskaja-Babu\v{s}ka-Brezzi (LBB) condition,
	namely
	\begin{align}
		\sup_{\pmb{\phi}_h \in \mathbb{H}_h} \frac{\bigl(\div \pmb{\phi}_h,q_h\bigr)}{\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}} \geq C\|q_h\|_{L^2}\qquad\forall q_h\in L_h,
	\end{align}
	where the constant $C>0$ is independent of $h$ and $k$.

	Next, we define the space of weakly divergent-free velocity field as follows:
	\begin{align*}
		\mathbb{V}_h := \bigl\{\pmb{\phi}_h \in \mathbb{H}_h;\, \bigl(\div \pmb{\phi}_h, q_h\bigr) =0\qquad\forall q_h \in L_h\bigr\}.
	\end{align*}
	In general, $\mathbb{V}_h$ is not a subspace of $\mathbb{V}$.
	Let ${\bf Q}_h: {\bf L}^2 \rightarrow \mathbb{V}_h$ denote the $L^2$-orthogonal projection, which is defined by
	\begin{align}
		\bigl({\bf v} - {\bf Q}_h{\bf v}, \pmb{\phi}_h\bigr) = 0\qquad\forall \pmb{\phi}_h\in \mathbb{V}_h.
	\end{align}
	It is well-known \cite{Girault_Raviart86} that ${\bf Q}_h$ satisfies the following estimates:
	\begin{align}
		\|{\bf v} - {\bf Q}_h{\bf v}\|_{{\bf L}^2} + h\|\nabla({\bf v} - {\bf Q}_h{\bf v})\|_{{\bf L}^2} &\leq C h^2\|{\bf A}{\bf v}\|_{{\bf L}^2}\qquad\forall {\bf v}\in \mathbb{V}\cap{\bf H}^2(D),\\
		\|{\bf v} - {\bf Q}_h{\bf v}\|_{{\bf L}^2} &\leq Ch\|\nabla{\bf v}\|_{{\bf L}^2}\qquad\forall {\bf v}\in \mathbb{V}\cap{\bf H}^1(D).
	\end{align}
	Similarly, let $P_h: L^2_{per}(D) \rightarrow L_h$ be the scalar $L^2$-orthogonal projection, defined by
	\begin{align}
		\bigl(\psi - P_h\psi, q_h\bigr) = 0\qquad\forall q_h \in L_h,
	\end{align}
	then there also holds
	\begin{align}
		\|\psi - P_h\psi\|_{L^2} \leq Ch\|\nabla \psi\|_{{\bf L}^2}\qquad\forall \psi \in L^2_{per}(D)\cap H^1(D).
	\end{align}
	Next, we also introduce the well-known trilinear form \cite{Temam}
	\begin{align}
		\tilde{b}({\bf u},{\bf v},{\bf w}) = \bigl({\bf u}\cdot\nabla{\bf v},{\bf w}\bigr) + \frac12 \bigl([\div {\bf u}]{\bf v},{\bf w}\bigr)\qquad\forall {\bf u}, {\bf v}, {\bf w} \in {\bf H}^1(D),
	\end{align}
	which is anti-symmetric in the sense that
	\begin{align}
		\tilde{b}({\bf u},{\bf v},{\bf w}) = - \tilde{b}({\bf u},{\bf w},{\bf v}) \qquad\forall{\bf u},{\bf v},{\bf w} \in {\bf H}^1(D).
	\end{align}
	Therefore, $$\displaystyle \tilde{b}({\bf u}, \pmb{\phi},\pmb{\phi}) = 0\qquad\forall{\bf u},\pmb{\phi} \in {\bf H}^1(D).$$

	Our fully discrete mixed finite method is defined by the following algorithm.

	\textbf{Algorithm 2.} Let ${\bf u}_h^0$ be a given $\mathbb{H}_h$-valued random variable. Find $\displaystyle \bigl({\bf u}_h^{n+1},p_h^{n+1}\bigr) \in L^2(\Omega;\mathbb{H}_h\times L_h)$ such that there holds $\mathbb{P}$-a.s.
	\begin{align}
		\label{equu4.5}	&\bigl({\bf u}_h^{n+1} - {\bf u}_h^n,\pmb{\phi}_h\bigr) + \nu k \bigl(\nabla {\bf u}_h^{n+1},\nabla\pmb{\phi}_h\bigr) + k\bigl({\bf u}^{n+1}_h\cdot\nabla{\bf u}_h^{n+1},\pmb{\phi}_h\bigr) \\\nonumber
		&\qquad+ \frac{k}{2} \bigl([\div{\bf u}_h^{n+1}]{\bf u}_h^{n+1},\pmb{\phi}_h\bigr) - k\bigl(p_h^{n+1},\div \pmb{\phi}_h\bigr) = \bigl({\bf g}_h^n\Delta W_{n+1},\pmb{\phi}_h\bigr),\\
		&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\bigl(\div{\bf u}_h^n,\psi_h\bigr) =0,
	\end{align}
	for all $\pmb{\phi}_h \in \mathbb{H}_h$ and $\psi_h\in L_h$. ${\bf g}_h^n$ denotes an approximation of ${\bf g}^n$, which satisfies $\|{\bf g}^n - {\bf g}^n_h\|_{{\bf L}^2} \leq Ch^2$.

	\smallskip
	Next, we recall the following stability estimates for $\{{\bf u}_h^n\}$ from \cite[Lemma 3.1]{BCP12}.
	\begin{lemma}\label{stability_FEMs}
		Let $1 \leq q < \infty$ and ${\bf u}_h^0 \in L^{2^q}(\Omega;\mathbb{H}_h)$ such that $\displaystyle \mathbb{E}\bigl[\|{\bf u}_h^0\|^{2^q}_{{\bf L}^2}\bigr] \leq C$. Then, there exists a pair $\bigl\{({\bf u}_h^{n},p^{n}_h)\bigr\}_{n=1}^M \subset L^2(\Omega; \mathbb{H}_h\times L_h)$ that solves Algorithm 2 and satisfies
		\begin{enumerate}[{\rm (i)}]
			\item $\displaystyle \mathbb{E}\biggl[\max_{1\leq n \leq M}\|{\bf u}_h^n\|^{2^q}_{{\bf L}^2} + \nu k\sum_{n=1}^M \|{\bf u}_h^n\|^{2^{q}-2}_{{\bf L}^2}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\biggr] \leq C_{T,q},$
			\item $\displaystyle \mathbb{E}\biggl[\biggl(k\sum_{n=1}^M \|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\biggr)^{2^{q-1}}\,\biggr] \leq C_{T,q},$

		\end{enumerate}
		where $C_{T,q} = C_{T,q}(D_T,q,{\bf u}_h^0)$.
	\end{lemma}

	\subsection{High moment and pathwise error estimates for the fully discrete velocity approximation}
	In this subsection, we state and prove another  main result of this paper.

	\begin{theorem}\label{theorem_fully_chapter5}
		Let $2\leq q < \infty$ and ${\bf u}_h^0 = {\bf Q}_h {\bf u}_0$.  Assume that ${\bf u}_0 \in L^{q}(\Omega;\mathbb{V})$. Let $\{({\bf u}^n, p^n)\}$ and $\{{\bf u}^n_h,p^n_h\}$ be the velocity and pressure approximations generated by Algorithm 1 and 2, respectively. Then there holds
		\begin{align*}
			\bigl(\mathbb{E}\bigl[\max_{1\leq n \leq M}\|{\bf u}^n - {\bf u}_h^n\|^{q}_{{\bf L}^2}\bigr]\bigr)^{\frac1q} &+ \biggl(\mathbb{E}\biggl[\Bigl(\nu k \sum_{n=1}^M\|\nabla({\bf u}^n - {\bf u}_h^n)\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\biggr]\biggr)^{\frac1q}\\\nonumber
			&\leq C_q\biggl( k^{\frac12} + h + k^{-\frac{1}{2}}h\,\biggr),
		\end{align*}
		where $C_q = C(T,q,{\bf u}_0,{\bf f})>0$ is independent of $k$ and $h$.
	\end{theorem}

	\bigskip

	The proof of Theorem \ref{theorem_fully_chapter5} has several steps and will be accomplished through the following two lemmas.
	\begin{lemma}\label{lemma4.3} Under the assumptions of Theorem \ref{theorem_fully_chapter5}, there holds $\mathbb{P}$-a.s.
		\begin{align*}
			\max_{1\leq \ell \leq M}&\biggl[\|{\bf Q}_h({\bf u}^{\ell} - {\bf u}_h^{\ell})\|^2_{{\bf L}^2} + k\sum_{n=1}^{\ell}\|\nabla({\bf u}^{n} - {\bf u}_h^n)\|^2_{{\bf L}^2}\biggr] \leq \bigl(A_h + B_h\bigr)\exp\biggl(k\sum_{n=1}^{M} D^n_h\biggr),
		\end{align*}
		where
		\begin{align}\label{equu4.13}
			A_h &:= \|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{M}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{M}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}   \\\nonumber
			&\qquad+ Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\\\nonumber
			&\qquad +\frac{16C_L^2k}{\nu}\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{M}\|\nabla p^{n}\|^2_{{\bf L}^2},\\\nonumber
			B_h &:= \sum_{n=1}^{M} \|({\bf g}^{n-1}- {\bf g}_h^{n-1})\Delta W_{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+  \max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\Bigr|,\\\nonumber
			D_h^n &:= \frac{16C_L^2}{\nu}\|\nabla{\bf u}^n\|^2_{{\bf L}^2} + \frac{2\tilde{C}}{\nu}\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}.
		\end{align}
	\end{lemma}

	\begin{proof} Let ${\bf E}^n := {\bf u}^n - {\bf u}_h^n$ for all $0 \leq n \leq M$. Subtracting \eqref{equu3.1} from \eqref{equu4.5}, we obtain the following error equation
		\begin{align}\label{equu4.14}
			\bigl({\bf E}^{n+1} &- {\bf E}^{n},\pmb{\phi}_h\bigr) + \nu k \bigl(\nabla{\bf E}^{n+1},\nabla\pmb{\phi}_h\bigr) + k\,\tilde{b}({\bf u}^{n+1},{\bf u}^{n+1},\pmb{\phi}_h) \\\nonumber
			&\qquad- k \,\tilde{b}({\bf u}_h^{n+1},{\bf u}_h^{n+1},\pmb{\phi}_h)
			- k\bigl(p^{n+1} - p^{n+1}_h,\div\pmb{\phi}_h\bigr) \\\nonumber
			&\qquad= \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},\pmb{\phi}_h\bigr)\qquad\forall\pmb{\phi}_h\in\mathbb{H}_h.
		\end{align}
		Setting $\pmb{\phi}_h = {\bf Q}_h {\bf E}^{n+1} \in \mathbb{V}_h$ in \eqref{equu4.14}, using the orthogonality of the ${\bf L}^2$-projection and the fact that $\bigl(p^{n+1}_h,\div{\bf Q}_h{\bf E}^{n+1}\bigr) = 0$,  we obtain
		\begin{align}\label{equu4.15}
			\bigl({\bf Q}_h{\bf E}^{n+1} &- {\bf Q}_h{\bf E}^{n},{\bf Q}_h{\bf E}^{n+1}\bigr) + \nu k\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2}   \\\nonumber
			&= \nu k \bigl(\nabla{\bf E}^{n+1}, \nabla({\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1})\bigr) - k \,\tilde{b}({\bf u}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}) \\\nonumber
			&\qquad+ k \,\tilde{b}({\bf u}_h^{n+1},{\bf u}_h^{n+1},{\bf Q}_h{\bf E}^{n+1}) + 	 k\bigl(p^{n+1},\div{\bf Q}_h{\bf E}^{n+1}\bigr) \\\nonumber
			&\qquad+ \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr).
		\end{align}
		We use the identity $2(a,a-b) = \|a\|^2 - \|b\|^2 + \|a-b\|^2$ to rewrite the left-hand side of \eqref{equu4.15} as
		\begin{align}\label{equu4.16}
			\frac12&\bigl[\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^2} - \|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2} + \|{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\|^2_{{\bf L}^2}\bigr] + \nu k\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2}   \\\nonumber
			&= \nu k \bigl(\nabla{\bf E}^{n+1}, \nabla({\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1})\bigr) -\bigl[ k \,\tilde{b}({\bf u}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}) \\\nonumber
			&\qquad - k \,\tilde{b}({\bf u}_h^{n+1},{\bf u}_h^{n+1},{\bf Q}_h{\bf E}^{n+1})\bigr] + 	 k\bigl(p^{n+1},\div{\bf Q}_h{\bf E}^{n+1}\bigr) \\\nonumber
			&\qquad+ \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr)\\\nonumber
			&= {\tt I + II + III + IV}.
		\end{align}
		We note that the term {\tt II} needs a special treatment but other terms can be easily controlled as follows:
		\begin{align}
			{\tt I} \leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ck h^2\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}.
		\end{align}
		On noting that $\bigl(P_h p^{n+1}, \div {\bf Q}_h{\bf E}^{n+1}\bigr) = 0$, we have
		\begin{align}
			{\tt III} &= k\bigl(p^{n+1},\div{\bf Q}_h{\bf E}^{n+1}\bigr)\\\nonumber
			&=k\bigl(p^{n+1} - P_h p^{n+1}, \div{\bf Q}_h{\bf E}^{n+1}\bigr)\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla {\bf E}^{n+1}\|^2_{{\bf L}^2} + Ck h^2\|\nabla p^{n+1}\|^2_{{\bf L}^2}.
		\end{align}
		\begin{align}
			{\tt IV} &= \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h({\bf E}^{n+1}-{\bf E}^n)\bigr) + \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^n\bigr) \\\nonumber
			&\leq 2\|({\bf g}^n - {\bf g}^n_h)\Delta W_{n+1}\|^2_{{\bf L}^2} + \frac18\|{\bf Q}_h({\bf E}^{n+1} -{\bf E}^n)\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^n\bigr)
		\end{align}
		To bound the term {\tt II}, we first rewrite it as
		\begin{align}\label{eqII}
			{\tt II} &= -k\bigl[ \tilde{b}({\bf u}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}) - \tilde{b}({\bf u}^{n+1}_h, {\bf u}_h^{n+1},{\bf Q}_h{\bf E}^{n+1})\bigr]\\\nonumber
			&= -k\tilde{b}({\bf u}^{n+1},{\bf E}^{n+1},{\bf Q}_h{\bf E}^{n+1}) - k\tilde{b}({\bf E}^{n+1},{\bf u}_h^{n+1},{\bf Q}_h{\bf E}^{n+1})\\\nonumber
			&= -k\tilde{b}({\bf u}^{n+1},{\bf E}^{n+1},{\bf E}^{n+1}) + k\tilde{b}({\bf u}^{n+1}, {\bf E}^{n+1}, {\bf u}^{n+1} - {\bf Q}_h {\bf u}^{n+1}) \\\nonumber
			&\qquad+ k\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf Q}_h{\bf E}^{n+1}) - k\tilde{b}({\bf E}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1})\\\nonumber
			&= -k\tilde{b}({\bf u}^{n+1},{\bf E}^{n+1},{\bf E}^{n+1}) + k\tilde{b}({\bf u}^{n+1}, {\bf E}^{n+1}, {\bf u}^{n+1} - {\bf Q}_h {\bf u}^{n+1}) \\\nonumber
			&\qquad+ k\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf E}^{n+1}) - k\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}) \\\nonumber
			&\qquad\qquad- k\tilde{b}({\bf E}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1})\\\nonumber
			&= k\tilde{b}({\bf u}^{n+1}, {\bf E}^{n+1}, {\bf u}^{n+1} - {\bf Q}_h {\bf u}^{n+1}) - k\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}) \\\nonumber
			&\qquad\qquad- k\tilde{b}({\bf E}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1})\\\nonumber
			& := {\tt II_1 + II_2 + II_3}.
		\end{align}
		Here, in order to obtain the last line of \eqref{eqII}, we have used the facts that
		$\tilde{b}({\bf u}^{n+1},{\bf E}^{n+1},{\bf E}^{n+1}) = 0$ and $\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf E}^{n+1}) = 0$.
		Then we have
		\begin{align*}
			{\tt II_1} &= k\,\tilde{b}({\bf u}^{n+1},{\bf E}^{n+1}, {\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1})\\\nonumber
			&= k\bigl({\bf u}^{n+1}\cdot\nabla{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\bigr)+ \frac12 k\bigl([\div{\bf u}^{n+1}]{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\bigr)\\\nonumber
			&\leq k\|{\bf u}^{n+1}\|_{{\bf L}^4}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\|_{{\bf L}^4}\\\nonumber
			&\leq \frac{\nu k}{16} \|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ck\|{\bf u}^{n+1}\|^2_{{\bf L}^4}\|{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\|^2_{{\bf L}^4}\\\nonumber
			&\leq \frac{\nu k}{16} \|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}. \\

			{\tt II_2} &= -k\,\tilde{b}({\bf E}^{n+1},{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1})\\\nonumber
			&= -k\bigl({\bf E}^{n+1}\cdot\nabla{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\bigr) -\frac{k}{2}\bigl([\div {\bf E}^{n+1}]{\bf E}^{n+1},{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\bigr)\\\nonumber
			&\leq 2k\|{\bf E}^{n+1}\|_{{\bf L}^4}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\|_{{\bf L}^4}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ck\|{\bf E}^{n+1}\|^2_{{\bf L}^4}\|{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\|^2_{{\bf L}^4}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|{\bf E}^{n+1}\|^2_{{\bf L}^4}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|{\bf u}^{n+1}\|^2_{{\bf L}^4}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|{\bf u}_h^{n+1}\|^2_{{\bf L}^4}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} + Ckh^2\|{\bf u}_h^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}.
		\end{align*}
		Here, we have used the inverse inequality
		$\|{\bf u}_h^{n+1}\|_{{\bf L}^4} \leq Ch^{-\frac12}\|{\bf u}_h^{n+1}\|_{{\bf L}^2}$ to get the last inequality.
		\begin{align}
			{\tt II_3} &= -k\,\tilde{b}({\bf E}^{n+1},{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1})\\\nonumber
			&=-k\bigl({\bf E}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr) - \frac{k}{2}\bigl([\div{\bf E}^{n+1}]{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr)\\\nonumber
			&= -k\bigl({\bf E}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr) - \frac{k}{2}\bigl([\div{\bf E}^{n+1}]{\bf u}^{n+1},{\bf Q}_h({\bf E}^{n+1} - {\bf E}^{n})\bigr) \\\nonumber
			&\qquad+ \frac{k}{2}\bigl([\div{\bf E}^{n+1}]{\bf u}^{n+1},{\bf Q}_h{\bf E}^n\bigr) \\\nonumber
			&={\tt II_{3a} + II_{3b} + II_{3c}}.
		\end{align}
		We can bound each term on the right-hand side of ${\tt II_3}$ as follows.
		\begin{align*}
			{\tt II_{3a}} &= -k\bigl({\bf Q}_h{\bf E}^{n+1}\cdot\nabla{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr)- k\bigl(({\bf u}^{n+1}-{\bf Q}_h{\bf u}^{n+1})\cdot\nabla{\bf u}^{n+1},{\bf Q}_h{\bf E}^{n+1}\bigr)\\\nonumber
			&\leq k\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^4}\|\nabla{\bf u}^{n+1}\|_{{\bf L}^2} + k\|{\bf u}^{n+1}-{\bf Q}_h{\bf u}^{n+1}\|_{{\bf L}^4}\|\nabla{\bf u}^{n+1}\|_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n+1}\|_{{\bf L}^4}\\\nonumber
			&\leq C_Lk\|{\bf Q}_h{\bf E}^{n+1}\|_{{\bf L}^2}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|\nabla{\bf u}^{n+1}\|_{{\bf L}^2}\\\nonumber
			&\qquad+Ck\|{\bf u}^{n+1}-{\bf Q}_h{\bf u}^{n+1}\|_{{\bf L}^4}\|\nabla{\bf u}^{n+1}\|_{{\bf L}^2}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\\\nonumber
			&\leq \frac{8C_L^2}{\nu}k\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck\|{\bf u}^{n+1} - {\bf Q}_h{\bf u}^{n+1}\|^2_{{\bf L}^4}\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2} + \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\leq \frac{8C_L^2}{\nu} k\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} + \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2}.\\

			{\tt II_{3b}} &= -\frac{k}{2}\bigl([\div{\bf E}^{n+1}]{\bf u}^{n+1},{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\bigr)\\\nonumber
			&\leq \frac{k}{2}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|{\bf u}^{n+1}\|_{{\bf L}^4}\|{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\|_{{\bf L}^4}\\\nonumber
			&\leq C k \|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|\nabla{\bf u}^{n+1}\|_{{\bf L}^2}\|{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\|^{1/2}_{{\bf L}^2}\|\nabla ({\bf E}^{n+1} - {\bf E}^n)\|^{1/2}_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + Ck\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\|_{{\bf L}^2}\|\nabla ({\bf E}^{n+1} - {\bf E}^n)\|_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + \frac{1}{8}\|{\bf Q}_h({\bf E}^{n+1} - {\bf E}^n)\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck^2\|\nabla{\bf u}^{n+1}\|^4_{{\bf L}^2}\|\nabla({\bf E}^{n+1} - {\bf E}^n)\|^2_{{\bf L}^2}.\\\nonumber
			{\tt II_{3c}} &= \frac{k}{2}\bigl([\div{\bf E}^{n+1}]{\bf u}^{n+1},{\bf Q}_h{\bf E}^n\bigr)\\\nonumber
			&\leq \frac{k}{2}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|{\bf u}^{n+1}\cdot {\bf Q}_h{\bf E}^n\|_{{\bf L}^2}\\\nonumber
			&\leq \frac{k}{2}\|\nabla{\bf E}^{n+1}\|_{{\bf L}^2}\|{\bf u}^{n+1}\|_{{\bf L}^{\infty}}\|{\bf Q}_h{\bf E}^n\|_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + \frac{k}{\nu}\|{\bf u}^{n+1}\|^2_{{\bf L}^{\infty}}\|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2}\\\nonumber
			&\leq \frac{\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2} + \frac{\tilde{C}}{\nu} k\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2},
		\end{align*}
		where the last inequality was obtained by using the Sobolev embedding: $\|{\bf u}^{n+1}\|^2_{{\bf L}^{\infty}} \leq \tilde{C}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}$ for some constant $\tilde{C}>0$.

		Now, substituting the above estimates for the terms {\tt I}--{\tt IV} into the right-hand side of  the right-hand side of \eqref{equu4.16}, we obtain
		\begin{align}\label{equu4.31}
			\frac12\bigl[&\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^2} - \|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2}\bigr] +\frac14\|{\bf Q}_h({\bf E}^{n+1} - {\bf Q}_h{\bf E}^n)\|^2_{{\bf L}^2} \\\nonumber
			&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+ \frac{9\nu k}{16}\|\nabla{\bf E}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\leq Ckh^2\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} + Ckh^3\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ckh^2\|{\bf u}_h^{n+1}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2} + \frac{8C_L^2}{\nu} k\|\nabla{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ck^2\|\nabla{\bf u}^{n+1}\|^4_{{\bf L}^2}\|\nabla({\bf E}^{n+1} - {\bf E}^n)\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ \frac{\tilde{C}k}{\nu}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2} + Ck h^2\|\nabla p^{n+1}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \|({\bf g}^n - {\bf g}^n_h)\Delta W_{n+1}\|^2_{{\bf L}^2} + \bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^n\bigr).
		\end{align}

		Lowering the index $n$ in \eqref{equu4.31} by $1$ and taking the summation operator $\displaystyle \sum_{n=1}^{\ell}$ on both sides yield
		\begin{align}\label{equu4.32}
			\frac12&\|{\bf Q}_h{\bf E}^{\ell}\|^2_{{\bf L}^2} + \frac14\sum_{n=1}^{\ell}\|{\bf Q}_h({\bf E}^{n} - {\bf E}^{n-1})\|^2_{{\bf L}^2} + \frac{9\nu k}{16} \sum_{n=1}^{\ell}\|\nabla{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\leq \frac12\|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{\ell}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{\ell}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + \frac{8C_L^2}{\nu} k\sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck^2 \sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\bigl(\|\nabla({\bf u}_h^{n} - {\bf u}_h^{n-1})\|^2_{{\bf L}^2} + \|\nabla({\bf u}^n - {\bf u}^{n-1})\|^2_{{\bf L}^2}\bigr)\\\nonumber
			&\qquad+ \frac{\tilde{C}}{\nu} k  \sum_{n=1}^{\ell} \|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n-1}\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{\ell}\|\nabla p^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \sum_{n=1}^{\ell} \|({\bf g}^{n-1}- {\bf g}_h^{n-1})\Delta W_{n}\|_{{\bf L}^2}^2 + \sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr).
		\end{align}

		Bounding the sixth term on the right-hand side of \eqref{equu4.32} as follows
		\begin{align*}
			&Ck^2 \sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\bigl(\|\nabla({\bf u}_h^{n} - {\bf u}_h^{n-1})\|^2_{{\bf L}^2} + \|\nabla({\bf u}^n - {\bf u}^{n-1})\|^2_{{\bf L}^2}\bigr)\\\nonumber
			&\leq Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr),
		\end{align*}
		then, \eqref{equu4.32} becomes
		\begin{align}\label{equu4.322}
			\frac12&\|{\bf Q}_h{\bf E}^{\ell}\|^2_{{\bf L}^2} + \frac14\sum_{n=1}^{\ell}\|{\bf Q}_h({\bf E}^{n} - {\bf E}^{n-1})\|^2_{{\bf L}^2} + \frac{9\nu k}{16} \sum_{n=1}^{\ell}\|\nabla{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\leq \frac12\|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{\ell}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{\ell}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + \frac{8C_L^2}{\nu} k\sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\\\nonumber
			&\qquad+ \frac{\tilde{C}}{\nu} k  \sum_{n=1}^{\ell} \|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n-1}\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{\ell}\|\nabla p^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \sum_{n=1}^{\ell} \|({\bf g}^{n-1}- {\bf g}_h^{n-1})\Delta W_{n}\|_{{\bf L}^2}^2 + \sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\\\nonumber
			&\leq \frac12\|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{\ell}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{\ell}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{\ell}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+  k\sum_{n=1}^{\ell-1}\Bigl(\frac{8C_L^2}{\nu}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + \frac{\tilde{C}}{\nu}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\Bigr)\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\\\nonumber
			&\qquad+ \frac{8C_L^2k}{\nu} \max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{\ell}\|\nabla p^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \sum_{n=1}^{\ell} \|({\bf g}^{n-1}- {\bf g}_h^{n-1})\Delta W_{n}\|_{{\bf L}^2}^2 + \sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\\\nonumber
			&\leq \frac12\|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{M}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{M}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+  k\sum_{n=1}^{\ell-1}\Bigl(\frac{8C_L^2}{\nu}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + \frac{\tilde{C}}{\nu}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\Bigr)\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+ Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\\\nonumber
			&\qquad+ \frac{8C_L^2k}{\nu} \max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^n\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{M}\|\nabla p^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ \sum_{n=1}^{M} \|({\bf g}^{n-1}- {\bf g}_h^{n-1})\Delta W_{n}\|_{{\bf L}^2}^2 \\\nonumber
			&\qquad+ \max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\Bigr|,
		\end{align}
		which is equivalent to (see the notations defined in \eqref{equu4.16})
		\begin{align}\label{equu4.34}
			&\|{\bf Q}_h{\bf E}^{\ell}\|^2_{{\bf L}^2} + \sum_{n=1}^{\ell}\|{\bf Q}_h({\bf E}^{n} - {\bf E}^{n-1})\|^2_{{\bf L}^2} + \nu k \sum_{n=1}^{\ell}\|\nabla{\bf E}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\leq A_h +B_h + k\sum_{n=1}^{\ell-1}D^n_h\max_{1 \leq \ell \leq n}\|{\bf Q}_h{\bf E}^{\ell}\|^2_{{\bf L}^2}\\\nonumber
			& \leq \bigl(A_h + B_h\bigr) \exp\Bigl( k\sum_{n=1}^{m-1}D^n_h\Bigr).
		\end{align}

		The proof is completed by replacing $m$ by $M$ in \eqref{equu4.34}.
	\end{proof}

	\smallskip
	\begin{lemma}\label{lemma4.4} Under the assumptions of Theorem \ref{theorem_fully_chapter5} and Lemma \ref{lemma4.3}, there holds
		\begin{align}
			\biggl(\mathbb{E}\biggl[\bigl(A_h +& B_h\bigr)^{\frac{q}{2}}\exp\Bigl(\frac{qk}{2}\sum_{n=1}^M D^n_h\Bigr)\biggr]\biggr)^{\frac1q}  \leq C_q \biggl(k^{\frac12} + h + k^{-\frac{1}{2}}h\,\biggr).
		\end{align}
	\end{lemma}
	\begin{proof}
		First, by using Cauchy-Schwarz inequality, we have:
		\begin{align}\label{equu4.35}
			\mathbb{E}\biggl[\bigl(A_h +& B_h\bigr)^{\frac{q}{2}}\exp\Bigl( \frac{qk}{2}\sum_{n=1}^M D_h^n\Bigr)\biggr] \\\nonumber
			&\leq C_q\bigl(\mathbb{E}\bigl[A_h^q\bigr]\bigr)^{\frac{1}{2}}\biggl(\mathbb{E}\Bigl[\exp\Bigl(qk\sum_{n = 1}^MD^n_h\Bigr)\Bigr]\biggr)^{\frac12} \\\nonumber
			&\qquad+ C_q\bigl(\mathbb{E}\bigl[B_h^q\bigr]\bigr)^{\frac{1}{2}}\biggl(\mathbb{E}\Bigl[\exp\Bigl(qk\sum_{n = 1}^MD^n_h\Bigr)\Bigr]\biggr)^{\frac12}.
		\end{align}
		By using Lemma \ref{lemma_exp_discrete}, we can easily control the factor $\biggl(\mathbb{E}\Bigl[\exp\Bigl(q k\sum_{n=1}^M D^n_h\Bigr)\Bigr]\biggr)^{\frac12}$. In fact,
		\begin{align*}
			\mathbb{E}\Bigl[\exp\Bigl(qk\sum_{n=1}^M D^n_h\Bigr)\Bigr] &= \mathbb{E}\Bigl[\exp\Bigl( k\sum_{n=1}^M \frac{32qC_L^2}{\nu}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\frac{4q\tilde{C}}{\nu}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\Bigr)\Bigr] \\\nonumber
			&\leq \mathbb{E}\Bigl[\exp\Bigl( \frac{32qC_L^2T}{\nu}\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\frac{4q\tilde{C}}{\nu}\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\Bigr)\Bigr] \\\nonumber
			&\leq \mathbb{E}\Bigl[\exp\Bigl( \sigma\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2} + \sigma k\sum_{n = 1}^M\|{\bf A}{\bf u}^{n+1}\|^2_{{\bf L}^2}\Bigr)\Bigr] \leq C_2,
		\end{align*}
		where $\sigma_q = \max\bigl\{\frac{32qC_L^2 T}{\nu}, \frac{4q\tilde{C}}{\nu}\bigr\}$.

		Therefore, it remains to bound $\mathbb{E}[A_h^q]$ and $\mathbb{E}[B_h^q]$. By the assumptions on ${\bf g}$, we have
		\begin{align}\label{equu4.36}
			\mathbb{E}[B_h^q] &= \mathbb{E}\Bigl[\Bigl(\sum_{n=1}^{M} \|({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n}\|^2_{{\bf L}^2} \\\nonumber
			&\qquad+  \max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^n\bigr)\Bigr|\Bigr)^q\Bigr]\\\nonumber
			&\leq C_q\mathbb{E}\Bigl[\Bigl(\sum_{n = 1}^M\|{\bf g}^{n-1}-{\bf g}_h^{n-1}\|^2_{{\bf L}^2}|\Delta W_n|^2\Bigr)^q\Bigr]\\\nonumber
			&\qquad+ C_q\mathbb{E}\Big[\max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n+1},{\bf Q}_h{\bf E}^n\bigr)\Bigr|^q\Big]\\\nonumber
			&\leq C_qh^{2q}\mathbb{E}\Bigl[\Bigl(\sum_{n = 1}^M|\Delta W_n|^2\Bigr)^q\Bigr] \\\nonumber
			&\qquad+ C_q\mathbb{E}\Big[\max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\Bigr|^q\Big]\\\nonumber
			&= B_{h,1} + B_{h,2}.
		\end{align}
		It follows from the discrete H\"older inequality and \eqref{mean_wiener} that
		\begin{align}
			B_{h,1} &= C_qh^{2q}\mathbb{E}\Bigl[\Bigl(\sum_{n = 1}^M|\Delta W_n|^2\Bigr)^q\Bigr]\\\nonumber
			&\leq C_qh^{2q}\, M^{q-1}\sum_{n = 1}^M\mathbb{E}\bigl[|\Delta W_n|^{2q}\bigr]\leq C_qh^{2q} M^{q} k^q \leq C_qh^{2q}.
		\end{align}
		To estimate $B_{h,2}$, we use the Burkholder--Davis--Gundy inequality and the assumptions on ${\bf g}$ to obtain
		\begin{align}
			B_{h,2} &= C_q\mathbb{E}\Big[\max_{1\leq \ell \leq M}\Bigl|\sum_{n = 1}^{\ell}\bigl(({\bf g}^n - {\bf g}_h^n)\Delta W_{n},{\bf Q}_h{\bf E}^{n-1}\bigr)\Bigr|^q\Big]\\\nonumber
			&\leq C_q\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|{\bf g}^{n-1} - {\bf g}^{n-1}_h\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n-1}\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}} \Bigr]\\\nonumber
			&\leq C_q h^{2q} \mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|{\bf E}^{n-1}\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigr] \leq C_qh^{2q}.
		\end{align}
		In addition, using the stability estimates from Lemmas \ref{stability_means}, \ref{stability_FEMs} and \ref{stability_pressures} and the assumption $\mathbb{E}[\|{\bf Q}_h{\bf E}^0\|^q_{{\bf L}^2}] \leq Ch^q$, we obtain
		\begin{align}\label{equu4.37}
			\mathbb{E}[A_h^q] &= \mathbb{E}\bigg[\biggl(\|{\bf Q}_h{\bf E}^0\|^2_{{\bf L}^2} + Ch^2  k\sum_{n=1}^{M}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2} + Ch^3 k\sum_{n=1}^{M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\\\nonumber
			&\qquad+ Ch^2 k\sum_{n=1}^{M}\|{\bf u}_h^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}   \\\nonumber
			&\qquad+ Ck\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^4_{{\bf L}^2}\Bigl(k\sum_{n = 1}^{M}\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2} + k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)\\\nonumber
			&\qquad +\frac{16C_L^2k}{\nu}\max_{1 \leq n \leq M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf Q}_h{\bf E}^{n}\|^2_{{\bf L}^2} + Ch^2 k \sum_{n=1}^{M}\|\nabla p^{n}\|^2_{{\bf L}^2}\biggr)^q\bigg]\\\nonumber
			&\leq C_qh^{2q} + C_qh^{2q}\mathbb{E}\Bigl[\Bigl(k\sum_{n=1}^{M}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\Bigr)^q\Bigr]\\\nonumber
			&\qquad+ Ch^{3q}\mathbb{E}\Bigl[\Bigl(k\sum_{n=1}^{M}\|\nabla{\bf u}^{n}\|^2_{{\bf L}^2}\|{\bf A}{\bf u}^{n}\|^{2}_{{\bf L}^2}\Bigr)^q\Bigr] \\\nonumber
			&\qquad+C_qh^{2q}\mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|{\bf u}_h^{n}\|^{2q}_{{\bf L}^2}\Bigl(k\sum_{n=1}^{M}\|{\bf A}{\bf u}^{n}\|^2_{{\bf L}^2}\Bigr)^{q}\Bigr]\\\nonumber
			&\qquad+C_qk^q\mathbb{E}\biggl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{4q}_{{\bf L}^2}\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}^n_h\|^2_{{\bf L}^2}\Bigr)^q\biggr] \\\nonumber
			&\qquad+  C_qk^q\mathbb{E}\biggl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{4q}_{{\bf L}^2}\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)^{q}\biggr]\\\nonumber
			&\qquad + C_qk^q\mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{2q}_{{\bf L}^2}\|{\bf E}^n\|^{2q}_{{\bf L}^2}\Bigr] + C_qh^{2q}\mathbb{E}\Bigl[\Bigl(k \sum_{n=1}^{M}\|\nabla p^{n}\|^2_{{\bf L}^2}\Bigr)^{q}\Bigr]\\\nonumber
			&\leq C_qh^{2q} + C_qh^{2q} \Bigl(\mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|{\bf u}_h^n\|^{4q}_{{\bf L}^2}\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|{\bf A}{\bf u}^n\|^2_{{\bf L}^2}\Bigr)^{2q}\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\qquad+ C_qk^q \Bigl(\mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{8q}_{{\bf L}^2}\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\Bigr)^{2q}\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\qquad+ C_qk^q \Bigl(\mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{8q}_{{\bf L}^2}\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}^n\|^2_{{\bf L}^2}\Bigr)^{2q}\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\qquad+ C_qk^q \mathbb{E}\Bigl[\max_{1 \leq n \leq M}\|\nabla{\bf u}^n\|^{2q}_{{\bf L}^2}\bigl(\|{\bf u}^n\|^{2q}_{{\bf L}^2} + \|{\bf u}_h^n\|^{2q}_{{\bf L}^2}\bigr)\Bigr]\\\nonumber
			&\qquad+ C_qh^{2q}\mathbb{E}\Bigl[\Bigl(k \sum_{n=1}^{M}\|\nabla p^{n}\|^2_{{\bf L}^2}\Bigr)^q\Bigr]\leq C_qh^{2q} + C_qk^q + C_q k^{-q}h^{2q}.
		\end{align}
		We note that the last term $C_qk^{-q}h^{2q}$ is a consequence of allowing general noises which clearly impact on the stability estimate of the time-semidiscrete pressure approximation as seen from Lemma \ref{stability_pressures}.

		In summary, substituting the estimates for $A_h$ and $B_h$ into \eqref{equu4.35}, we obtain
		\begin{align}
			\mathbb{E}\biggl[\bigl(A_h +& B_h\bigr)^{\frac{q}{2}}\exp\Bigl(\frac{qk}{2}\sum_{n=1}^M D_h^n\Bigr)\biggr] \leq C_q\Bigl(k^{\frac{q}{2}} + h^q + k^{-\frac{q}{2}}h^q\Bigr).
		\end{align}
		The proof is complete.
	\end{proof}

	We conclude this section by stating a pathwise error estimate for the velocity approximation generated by Algorithm 2 which is a direct corollary of the Kolmogorov Theorem (cf. Theorem \ref{kolmogorov}) and
	the high moment error estimates of Theorem \ref{theorem_fully_chapter5}.
	\begin{theorem}\label{theorem_pathwise_fully_chapter5}
		Assume that the assumptions of Theorem \ref{theorem_fully_chapter5} hold. Let $2\leq q <\infty$ and $0 < \gamma_2 < 1 - \frac{1}{q}$. Then, there exists a random variable $K = K(\omega)$ with $\mathbb{E}\bigl[|K|^q\bigr] < \infty$ such that there holds $\mathbb{P}$-a.s.
		\begin{align}
			\max_{1\leq n \leq M}\|{\bf u}^n - {\bf u}^n_h\|_{{\bf L}^2} + \Bigl(\nu k \sum_{n=1}^M\|\nabla({\bf u}^n &- {\bf u}_h^n)\|^2_{{\bf L}^2}\Bigr)^{1/2}  \\\nonumber
			&\leq K\bigl(k^{\frac{\gamma_2}{2}} + h^{\gamma_2} + k^{-\frac{\gamma_2}{2}}h^{\gamma_2}\bigr).
		\end{align}
	\end{theorem}

	\subsection{High moment and pathwise error estimates for the fully discrete pressure approximation}
	In this subsection, we establish high moment and pathwise error estimates for the pressure approximation generated by Algorithm 2.
	\begin{theorem}\label{fully_pressure_error}
		Let $2\leq q <\infty$, under the assumptions of Theorem \ref{theorem_fully_chapter5}, there holds
		\begin{align}\label{equu4.24}
			\biggl(\mathbb{E}\biggl[\biggl\|k\sum_{n=1}^M\big(p^n - p^n_h\big)\biggr\|^{q}_{L^2}\biggr]\biggr)^{\frac1q} \leq C_q\bigl(k^{\frac12} + h + k^{-\frac{1}{2}}h\bigr),
		\end{align}
		where $C_q = C(\beta_1, q,D_T, {\bf f},{\bf u}_0)>0$ and independent of $k$ and $h$.
	\end{theorem}

\begin{proof}
		The proof of this theorem is similar to that of Theorem \ref{semi_pressure_error}.
		We recall the discrete inf-sup (LBB) condition: There exists a $\beta_1 > 0$ such that
		\begin{align}
			\beta_1\|\xi_h\|_{L^2} \leq \sup_{\pmb{\phi}_h\in \mathbb{H}_h}\frac{\bigl(\xi_h,\div \pmb{\phi}_h\bigr)}{\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}}\,\,\,\qquad\forall \xi_h\in L_h.
		\end{align}
		Set ${\bf E}^m := {\bf u}^m - {\bf u}^m_h$ and $E_p^m := k\sum_{n = 1}^m (p^n - p^n_h)$.

		Summing \eqref{equu4.14} from $1$ to $M$, we obtain
		\begin{align}
			\bigl({\bf E}^M - {\bf E}^0,\pmb{\phi}_h\bigr) &+ \nu k\sum_{n = 1}^M\bigl(\nabla{\bf E}^n,\nabla\pmb{\phi}_h\bigr) + k\sum_{n=1}^M\bigl({\bf u}^n\cdot\nabla{\bf u}^n - {\bf u}_h^n\cdot\nabla{\bf u}_h^n, \pmb{\phi}_h\bigr) \\\nonumber
			&-\frac{k}{2}\sum_{n=1}^M\bigl([\div{\bf u}^n_h]{\bf u}_h^n,\pmb{\phi}_h\bigr) -\bigl(E_p^M,\div \pmb{\phi}_h\bigr) \\\nonumber
			&= \Bigl(\sum_{n=1}^M({\bf g}^{n-1} - {\bf g}^{n-1}_h)\Delta W_{n},\pmb{\phi}_h\Bigr)\,\,\qquad\forall\pmb{\phi}_h\in \mathbb{H}_h.
		\end{align}
		Thus,
		\begin{align}
			\bigl(E_p^M,\div \pmb{\phi}_h\bigr) &= \bigl({\bf E}^M - {\bf E}^0,\pmb{\phi}_h\bigr) + \nu k\sum_{n = 1}^M\bigl(\nabla{\bf E}^n,\nabla\pmb{\phi}_h\bigr) \\\nonumber
			&\qquad+ k\sum_{n=1}^M\bigl({\bf u}^n\cdot\nabla{\bf u}^n - {\bf u}_h^n\cdot\nabla{\bf u}_h^n, \pmb{\phi}_h\bigr) \\\nonumber
			&\qquad-\frac{k}{2}\sum_{n=1}^M\bigl([\div{\bf u}^n_h]{\bf u}_h^n,\pmb{\phi}_h\bigr) - \Bigl(\sum_{n=1}^M({\bf g}^{n-1} - {\bf g}^{n-1}_h)\Delta W_{n},\pmb{\phi}_h\Bigr)\\\nonumber
			&= \bigl({\bf E}^M - {\bf E}^0,\pmb{\phi}_h\bigr) + \nu k\sum_{n = 1}^M\bigl(\nabla{\bf E}^n,\nabla\pmb{\phi}_h\bigr) \\\nonumber
			&\qquad+ k\sum_{n=1}^M\bigl({\bf E}^n\cdot\nabla{\bf u}^n + {\bf u}_h^n\cdot\nabla{\bf E}^n, \pmb{\phi}_h\bigr) \\\nonumber
			&\qquad+ \frac{k}{2}\sum_{n = 1}^M \bigl([\div {\bf E}^n]{\bf u}_h^n,\pmb{\phi}_h\bigr) - \Bigl(\sum_{n=1}^M({\bf g}^{n-1} - {\bf g}^{n-1}_h)\Delta W_{n},\pmb{\phi}_h\Bigr)\\\nonumber
			&\leq C(\|{\bf E}^M\|_{{\bf L}^2} + \|{\bf E}^0\|_{{\bf L}^2})\|\nabla\pmb{\phi}_h\|_{{\bf L}^2} + \nu k \sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}\\\nonumber
			&\qquad+ k\sum_{n=1}^M(\|{\bf E}^n\|_{{\bf L}^4}\|\nabla{\bf u}^n\|_{{\bf L}^2} + \|{\bf u}_h^n\|_{{\bf L}^4}\|\nabla{\bf E}^n\|_{{\bf L}^2})\|\pmb{\phi}_h\|_{{\bf L}^4}\\\nonumber
			&\qquad+ \frac{k}{2}\sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\|{\bf u}^n_h\|_{{\bf L}^4}\|\pmb{\phi}_h\|_{{\bf L}^4}\\\nonumber
			&\qquad+ C\Bigl\|\sum_{n = 1}^M({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_n\Bigr\|_{{\bf L}^2}\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}\\\nonumber
			&\leq C(\|{\bf E}^M\|_{{\bf L}^2} + \|{\bf E}^0\|_{{\bf L}^2})\|\nabla\pmb{\phi}_h\|_{{\bf L}^2} + \nu k \sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}\\\nonumber
			&\qquad+Ck\sum_{n=1}^M(\|\nabla{\bf E}^n\|_{{\bf L}^2}\|\nabla{\bf u}^n\|_{{\bf L}^2} \\\nonumber
			&\qquad+ \|\nabla{\bf u}_h^n\|_{{\bf L}^2}\|\nabla{\bf E}^n\|_{{\bf L}^2})\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}\\\nonumber
			&\qquad+ \frac{Ck}{2}\sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\|\nabla{\bf u}^n_h\|_{{\bf L}^2}\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}\\\nonumber
			&\qquad+ C\Bigl\|\sum_{n = 1}^M({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_n\Bigr\|_{{\bf L}^2}\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}.
		\end{align}
		Using the discrete inf-sup condition, we obtain
		\begin{align}\label{equu4.45}
			\beta_1\|E_p^M\|_{L^2} &\leq \sup_{\pmb{\phi}_h\in \mathbb{H}_h}\frac{\bigl(E_p^M,\div \pmb{\phi}_h\bigr)}{\|\nabla\pmb{\phi}_h\|_{{\bf L}^2}}\\\nonumber
			&\leq C(\|{\bf E}^M\|_{{\bf L}^2} + \|{\bf E}^0\|_{{\bf L}^2}) + \nu k \sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\\\nonumber
			&\qquad+Ck\sum_{n=1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}(\|\nabla{\bf u}^n\|_{{\bf L}^2} + \|\nabla{\bf u}_h^n\|_{{\bf L}^2})\\\nonumber
			&\qquad+ \frac{Ck}{2}\sum_{n = 1}^M\|\nabla{\bf E}^n\|_{{\bf L}^2}\|\nabla{\bf u}^n_h\|_{{\bf L}^2} \\\nonumber
			&\qquad+ C\Bigl\|\sum_{n = 1}^M({\bf g}^{n-1} - {\bf g}_h^{n-1})\Delta W_n\Bigr\|_{{\bf L}^2}.
		\end{align}
		Taking the $qth$-power, expectation, and then using Theorem \ref{theorem_fully_chapter5} and the stability estimates of Lemmas \ref{stability_means} and \ref{stability_FEMs}, we obtain
		\begin{align}
			\nonumber		\beta_1\mathbb{E}\bigl[\|E_p^M\|^q_{L^2}\bigr] &\leq C_q\Bigl(k^{\frac{q}{2}} + h^q + k^{-\frac{q}{2}}h^{q}\Bigr) \\\nonumber
			&\qquad+ C_q\mathbb{E}\biggl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf E}^n\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigl(k\sum_{n = 1}^M\bigl(\|\nabla{\bf u}^n\|^2_{{\bf L}^2} + \|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\bigr)\Bigr)^{\frac{q}{2}}\biggr]\\\nonumber
			&\qquad+ C_q\mathbb{E}\biggl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf E}^n\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\biggr]\\\nonumber
			&\qquad+ C_q\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|{\bf g}^{n-1} - {\bf g}^{n-1}_h\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\Bigr]\\\nonumber
			&\leq C_q\Bigl(k^{\frac{q}{2}} + h^q + k^{-\frac{q}{2}}h^{q}\Bigr) \\
			&\qquad+ C_q\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf E}^n\|^2_{{\bf L}^2}\Bigr)^q\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\bigl(\|\nabla{\bf u}^n\|^2_{{\bf L}^2} + \|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\bigr)\Bigr)^q\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\qquad+ C_q\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf E}^n\|^2_{{\bf L}^2}\Bigr)^q\Bigr]\Bigr)^{\frac12}\Bigl(\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|\nabla{\bf u}_h^n\|^2_{{\bf L}^2}\Bigr)^q\Bigr]\Bigr)^{\frac12}\\\nonumber
			&\qquad+ C_q\mathbb{E}\Bigl[\Bigl(k\sum_{n = 1}^M\|{\bf g}^{n-1} - {\bf g}^{n-1}_h\|^2_{{\bf L}^2}\Bigr)^{\frac12}\Bigr]\\\nonumber
			&\leq C_q\Bigl(k^{\frac{q}{2}} + h^q + k^{-\frac{q}{2}}h^{q}\Bigr).
		\end{align}
		Finally, the proof is complete after dividing both sides by $\beta_1$.
	\end{proof}

	An immediate consequence of the above high moment error estimates is the following pathwise error estimate for the pressure approximation $\{p_h^n\}$, it follows immediately from an application of Theorem \ref{kolmogorov}.

	\begin{theorem}\label{cor_4.5_chapter5}
		Assume that the assumptions of Theorem \ref{fully_pressure_error} hold. Let $2\leq q <\infty$ and $0 < \gamma_2 < 1 - \frac{1}{q}$. Then, there exists a random variable $K = K(\omega)$ with $\mathbb{E}\bigl[|K|^q\bigr] < \infty$ such that there holds $\mathbb{P}$-a.s.
		\begin{align}
			\Bigl\|k\sum_{n=1}^M\bigl(p^n - p^n_h\bigr)\Bigr\|_{L^2} \leq K\biggl(k^{\frac{\gamma_2}{2}} + h^{\gamma_2} + k^{-\frac{\gamma_2}{2}}h^{\gamma_2}\biggr).
		\end{align}
	\end{theorem}

	 \smallskip
	We conclude this section by stating the global error estimates for our fully discrete numerical
	solution generated by Algorithm 2, which follows from combining the above temporal and spatial error estimates.

	\begin{theorem} \label{thm4.8}
		Let $2 \leq q <\infty$ and $0< \gamma <\frac12$, under the assumptions of Theorem \ref{theorem_semi_chapter5} and Theorem \ref{theorem_fully_chapter5}, there exists a constant $C_q = C(D_T, {\bf u}_0, q, {\bf f})>0$ such that
		\begin{align}\label{eq4.47}
			\bigl(\mathbb{E}\bigl[\max_{1\leq n \leq M}\|{\bf u}(t_n) - {\bf u}^n_h\|^q_{{\bf L}^2}\bigr]\bigr)^{\frac1q} &+ \biggl(\mathbb{E}\biggl[\Bigl(\nu k \sum_{n=1}^{M}\|\nabla\bigl({\bf u}(t_n) - {\bf u}^n_h\|\bigr)\|^2_{{\bf L}^2}\Bigr)^{\frac{q}{2}}\biggr]\biggr)^{\frac1q} \\\nonumber
			&\leq C_q\Bigl( k^{\frac12 - \gamma} + h + k^{-\frac{1}{2}}h\Bigr).
		\end{align}
		In addition, let $2<q<\infty$ and $0<\gamma < \frac12$ such that $\frac12 - \gamma - \frac{1}{q}>0$. Then, for any $0<\gamma_1 < \frac12 - \gamma - \frac{1}{q}$ and $0<\gamma_2 < 1 - \frac{1}{q}$, there exists a random variable $K$ with $\mathbb{E}[|K|^q] < \infty$ such that there holds $\mathbb{P}-a.s.$
		\begin{align}\label{equu5.2}
			\max_{1\leq n \leq M}\|{\bf u}(t_n) - {\bf u}^n_h\|_{{\bf L}^2} &+ \Bigl(\nu k \sum_{n=1}^{M}\|\nabla\bigl({\bf u}(t_n) - {\bf u}^n_h\bigr)\|^2_{{\bf L}^2}\Bigr)^{\frac12} \\\nonumber
			&\leq K\Bigl(k^{\gamma_1} + h^{\gamma_2} + k^{-\frac{\gamma_2}{2}}h^{\gamma_2}\Bigr).
		\end{align}
	\end{theorem}

	\begin{theorem} \label{thm4.9}
		Let $2 \leq q <\infty$ and $0< \gamma <\frac12$. Under the assumptions of Theorem \ref{semi_pressure_error} and Theorem \ref{fully_pressure_error}, there exists a constant $C_q = C(D_T, {\bf u}_0, q, {\bf f},\beta_0,\beta_1)>0$ such that for $1 \leq \ell \leq M$
		\begin{align}\label{equu4.33}
			\biggl(\mathbb{E}\biggl[\biggl\|P(t_{\ell}) - k\sum_{n=1}^{\ell} p^n_h\biggr\|^q_{L^2}\biggr]\biggr)^{\frac1q} \leq C_q\Bigl( k^{\frac12 - \gamma} + h + k^{-\frac{1}{2}}h\Bigr),
		\end{align}
		In addition, let $2<q<\infty$ and $0<\gamma < \frac12$ such that $\frac12 - \gamma - \frac{1}{q}>0$. Then, for any $0<\gamma_1 < \frac12 - \gamma - \frac{1}{q}$ and $0<\gamma_2 < 1 - \frac{1}{q}$, there exists a random variable $K$ with $\mathbb{E}[|K|^q] < \infty$ such that there holds $\mathbb{P}-a.s.$
		\begin{align}\label{eq5.104}
			\biggl\|P(t_{\ell}) - k\sum_{n=1}^{\ell} p^n_h\biggr\|_{L^2} \leq K\Bigl(k^{\gamma_1} + h^{\gamma_2} + k^{-\frac{\gamma_2}{2}}h^{\gamma_2}\Bigr).
		\end{align}
	\end{theorem}

	\section{Numerical experiments} \label{sec-5}
		In this section, we present three numerical experiments to validate our error estimates in \eqref{eq4.47}--\eqref{eq5.104}.
		In all our experiments,  we set $D = (0,1)^2\subset \mathbb{R}^2$, $T =1, \nu =1$, $g \equiv 10$, the body force is ${\bf f} = (f_1,f_2)$ with
		\begin{align*}
			f_1(x,y) &= \pi\cos(t)\sin(2\pi y)\sin(\pi x)\sin(\pi x)
			- 2\pi^3\sin(t)\sin(2\pi y)(2\cos(2\pi x)-1) \\\nonumber
			&\qquad- \pi\sin(t)\sin(\pi x)\sin(\pi y),\\\nonumber
			f_2(x,y) &= -\pi\cos(t)\sin(2\pi x)\sin(\pi y)\sin(\pi y)
			- 2\pi^3\sin(t)\sin(2\pi x)(1-2\cos(2\pi y)) \\\nonumber
			&\qquad+ \pi\sin(t)\cos(\pi x)\cos(\pi y).
		\end{align*}
		We choose ${\bf W}$ in \eqref{eq1.1} to be a ${\mathbb R}^J$-valued Wiener process, with  increment
		\begin{align}\label{used_noise}
			{\bf W}^J(t_{n+1},{\bf x}) - {\bf W}^J(t_n,{\bf x}) = \sqrt{k_0} \sum\limits_{j_1=1}^J\sum\limits_{j_1=1}^J \sqrt{\lambda_{j_1,j_2}}{\bf e}_{j_1,j_2}({\bf x})\xi_{j_1,j_2}^n\,,
		\end{align}
		where ${\bf x}=(x,y) \in D, \, \xi_{j_1,j_2}^n \sim N(0,1)$, $\lambda_{j_1,j_2} = \frac{1}{(j_1+j_2)^2}$, and
		\begin{align}
			{\bf e}_{j_1,j_2}({\bf x}) = \bigl(\sin(j_1\pi x)\sin(j_2 \pi y),\sin(j_1\pi x)\sin(j_2 \pi y)\bigr)\, .
		\end{align}
		In addition, we use the Taylor-Hood mixed finite element method for the spatial discretization and the homogeneous Dirichlet boundary condition is imposed on ${\bf u}$.
		Moreover, we choose $J = 4$ and $N_p =300$ as number of samples in our Monte Carlo simulations to compute the expectations in each test.
		We implement Algorithm 2 and evaluate the errors between the computed velocity (resp. pressure) and the reference velocity (resp. pressure) which are generated by Algorithm 2 on a very fine mesh with mesh sizes $k_0 = \frac{1}{1024}, h =\frac{1}{40}$. Furthermore, to evaluate errors in strong norms, we use the following numerical
		integration formulas: for any $0 \leq m \leq M$ and any integer $2\leq q <\infty$,
		\begin{align*}
			\pmb{\mathcal{E}}_{{\bf u},q}^m &:=\biggl(\mathbb{E}\Bigl[\|{\bf u}(t_m) -{\bf u}_h^m(k)\|_{{\bf L}^2}^q \Bigr]\biggr)^{\frac1q} \approx
			\biggl(\frac{1}{N_p}\sum_{\ell=1}^{N_p} \bigl\|{\bf u }_h^m(k_0,\omega_\ell)
			-{\bf u}_h^m(k,\omega_\ell) \bigr\|_{{\bf L}^2}^q \biggr)^{\frac1q}\, ,\\\nonumber
			{\mathcal{E}}_{P,q}^m &:=\biggl(\mathbb{E}\Bigl[\Big\|P(t_m) -k\sum_{n=1}^m p_h^n(k)\Big\|_{{\bf L}^2}^q \Bigr]\biggr)^{\frac1q}
			\\\nonumber
			&\,\,\approx \biggl(\frac{1}{N_p}\sum_{\ell=1}^{N_p}\Bigl\|k\sum_{n=1}^{\frac{t_m}{k_0}} p_{h}^{n}(k_0,\omega_\ell) -k\sum_{n=1}^{\frac{t_m}{k}}p_h^n(k,\omega_\ell)\Bigr\|_{{\bf L}^2}^q\Bigr) \biggr)^{\frac1q}\, .
		\end{align*}

		In addition, in order to compute the quantities of stochastic interest, we want to generate as many samples as possible in the Monte Carlo simulations. However, this approach is very expensive if Newton's method is used at each time step to solve the nonlinear system in Algorithm 2. Therefore, we introduce a cheaper fixed point iteration given below  in Algorithm 3  to solve the nonlinear system.

		\medskip
		{\bf Algorithm 3.} Given ${\bf u}_h^{n+1,\ell-1} = {\bf u}_h^n$ for $n = 0, 1,\cdots, M-1$; $\ell = 1,\cdots,L$. Find $({\bf u}_h^{n+1,\ell}, p_h^{n+1,\ell}) \in L^2(\Omega;\mathbb{H}_h\times L_h)$ such that $\mathbb{P}$-a.s.
		\begin{align*}
			\bigl({\bf u}_h^{n+1,\ell},\pmb{\phi}_h\bigr) &+ \nu k\bigl(\nabla{\bf u}_h^{n+1,\ell},\pmb{\phi}_h\bigr) - k\bigl(p_h^{n+1,\ell},\div\pmb{\phi}_h\bigr) \\\nonumber
			&= \bigl({\bf u}_h^{n},\pmb{\phi}_h\bigr) - k\tilde{b}\bigl({\bf u}_h^{n+1,\ell-1},{\bf u}_h^{n+1,\ell-1},\pmb{\phi}_h\bigr) + \bigl(g_h^n\Delta {\bf W}_{n+1},\pmb{\phi}_h\bigr)\,\,\forall \pmb{\phi}_h\in\mathbb{H}_h,\\
			\bigl(\div{\bf u}_h^{n+1,\ell},\psi_h\bigr) &=0\qquad\forall \psi_h\in L_h.
		\end{align*}

		\medskip
		{\bf Test 1.} In this test, we want to verify the convergence rates in the error estimates \eqref{eq4.47}, \eqref{equu4.33}. To do that, the high moment errors $\mathcal{E}_{{\bf u},q}^M$ and $\mathcal{E}_{P,q}^M$ with $q = 2,4,8$ of the velocity and pressure approximations generated by Algorithm 2 are computed and the numerical results are shown in Table \ref{table5.1} and  \ref{table5.2}. The numerical results verify almost a half order convergence rate as predicted by our error estimate results.
		\begin{table}[tbhp]
			\begin{center}
				\begin{tabular}{ |c|c|c|c|c|c|c|}
					\hline
					\bf $k$ & $\pmb{\mathcal{E}}_{{\bf u},2}^M$  & \mbox{order} & $\pmb{\mathcal{E}}_{{\bf u},4}^M$ &  \mbox{order}&$\pmb{\mathcal{E}}_{{\bf u},8}^M$&order\\
					\hline
					$1/64$  & 0.01031 &  & 0.0136263&&0.0189242&\\
					\hline
					$1/128$  & 0.00873561 &0.2357  & 0.0118123&0.2068&0.0169163&0.1637\\
					\hline
					$1/256$  & 0.00618757 & 0.4975 & 0.007929&0.5751&0.010652&0.6673\\
					\hline
					$1/512$  & 0.0045004 & 0.4594 & 0.00576594&0.4586&0.00776608&0.4555\\
					\hline
				\end{tabular}
				\caption{Algorithm 2: Errors of the computed velocity $\{ {\bf u}^n_h\}_n$.}
				\label{table5.1}
			\end{center}
		\end{table}
		\begin{table}[tbhp]
			\begin{center}
				\begin{tabular}{ |c|c|c|c|c|c|c|}
					\hline
					\bf $k$ & ${\mathcal{E}}_{P,2}^M$  & \mbox{order} & ${\mathcal{E}}_{P,4}^M$ &  \mbox{order}&${\mathcal{E}}_{P,8}^M$&order\\
					\hline
					$1/64$  & 0.048219 &  & 0.0565676&&0.0648501&\\
					\hline
					$1/128$  & 0.0335814 & 0.5219 & 0.039123&0.5320&0.0446196&0.5397\\
					\hline
					$1/256$  & 0.02364 & 0.5064 & 0.0274858&0.5093&0.0312931&0.5118\\
					\hline
					$1/512$  &  0.0155695& 0.6025 &0.0180596 &0.6060&0.020537&0.5758\\
					\hline
				\end{tabular}
				\caption{Algorithm 2: Errors of the computed pressure $\{ {p}^n_h\}_n$.}
				\label{table5.2}
			\end{center}
		\end{table}

		{\bf Test 2.} The constant $C_q$ in Theorem \ref{thm4.8} and Theorem \ref{thm4.9} depends on $q$. However, no analytical formula for $C_q$ is known. The goal of this test is to find out numerically how the error constant $C_q$ depends on $q$. To the end, we fix $h = \frac{1}{40}$ and choose the time step $k = \frac{1}{16}$ to compute the errors $\mathcal{E}_{{\bf u},q}^M$ and $\mathcal{E}_{P,q}^M$ for different values of $q$. The numerical results are given  in Figure \ref{figure5.1} and \ref{figure5.2} (left). We observe that the numerical results  suggest the constant $C_q$ is increasing (and blows up) in $q$. For a close-up view, we zoom in at large
		$q$ to see the behaviors of both velocity and pressure errors in Figure \ref{figure5.1} and \ref{figure5.2} (right). We still see the increasing  of the errors in $q$ although the  growth becomes slower for large $q$. A consequence of this analysis also shows that we can not simply take limit as $p\to \infty$ in the high moment error estimates of Theorem \ref{thm4.8} and \ref{thm4.9} to derive pathwise error estimates \eqref{equu5.2} and \eqref{eq5.104}, and using Kolmogorov's Theorem is still the only viable approach for obtaining pathwise error estimates.
		\begin{figure}[thb]
			\begin{center}
				\centerline{
					\includegraphics[scale=0.4]{NS/VEQ_NS}
					\includegraphics[scale=0.4]{NS/VEQ_NS_Zoom}
				}
				\caption{Errors of the velocity approximation (left) in $\pmb{\mathcal{E}}^M_{{\bf u},q}$ norm and its behavior (right) for large $q$.}\label{figure5.1}
			\end{center}
		\end{figure}
		\begin{figure}[thb]
			\begin{center}
				\centerline{
					\includegraphics[scale=0.4]{NS/VEP_NS}
					\includegraphics[scale=0.4]{NS/VEP_NS_Zoom}
				}
				\caption{Errors of the pressure approximation (left) in $\mathcal{E}_{P,q}^M$ norm and its  behavior (right)  for large $q$.}\label{figure5.2}
			\end{center}
		\end{figure}

		\medskip
		{\bf Test 3.} In this test, we verify the pathwise error estimates \eqref{equu5.2} and \eqref{eq5.104}. We choose five sample paths and compute their $L^2$-norm errors for both velocity and pressure approximations. The numerical results are displayed in Figure \ref{fig5.3} which verify that the pathwise convergence of the velocity and pressure approximations is almost of order $O(k^{\frac12})$ as predicted by our pathwise error estimates.
		\begin{figure}[thb]
			\begin{center}
				\centerline{
					\includegraphics[scale=0.4]{NS/Pathwise_NS_V}
					\includegraphics[scale=0.4]{NS/Pathwise_NS_P}
				}
				\caption{Five sample pathwise errors of the velocity (left) and pressure approximations (right).}\label{fig5.3}
			\end{center}
		\end{figure}