APPROXIMATION OF STATIONARY STATISTICAL PROPERTIES OF THE THREE DIMENSIONAL AUTONOMOUS PLANETARY GEOSTROPHIC EQUATIONS OF LARGE-SCALE OCEAN CIRCULATION

. The main objective of this paper is to study the semi-implicit semi-discrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation. We prove the global attractor and stationary statistical properties of the semi-implicit semi-discrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation converge to those of the three di- mensional autonomous planetary geostrophic equations of large-scale ocean circulation as the time step goes to zero.

1. Introduction. In this paper, we mainly consider the temporal approximation of the following three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation (see [18,20]) in the domain where M ⊂ R 2 is a bounded domain with smooth boundary ∂M, h > 0 is the depth of the ocean. The unknown functions are the velocity field (v, w) = (v 1 , v 2 , w) = (v 1 (x, y, z, t), v 2 (x, y, z, t), w(x, y, z, t)), the pressure function p(x, y, z, t) and the temperature function T (x, y, z, t). The given function f = f 0 (β + y) is the Coriolis parameter, > 0 is a small constant, Q(x, y, z) is a heat source. Moreover, throughout this paper, we denote the two-dimensional horizontal gradient and Laplacian by ∇ and ∆, respectively, and the operators L 1 and L 2 are given by where A h , A ν are positive eddy viscosity constants and K h , K ν are positive conductivity constants. In order to state the boundary conditions of problem (1), for brevity, define Equations (1) are subject to the following boundary conditions with wind-driven on the top surface, no-penetration and non-flux on the side walls and bottom (see [15,20]) v · n| Γ l = 0, ∂v ∂ n × n Γ l = 0, ∂T ∂ n Γ l = 0 (2) and initial data T (x, y, z, 0) = T 0 (x, y, z), where τ (x, y) is the wind stress, n is the unit outward normal vector on Γ l , α > 0 is a positive constant and T * (x, y) is the typical temperature of the top surface, which satisfies the compatibility boundary condition The inviscid planetary geostrophic equations are derived from the Boussinesq equations for the planetary scale ocean using standard scale analysis (see [14,15,16,17,32]). During the past several decades, many authors have considered the well-posedness and the long-time behavior of solutions for the three dimensional planetary geostrophic equations of large-scale ocean circulation (see [1,18,19,33,34,35,36]). In particular, the authors in [1] proved the well-posedness of solutions and the existence of a finite dimensional global attractor in H. In [34,36], the authors proved the existence of pullback attractors for the three dimensional non-autonomous planetary geostrophic equations of large-scale ocean circulation. Moreover, the existence of random attractors for the three dimensional stochastic planetary geostrophic viscous equations of large-scale ocean circulation and its upper semi-continuity properties were established in [33,35].
Many dissipative dynamical systems arising in physical applications possess very complex behavior with abundant instability and sensitive dependence on initial data as well as parameters (see [24]). It is well-known that statistical properties of these kinds of systems are much more important, physically relevant and stable than single trajectories (see [7,10,11,12,13,28]).
For a given abstract autonomous continuous in time dynamical system determined by a semigroup {S(t)} t≥0 on a separable metric space H, we known that if the system reaches a statistical equilibrium in the sense that the statistics are time-independent stationary statistical properties, the probability measure µ on H describing the stationary statistical properties can be characterized via either the strong (pull-back) or weak (push-forward) formulation (see [7,11,12]).
Let {S(t)} t≥0 be a continuous semigroup on a metric space H which associates with a dynamical system on H. A Borel probability measure µ on H is called an invariant measure (stationary statistical solution) of the dynamical system, if where B(H) represents the σ-algebra of all Borel sets on H. Equivalently, the invariant measure µ can be characterized through the following push-forward weak invariance formulation for all bounded continuous test functionals Φ.
Invariant measure (stationary statistical solution) for a discrete dynamical system generated by a map S discrete on a metric space H is defined in a similar fashion with the continuous time t replaced by discrete time n = 0, 1, 2, · · · .
The stationary statistical properties of the discrete dynamical system converge to those of the continuous dynamical system if the invariant measures converge in the weak sense. We are usually interested in the statistical average H Φ(u)dµ(u) for various test functionals Φ, these averaged quantities are also called observables in physics literatures. Due to the presumed complexity of the dynamics, the physically interesting stationary statistical properties need to be calculated using numerical methods in generic case. Even under the ergodicity assumption, it is not at all clear that classical numerical schemes which provide accurate approximations on finite time intervals will remain meaningful for stationary statistical properties (long time properties), since small errors will be amplified and accumulated over long time, except in the case that the underlying dynamics is asymptotically stable, where statistical approach is not necessary since there is no chaos. Addressing issues like this is of great importance in many real life applications such as numerical study of climate change, since the climate is the long time statistical property of the underlying system. Therefore, it is central and of great challenge to search for numerical methods that are able to capture stationary statistical properties of infinite dimensional complex dynamical system. In a series of recent works, Professor X. Wang proposed a general framework of constructing temporal approximations of dissipative systems such as the 2D Rayleigh-Bénard convection system, such that the stationary statistical properties of the numerical scheme converge to those of the underlying Boussinesq system (see [2,3,26,29,30,31]). The main contribution of this article is the application of the general theory proposed in [30,31] to some kind of semi-implicit semi-discrete scheme in time of the three dimensional planetary geostrophic equations of large-scale ocean circulation. In the future, we may consider the case of fully discrete scheme in space and time.
The main purpose of this paper is to study some kind of semi-implicit semidiscrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation. Motivated by the ideas of preservation of dissipativity proposed in [5,6,8,9,21,22,23,25,26,27], in this paper, we first construct the temporal approximations of the three dimensional planetary geostrophic equations of large-scale ocean circulation that guarantee the convergence of the stationary statistical properties is the preservation of the dissipativity in some appropriate sense, and then we prove the existence of global attractors for the semiimplicit semi-discrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation. Finally, we prove that the global attractor and stationary statistical properties of the semi-implicit semidiscrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation converge to those of the underlying three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation as the time step goes to zero by verifying the conditions of the abstract results for temporal approximations of stationary statistical properties for dissipative dynamical systems established in [30].
Throughout this paper, let X be a Banach space with norm · X , let u p be the L p (Ω)-norm of u and let C be positive constants which may be different from line to line.

New formulations.
Integrating the third equation of (1) in the z direction and combining (2), we get Integrating the second equation of (1) with respect to z, we obtain where p s (x, y, t) is a free function with M p s (x, y, t) dxdy = 0.
2.2. Some function spaces. To study problem (8), we will introduce the notations for function spaces on Ω as follows for any v ∈ V 1 , T ∈ V 2 and let X be the dual space of a Banach space X with the dual action ·, · . Let H 1 , H 2 be the closure of V 1 , V 2 with respect to the following norms In what follows, we recall two Lemmas used to perform some a priori estimates in the sequel.
If ∆t n G n ≤ 1 − δ for any n ≥ n 0 , then 3. The existence of global attractors.
3.1. The well-posedness of the semi-discrete scheme in time. We consider the following semi-implicit semi-discrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation DenoteT m = T m −T * , thenT m+1 satisfies the following semi-implicit semi-discrete scheme in time and boundary conditions as well as initial data where Q * = Q + K h ∆T * .
The well-posedness of solutions for the semi-implicit semi-discrete scheme in time of the three dimensional autonomous planetary geostrophic equations of large-scale ocean circulation (9) can be obtained by Lax-Milgram Theorem. Proof. We will prove this theorem by induction. Since T 0 = T 0 ∈ V 2 , assume that T m ∈ V 2 is known, we want to prove that problem (9) admits a weak solution T m+1 ∈ V 2 . Since T m ∈ V 2 , we deduce from the Lax-Milgram Theorem and the regularity theory of elliptic partial differential equations of second order that there exists a unique weak solution v m ∈ V 1 of the first equation of problem (9) with v m ∈ V 1 ∩ H 2 (Ω). Moreover, it satisfies the following estimates: and v m 2 Denote by T =T m+1 , then problem (10) can be written as for any T, ψ ∈ V 2 , then we have the following results: Moreover, for any ψ ∈ V 2 , we have which implies that F ∈ V 2 . Therefore, we infer from the Lax-Milgram Theorem that there exists a unique weak solution T ∈ V 2 such that for any ψ ∈ V 2 , which entails that problem (9) admits a unique weak solution T m+1 = T + T * ∈ V 2 . Moreover, we can define a discrete dynamical system S k defined on the Hilbert space H 2 given by T m+1 = S k (T m ).

3.2.
The existence of global attractors. In this subsection, we will carry out some a priori estimates of solutions for problem (9), which imply the existence of global attractors in H 2 and V 2 for problem (9).

H 2 estimates of T m .
Theorem 3.2. Assume that τ ∈ H 1 0 (M ), Q ∈ L 2 (Ω) and T * ∈ H 2 (M ). Then there exists a constant ρ 1 independent of the time step k such that problem (9) possesses an absorbing ball in H 2 with radius ρ 1 which absorbs all bounded sets in H 2 .
Proof. Taking the inner product of the third equation of (9) with T m+1 in H 2 and using the fourth equation of (9), we obtain where we use the following equality Combining inequality (14) with the Sobolev trace embedding Theorem, Young's inequality, Lemma 2.1, we obtain which implies that From Lemma 2.2, we deduce T m 2 2 ≤ ( which entails that there exists a N ∈ N such that for any m ≥ N, where Moreover, we obtain for all L ≥ N + 1, where ρ 2 = (k(K + 1) + 1)ρ 1 .

V 2 estimates of T m .
Theorem 3.3. Assume that τ ∈ H 1 0 (M ), Q ∈ L 2 (Ω) and T * ∈ H 2 (M ). Then there exists a constant ρ 2 independent of the time step k such that problem (9) possesses an absorbing ball in V 2 with radius 2ρ 3 + 2h T * 2 H 2 (M ) , which attracts all bounded sets in H 2 .
Moreover, we obtain for all L ≥ N + K 0 + 2, where Thanks to the compactness of V 2 ⊂ H 2 , we immediately obtain the existence of a global attractor in H 2 of problem (9). From inequality (21) and Theorem 3.4, we infer that Theorem 3.5. Assume that τ ∈ H 1 0 (M ), Q ∈ L 2 (Ω) and T * ∈ H 2 (M ). Then the set Σ = 0<k≤k0 A k is uniformly bounded in V 2 , where A k is the global attractor in H 2 for problem (9) established in Theorem 3.4. satisfies the following equations:
Remark 1. In fact, we can obtain the existence of a global attractor in H 2 (Ω) ∩ V 2 of problem (9) by asymptotical a priori estimates under the assumptions that τ ∈ H 1 0 (M ), Q ∈ L 2 (Ω) and T * ∈ H 2 (M ).

Convergence of global attractors and stationary statistical properties.
In this section, we prove the convergence of global attractors and stationary statistical properties for problem (9) by verifying the conditions of Theorem 1 and Proposition 1 in [30].

4.1.
Finite time uniform convergence. In this subsection, we prove the finite time uniform convergence for problem (9).
Proof. In the following, we rewrite problem (9) into the following form where For any fixed S > 0, T 0 ∈ Σ and (L + 1)k ≤ S, we deduce from inequalities (11), (12), (15) and (18) that for any L ≤ L, These facts and inequality (21) imply that T k (t) and θ k (t) are uniformly in k bounded in L ∞ (0, S; V ) and the bounds are uniform in k and T 0 ∈ Σ. We can reformulate the differential form of the scheme (31) as where For t ∈ [mk, (m + 1)k), we obtain y, z)).
Therefore, for T 0 ∈ Σ, we have and p(t) = p s (x, y, t) − q k (t), then (u(t), ξ(t), p(t)) satisfies the following equations Taking the inner product of the third equation of (33) with ξ in L 2 (Ω) and using the fifth equation of (33), we obtain It is easy to verify that and u H 2 (Ω) ≤ C( ξ + g k 2 ).
We infer from inequalities (37)-(40) that Combining Theorem 1 and Proposition 1 in [30] with Theorem 3.5, Theorem 4.1 and Theorem 4.2, we immediately obtain the following result. Theorem 4.3. Assume that τ ∈ H 1 0 (M ), Q ∈ L 2 (Ω) and T * ∈ H 2 (M ). Then the stationary statistical properties as well as the global attractor of problem (9) convergence to those of problem (8) in the sense of Theorem 1 and Proposition 1 in [30].