Global Well-Posedness of the Ocean Primitive Equations with Nonlinear Thermodynamics

Abstract

We consider the hydrostatic Boussinesq equations of global ocean dynamics, also known as the “primitive equations”, coupled to advection–diffusion equations for temperature and salt. The system of equations is closed by an equation of state that expresses density as a function of temperature, salinity and pressure. The equation of state TEOS-10, the official description of seawater and ice properties in marine science of the Intergovernmental Oceanographic Commission, is the most accurate equations of state with respect to ocean observation and rests on the firm theoretical foundation of the Gibbs formalism of thermodynamics. We study several specifications of the TEOS-10 equation of state that comply with the assumption underlying the primitive equations. These equations of state take the form of high-order polynomials or rational functions of temperature, salinity and pressure. The ocean primitive equations with a nonlinear equation of state describe richer dynamical phenomena than the system with a linear equation of state. We prove well-posedness for the ocean primitive equations with nonlinear thermodynamics in the Sobolev space ${{\mathcal {H}}^{1}}$. The proof rests upon the fundamental work of Cao and Titi (Ann. Math. 166:245–267, 2007) and also on the results of Kukavica and Ziane (Nonlinearity 20:2739–2753, 2007). Alternative and older nonlinear equations of state are also considered. Our results narrow the gap between the mathematical analysis of the ocean primitive equations and the equations underlying numerical ocean models used in ocean and climate science.

1 Introduction

During the last two decades important progress has been made in the field of mathematical analysis of the primitive equations of large-scale atmosphere and ocean dynamics and in the field of ocean thermodynamics. These two areas appear so remote to each other that any exchange between them is minimal if even existing. The purpose of this work is to remedy this situation and to combine the accomplishments of both areas. This allows us to prove a result that is relevant for ocean modelling, because it formulates a fundamental well-posedness result for the ocean primitive equations now with realistic oceanic thermodynamics and brings thereby the analysis of the ocean primitive equations closer to actual codes that are used in ocean and climate science.

We consider the primitive equations in an ocean specific setting, i.e. together with advection–diffusion equations for temperature and salinity. These equations are closed by an equation-of-state that expresses density as a function of temperature, salinity and pressure. It is this equation-of-state that we are concerned with. The thermodynamic behaviour of seawater shows a complex interplay between salinity, temperature and pressure that determines the density. For the oceanic equation of state no analytical expression is available such as the ideal gas law for atmospheric thermodynamics. As a substitute one uses a linear equation of state in which the density is a linear combination of temperature and salinity.

So far the mathematical analysis of the primitive equations has been restricted exclusively to a linear equation of state. The analysis of the primitive equations dates back to the work of Lions-Temam-Wang [25], using a linear equations of state they have proven the existence of weak solutions in $L^\infty ([0,{\mathcal {T}}], {{\mathcal {L}}^{2}}(\varOmega ))\cap L^2([0,{\mathcal {T}}], {{\mathcal {H}}^{1}}(\varOmega ))$ for square integrable initial conditions. In their paper Lions–Temam–Wang remark that a more general polynomial equations of state “lead to further mathematical difficulties” (p. 1013 in [25]). In a seminal paper Cao and Titi [2] were the first to prove global existence and uniqueness of strong solutions in $L^\infty ([0,{\mathcal {T}}], {{\mathcal {H}}^{1}}(\varOmega ))\cap L^2([0,{\mathcal {T}}], {{\mathcal {H}}^{2}}(\varOmega ))$ for initial data in ${{\mathcal {H}}^{1}}$. They used an equation set consisting of velocity and temperature but the extension to a linear equation of state including salinity is straightforward. A different proof was given by Kobelkov [21]. For a different set of boundary conditions well-posedness in ${{\mathcal {H}}^{1}}$ was established in Kukavica-Ziane [22]. The authors solved first the mathematical challenging problems of the velocity equation without any thermodynamics before they indicated how to include a linear equation of state.

Modifications of the primitive equations that use only partial viscosity or diffusivity in several combinations were studied by Cao–Li–Titi in a series of papers [3,4,5,6,7]. A $L^p$-approach for the primitive equation was developed in Hieber and Kashiwabara [17], following the Fujita–Kato approach for the Navier–Stokes equations and without considering thermodynamics. This methodology allows to weaken the differentiability assumption on the initial conditions. This was extended in Hieber et al. [18] to the primitive equation with temperature and salinity using a linear equation of state. By using techniques from the theory of maximal $L^p$-regularity it was shown in Giga and Hieber [11] that the solution to the primitive equations becomes analytic for positive time for initial data in a Besov space. The existence of an attractor of the primitive equations was shown by Chuesov [8] and Ju [20]. For reviews on mathematical results of the primitive equations we refer to [16, 24, 29].

A linear equation of state is insufficient to represent the density dynamics in the world ocean. Density changes in the ocean are small but they have significant dynamical consequences, and although density changes act on small scales they shape the large scale circulation on long time scales. The density becomes dynamically active via the pressure gradient as vertical integral of density. It plays a central role in fundamental relations of geophysical fluid dynamics such as the geostrophic balance or the hydrostatic relation. The global circulation of the ocean is not only wind-driven but also density-driven. At high latitudes, for example, the cooling by the atmosphere increases the density of surface waters and this leads to sinking of water masses and convective overturning. These water masses flow in great depths southwards and get on their way slowly mixed. Important mixing processes such as cabbeling and thermobaricity are not captured by a linear equation of state. The downwelling of water must be compensated by an upwelling elsewhere. This resulting thermohaline circulation that forms the backbone of the global ocean circulation requires an adequate description of the oceanic density dynamics.

An oceanic equation of state has to include ocean phenomenology and to represent actual measurements of temperature, salinity, pressure and density to the highest possible degree and at the same time it has to strive for thermodynamic consistency. Thus, on one hand ocean thermodynamics has to consider observable quantities of real-world phenomena, on the other hand it has to reproduce fundamental and theoretical relations between abstract concepts such as entropy and enthalpy that are hardly measurable.

In the formulation of the ocean thermodynamics, as documented in the international standard TEOS-10 [19] (see e.g. [27] for a brief description), the internal consistency of the thermodynamic description is implemented by deriving all thermodynamic quantities from a general Gibbs potential function ${\mathbf {G}}$ (see [10] or [32] and references therein). From the Gibbs function all thermodynamics quantities are derived as mixed partial Frechet derivatives. Density, for example, is defined as partial derivative of the Gibbs function with respect to pressure. For the Gibbs function ${\mathbf {G}}={\mathbf {G}}(T, S, p)$ an Ansatz is made as high-order rational function of its dependent variables temperature, salinity and pressure. The phenomenology of the ocean is incorporated by taking a combination of partial derivatives of the Gibbs function that correspond to an observable quantity and determine the coefficients of the Gibbs function through a fit to the observational data of this quantity. We emphasize that this equation of state is not a general one, but one that was constructed for a very specific fluid, namely seawater. The procedure behind TEOS-10 that is sketched above leads to a formulation of ocean thermodynamics that is consistent and at the same time represents accurately real-world ocean observations.

The TEOS-10 thermodynamics is for several reasons not immediately applicable in numerical ocean models. First, the TEOS-10 equations of state is computationally quite expensive to calculate. The computation of density involves the calculation of a high-order polynomial at each grid cell, and also related quantities such as the thermal expansion and the saline contraction functions that are derivatives of the density with respect to temperature and salinity have the same computational cost. Second, the thermodynamic state variables in TEOS-10 are in-situ temperature and practical salinity, but in-situ temperature is not a useful variable in ocean circulation models. The third reason is that the primitive equations apply simplifications such as the Boussinesq approximation that are not made in TEOS-10. An approximation to the TEOS-10 equations of state that is suitable for the use in ocean models was developed in [30] and became part of the TEOS-10 standard.

Since the description of important oceanic phenomena necessitate an accurate—and that means a nonlinear—equation of state, it is important to include such an equation of state into the mathematical analysis of the primitive equations. This is the purpose of this paper. A nonlinear equation of state renders the pressure gradient in the velocity equations nonlinear with respect to temperature and salinity. As a result the velocity equation comprises now two nonlinearities, the velocity advection and the pressure gradient, where the latter one depends through the vertical integral of the density in a nonlinear fashion on temperature and salinity. In our formulation of the ocean primitive equations three types of pressure variables occur: the surface pressure $p_s$ that is depth independent, the hydrostatic pressure $p_{hyd}$ that is determined via the vertical integral of the density, and the static pressure $p_{stat}$ that is used in the equation of state to determine the density and that is a function of depth. The thermodynamic pressure is determined by energetic consistency considerations of the Boussinesq approximation.

In this paper we prove a well-posedness result for the ocean primitive equations with nonlinear thermodynamics given by the equation of state TEOS-10. As a corollary we obtain also well-posed for the predecessor of TEOS-10, the equation of state UNESCO-80. We consider the velocity boundary conditions that were used in the work of Cao and Titi [2] as well as the boundary conditions of Kukavica and Ziane [22]. Their a priori estimates provide the basis of our proof. We emphasize that the mathematically challenging part of the well-posedness proofs are the a priori estimates for the velocity advection (and surface pressure in case of mixed boundary conditions), while the estimates of the hydrostatic pressure gradient for a nonlinear equation of state are established straightforwardly. Although we follow here the energy-based approaches of [2, 22] we believe that the $L^p$- and the maximal regularity of approach of [11, 18] can without problems be generalized to include a nonlinear equation of state. This is not pursued here.

The paper is organized as follows: In Sect. 2 we describe the ocean primitive equations and specify the functional analytic setup. Section 3 introduces the basic notions of ocean thermodynamics and in the final Sect. 4 we state and prove the well-posedness theorem of the primitive equations with the nonlinear equation of state.

2 The Primitive Equations

The primitive equations are a central set of model equations for atmosphere and ocean dynamics. These equations approximate the large-scale dynamics of ocean and atmosphere and describe ocean or atmosphere as thin layer of a rotating and stratified fluid, together with the Boussinesq approximation and the hydrostatic approximation. The primitive equations are derived as an approximation of the rotating Navier–Stokes equations [23]. The state vector of the primitive equations consists of a horizontal velocity field $\mathbf{u}=(u,v)$ and temperature T and salinity S. We will say more about the temperature and salinity notions used in oceanography in Sect. 3.1.

Denote by $\varOmega :=M\times (-h,0)$ a cylindrical domain with a flat bottom, where $h>0$ denotes the depth of the ocean and $M\subseteq {\mathbb {R}}^2$ is bounded, with a smooth boundary $\partial M$. The primitive equations are given by

$$\begin{aligned}&\partial _t \mathbf{u}+(\mathbf{u}\cdot \nabla ) \mathbf{u}+w\partial _z \mathbf{u}+\nabla p +f\mathbf {k}\times \mathbf{u}-\frac{1}{Re_1}\triangle \,\mathbf{u}- \frac{1}{Re_2}\partial ^2_{zz} \mathbf{u}=0, \end{aligned}$$

(1a)

$$\begin{aligned}&\partial _z p +g\rho =0, \end{aligned}$$

(1b)

$$\begin{aligned}&\nabla \cdot \mathbf{u}+\partial _zw=0, \end{aligned}$$

(1c)

$$\begin{aligned}&\partial _t {T}+ (\mathbf{u}\cdot \nabla ) {T}+w\partial _z {T}-\frac{1}{Rt_1}\triangle {T}-\frac{1}{Rt_2}\partial ^2_{zz}{T}=0, \end{aligned}$$

(1d)

$$\begin{aligned}&\partial _t {S}+ (v\cdot \nabla ) {S}+w\partial _z {S}-\frac{1}{Rs_1}\triangle {S}-\frac{1}{Rs_2}\partial ^2_{zz}{S}=0,, \end{aligned}$$

(1e)

$$\begin{aligned}&\rho =\rho ({T},{S}, p), \end{aligned}$$

(1f)

where w denotes the vertical velocity, f the Coriolis parameter, and $\rho$ denotes an equation of state with pressure p. The numbers $Re_1,Re_2>0$ denote the horizontal and vertical Reynolds number that represent the viscosity coefficients. The numbers $Rt_1,Rt_2>0$ and $Rs_1,Rs_2>0$ are the horizontal and vertical mixing coefficients for temperature and salinity.

The operators $\nabla , \nabla \cdot$ and $\triangle$ denote the horizontal gradient, divergence and Laplacian, respectively. The partial derivative with respect to the vertical direction is denoted by $\partial _z$.

The vertical velocity w in (1) is determined by the constraint (1c) and can by virtue of the boundary conditions be expressed as

$$\begin{aligned} w(x,y,z,t)=-\int _{-z}^0\nabla \cdot \mathbf{u}(x,y,\xi ,t)\, d\xi =-\nabla \cdot \int _{-z}^0 \mathbf{u}(x,y,\xi ,t)\, d\xi . \end{aligned}$$

(2)

The theme of this paper is the pressure decomposition. The pressure term can with the hydrostatic balance be reformulated as

$$\begin{aligned} \begin{aligned}&\nabla p(x,y,z,t)=\nabla p_s(x,y,t) + \nabla p_{hyd}(x,y,z,t) ,\\ \text {with }&p_{hyd}:=\int _{-h}^z \rho [{T},{S},p](x,y,\xi ,t)\, d\xi \end{aligned}\end{aligned}$$

(3)

where $p_s$ denotes the surface pressure that acts as a Lagrange multiplier and has to be determined. The hydrostatic pressure $p_{hyd}$ is calculated from density $\rho$ via an evaluation of the equation of state. Through this equation of state are temperature and salinity equations coupled to the velocity equation. The nature of this coupling changes when one progresses from a linear to a nonlinear equation of state.

Boundary Conditions The statement of the primitive equations is completed by the description of the imposed boundary conditions. In this work here the external forcing is applied via the boundary conditions and not via the right hand side of the respective equation. We use one set of boundary conditions for temperature and salinity and consider two sets of boundary conditions for the velocity. Including a nonlinear equation of state has no consequences on the boundary conditions. We consider both types of velocity boundary conditions to demonstrate that the two corresponding well-posedness theorems and their different proofs remain valid for nonlinear ocean thermodynamics.

The boundary $\partial \varOmega$ consists of

$$\begin{aligned}&\text{ a } \text{ surface } \varGamma _u:=\{(x,y,z)\in {\bar{\varOmega }}:z=0 \}\\ {}&\text{ a } \text{ lateral } \text{ boundary } \varGamma _s:=\{(x,y,z)\in {\bar{\varOmega }}:(x,y)\in \partial M \}\\ {}&\text{ a } \text{ bottom } \text{ boundary } \varGamma _b:=\{(x,y,z)\in {\bar{\varOmega }}:z=-h \} \end{aligned}$$

such that $\partial \varOmega =\varGamma _u\cup \varGamma _s\cup \varGamma _b$.

For temperature and salinity we impose the following boundary conditions

$$\begin{aligned} \text{ T-S } \text{ bound. } \text{ cond. } {\left\{ \begin{array}{ll} &{}{}on\ \varGamma _u: \ \nabla _3{T}\cdot n_3=-k_\theta ({T}-{T}^*), \nabla _3 {S}\cdot n_3=0,\\ &{}{}on\ \varGamma _s: \ \nabla _3{T}\cdot \mathbf {n}= \nabla _3 {S}\cdot \mathbf {n}=0, \\ &{}{}on\ \varGamma _b: \ \nabla _3{T}\cdot n_3=\nabla _3 {S}\cdot n_3=0,\ \end{array}\right. } \end{aligned}$$

(4)

where $\mathbf {n}$ is the normal vector at $\varGamma _s$, $n_3$ the vertical Cartesian unit vector, and ${T}^*$ is a given temperature field at the sea surface.

For velocity we consider two types of boundary conditions.

At the top we impose the boundary condition that the ocean is driven by wind stresses, for the vertical velocity we impose a rigid lid, on the sides we impose for horizontal velocity a no penetration boundary condition and a stress-free boundary condition (i.e. the tangent component of the stress vanishes at the boundary), at the bottom we also have the stress-free boundary condition for the horizontal velocity and for the vertical velocity a no penetration boundary condition:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}{}on\ \varGamma _u:\ \partial _z\mathbf {u}=h\tau , w=0,\\ &{}{}on\ \varGamma _s:\ (x,y)\in \partial M \}: \ \mathbf {u}\cdot \mathbf {n}=0, \\ &{}{}\qquad \partial _{\mathbf {n}}\mathbf {u}\times \mathbf {n}=0, w=0,\ \\ &{}{}on\ \varGamma _b:\ \partial _z\mathbf {u}=0, w=0, \end{array}\right. } \end{aligned}$$

(5)

where $\tau$ is a given 2D wind stress field. The second set of boundary conditions imposes also a rigid lid, but on the lateral sides a homogeneous Dirirchlet boundary condition, on the top the same wind-driven boundary condition as in (5) but on the bottom a homogeneous Dirichlet boundary condition:

$$\begin{aligned} {\left\{ \begin{array}{ll} &{}{}on\ \varGamma _u:\ \partial _z\mathbf {u}=h\tau , w=0,\\ &{}{}on\ \varGamma _s:\ (x,y)\in \partial M \}: \ \mathbf {u}=0,w=0,\\ &{}{}on\ \varGamma _b:\ \mathbf {u}=0, w=0, \end{array}\right. } \end{aligned}$$

(6)

The first set of velocity boundary conditions was considered in [2], the second one in [22].

At the ocean surface $\varGamma _u$ the boundary condition for velocity $\mathbf{u}$ and temperature ${T}$ can be modified by adding a $\tau$- and ${T}^*$-dependent contribution such that we can assume without loss of generality homogeneous boundary conditions at the ocean surface (see [25], Section 2.4., [2] p. 248, [1], p. 202)

$$\begin{aligned} \partial _z\mathbf {u}=0 \ \text{ and } \ \partial _z{T}=0,\quad \text{ on } \varGamma _u. \end{aligned}$$

(7)

We define in the following the fundamental function spaces for velocity, temperature and salinity

$$\begin{aligned} \begin{aligned} {\tilde{{V_{v}}}}:&=\bigg \{\mathbf{u}\in (\mathbf{C}^\infty ({\bar{\varOmega }}))^2 :\ \mathbf{u}\text { satisfies either boundary condition } (5)\text { or }(6), \\&\qquad \text { and }\nabla \cdot \mathbf{u}+\partial _z w=0 \bigg \},\\ {\tilde{{V_{\theta }}}}:&=\bigg \{{T}\in C^\infty ({\bar{\varOmega }}) :{T}\text { satisfies boundary condition} (4) \text { for temperature } \bigg \},\\ {\tilde{{V_{S}}}}:&=\bigg \{{S}\in C^\infty ({\bar{\varOmega }}) :\theta \text { satisfies boundary condition } (4)\text { for salinity, }\\&\qquad \text {and }\int _\varOmega S\, dxdydz=0 \bigg \}. \end{aligned}\end{aligned}$$

(8)

Denote by ${V_{v}}$, ${V_{\theta }}, {V_{S}}$ the closure of ${\tilde{{V_{v}}}}$ in the Sobolev space $\mathbf{H}^{1}(\varOmega )$ and of ${\tilde{{V_{\theta }}}}$ in ${H^{1}}(\varOmega )$. The norms are for $\mathbf{u}\in {V_{v}}$, $T\in {V_{\theta }}$ and $S\in {V_{S}}$ defined as

$$\begin{aligned} \begin{aligned}&||\mathbf{u}||_{\mathbf{H}^{1}}^2:=\int _\varOmega |\mathbf{u}(x,y,z)|^2dxdydz+\int _\varOmega |\nabla _3\mathbf{u}(x,y,z)|^2dxdydz,\\ \quad&||{T}||_{{H^{1}}}^2:=\int _\varOmega |{T}(x,y,z)|^2dxdydz+\int _\varOmega |\nabla _3{T}(x,y,z)|^2dxdydz\\&\qquad \qquad +k_{T}\int _{\varGamma _u}|\theta (x,y,0)|^2dxdy,\\ \text {and}\quad&||{S}||_{{H^{1}}}^2:=\int _\varOmega | {S}(x,y,z)|^2dxdydz+\int _\varOmega |\nabla _3 {S}(x,y,z)|^2dxdydz. \end{aligned}\end{aligned}$$

(9)

The product space denoted by

$$\begin{aligned} \begin{aligned} {\mathcal {V}}:={V_{v}}\times {V_{\theta }}\times {V_{S}}. \end{aligned} \end{aligned}$$

(10)

We denote the $L^2$-norm and the $H^2$-norm in the velocity and temperature, salinity spaces by $||\cdot ||_{\mathbf{L}^{2}}$, $||\cdot ||_{{L^{2}}}$ and by $||\cdot ||_{\mathbf{H}^{2}}$, $||\cdot ||_{{H^{2}}}$, respectively.

Definition 1

(Strong Solution). i) Let initial conditions $(\mathbf{u}_0,{T}_0, {S}_0)\in {\mathcal {V}}$ be given, and let $[0,\mathcal {T}]$, $\mathcal {T}>0$, be a time interval. Let the temperature and salinity boundary conditions be given by (4), and the velocity boundary condition either given by (5) or by (6). The state vector $(\mathbf{u},{T},{S})$ is called a strong solution of (1) on [0, T] if

$$\begin{aligned}\begin{aligned}&\mathbf{u}\in C([0,{T}],{V_{v}})\cap L^2([0,{T}], \mathbf{H}^{2}(\varOmega )),\\&{T},{S}\in C([0,{T}],{V_{\theta }})\cap L^2([0,{T}],{H^{2}}(\varOmega )),\\&\partial _t \mathbf{u}\in L^2([0,{T}], \mathbf{L}^{2}(\varOmega )),\\&\partial _t ({T},{S})\in L^2([0,{T}], L^2(\varOmega )), \end{aligned}\end{aligned}$$

and if it satisfies in the sense of distributions for all $\varPhi \in (C^\infty ({\bar{\varOmega }}))^2$ and $\phi \in C^\infty ({\bar{\varOmega }})$

$$\begin{aligned}&\int _\varOmega \partial _t \mathbf{u}\cdot \varPhi \, dx +\int _\varOmega [(\mathbf{u}\cdot \nabla ) \mathbf{u}]\cdot \varPhi \, dx+\int _\varOmega f\mathbf {k}\times \mathbf{u}\varPhi \, dxdydz +\int _{\varGamma _u}(\nabla p_s)\cdot \varPhi \,dxdydz\nonumber \\&-\int _\varOmega \big [ \big (\nabla \cdot \int _{-h}^z \mathbf{u}(x,y,\xi ,t)\, d\xi \big )\partial _z \mathbf{u}\big ]\cdot \varPhi \, dxdydz \nonumber \\&-\int _\varOmega \big [\nabla \int _{-h}^z g\rho (x,y,\xi ,t)\,d\xi \big ]\cdot \varPhi \, dxdydz +\int _\varOmega \frac{1}{Re_1}\nabla \mathbf{u}\cdot \nabla \varPhi + \frac{1}{Re_2}\partial _z \mathbf{u}\cdot \partial _z \varPhi \, dxdydz =0, \end{aligned}$$

(11a)

$$\begin{aligned} \int _\varOmega&\partial _t {T}\phi \,dxdydz + \int _\varOmega (\mathbf{u}\cdot \nabla ) {T}\phi \,dxdydz -\int _\varOmega \big [\nabla \cdot \int _{-h}^0 \mathbf{u}(x,y,\xi ,t)\, d\xi \big ]\partial _z {T}\, \phi \,dxdydz \nonumber \\&+\frac{1}{Rt_1}\int _\varOmega \nabla _3 {T}\cdot \nabla _3\phi \,dxdydz +\frac{1}{Rt_2}\int _\varOmega \partial _z{T}\partial _z\phi \,dxdydz\nonumber \\&+k_{T}\int _{\varGamma _u}{T}\phi \, dxdy=0, \end{aligned}$$

(11b)

$$\begin{aligned} \int _\varOmega&\partial _t {S}\phi \,dxdydz + \int _\varOmega (\mathbf{u}\cdot \nabla ) {S}\phi \,dxdydz -\int _\varOmega \big [\nabla \cdot \int _{-h}^0 \mathbf{u}(x,y,\xi ,t)\, d\xi \big ]\partial _z {S}\, \phi \,dxdydz \nonumber \\&+\frac{1}{Rs_1}\int _\varOmega \nabla _3 {S}\cdot \nabla _3\phi \,dxdydz +\frac{1}{Rs_2}\int _\varOmega \partial _z {S}\partial _z\phi \,dxdydz=0, \end{aligned}$$

(11c)

$$\begin{aligned} \text {with }&\ \rho :=\rho ({T},{S},p), \end{aligned}$$

(11d)

where w denotes the vertical velocity, f the Coriolis parameter. The numbers $Re_1,Re_2>0$ denote the horizontal and vertical Reynolds number that represent the viscosity coefficients. The numbers $Rt_1,Rt_2>0$ and $Rs_1,Rs_2>0$ are the horizontal and vertical mixing coefficients for temperature and salinity.

Theorem 1

(Global Well-Posedness [2, 22]). Let a time interval $[0, \mathcal {T} ]$, with $\mathcal {T} > 0$, be given. The boundary conditions for temperature and salinity are given by (4) and the boundary conditions for velocity are given either by (5) or by (6). Let the equation of state in (1f) given by

$$\begin{aligned} \rho :=\rho _0-\alpha ({T}-{\bar{{T}}})+\beta ({S}-{\bar{{S}}}), \end{aligned}$$

where $\rho _0, {\bar{{T}}},{\bar{{S}}}>0$ are a mean density, a mean temperature and a mean salinity, respectively. The numbers $\alpha ,\beta >0$ denote the thermal expansion and saline contraction coefficient. If the initial conditions for velocity, temperature and salinity satisfy $(\mathbf{u}_0, {T}_0,{S}_0) \in {\mathcal {V}}$, then there exists on the interval $[0, \mathcal {T} ]$ a unique strong solution $(\mathbf{u}, {T},{S})$ of the system (1) in the sense of Definition 1. This solution depends continuously on the initial data with respect to the $C([0,\mathcal {T}], L^2(\varOmega ))$-topology.

The regularity of solutions of the primitive equations, described in Theorem 1, still leaves room for a broad spectrum of spatio-temporal scales that comprises local and sudden events with large gradients such as fronts. For an analysis of the admissible scales in terms of Nusselt, Raleigh and Reynolds numbers we refer to [13, 14].

3 Ocean Thermodynamics

In what follows we need different notions of temperature and salinity. We collect our notation here. The different concepts are introduced below. For details and extensive description we refer to [10, 19].

Notation

• In-situ temperature: ${T_{is}}$

• Potential temperature: ${\theta _{pot}}$

• Conservative temperature: ${\varTheta }$

• Generic temperature variable (one of the three above): ${T}$

• Practical salinity: ${S_{prac}}$

• Absolute salinity: $\mathbf{S}$

• Generic salinity variable (one of the two above): ${S}$

3.1 Thermodynamic State Variables

Temperature: Oceanography uses several concepts of temperature such as in-situ temperature “${T_{is}}$”, potential temperature “${\theta _{pot}}$” or conservative temperature “${\varTheta }$”. In-situ temperature is the temperature that is actually measured at a certain position. It has, by definition, the advantage of being directly measurable but the disadvantage that is changes with pressure and does therefore not give an adequate picture of the heat content. The concept of potential temperature with respect to a reference pressure, resolves this problem. It is defined as the temperature a fluid parcel would have if it is transported without entropy change from its in-situ pressure to a reference pressure (typically the pressure at the ocean surface). If $\eta$ denotes the entropy, and ${T_{is}},p, {S}$ the in-situ temperature, salinity and pressure, respectively, then the potential temperature ${\theta _{pot}}$ is defined by the equation $\eta ({\theta _{pot}},{S},p_r)=\eta ({T_{is}},{S},p)$, where $p_r$ denotes the reference pressure. An even more accurate image of the heat content is given by the conservative temperature ${\varTheta }$, which is defined in terms of potential enthalpy H. Potential enthalpy is the heat content a fluid parcel would have if it is brought from its actual pressure to a reference pressure $p_r$ without energy exchange. The conservative temperature ${\varTheta }$ is now defined via the enthalpy as ${\varTheta }:=\frac{H(\theta ^*,S,p_r)}{c_p}$, with the heat capacity $c_p:=H(\theta ^*=25^\circ C, 35psu, p_r)/25\,^\circ \hbox {C}$. Conservative temperature is more accurate than potential temperature but the differences are in general relatively small. For ocean general circulation models the notion of potential temperature as well as of conservative temperature is appropriate. Potential temperature is the traditional notion used in ocean models, while conservative temperature is a more recent concept. To which extend the change from potential to conservative temperature improves ocean models in the sense that it reduces biases and decrease the difference to observations is subject of ongoing research

The evolution equation for potential temperature or conservative temperature have the same mathematical structure of an advection diffusion equation as in (1d). The equation for in-situ temperature is different (see (2.101) in [28]). This simplifies the analysis, but it is important to distinguish which temperature variable is used in TEOS-10.

Salinity: TEOS-10 use the concept of absolute salinity $\mathbf{S}$. Absolute Salinity is defined as the mass fraction of dissolved material in seawater. It has the advantage over other notions such as practical salinity “${S_{prac}}$” that it is influenced by the mass of dissolved constituents, whereas Practical Salinity depends on the conductivity of seawater. For a comprehensive discussion of salinity notions we refer to TEOS-10 [19].

The advection–diffusion equation for the different concepts of temperature and salinity outlined above all have the familiar structure shown in (1d) and (1e). The different definitions show up in different units and in different mixing coefficients. In this work we use nondimensional quantities $Rt_1,Rt_2$ and $Rs_1,Rs_2$ as mixing coefficients, such that there is no difference in the equations noticeable.

Pressure: Sea pressure p is defined to be the absolute pressure P, i.e. pressure relative to a vacuum, minus the absolute pressure of one standard atmosphere, $P_0 = 101325$ Pa, that is $p := P-P_0$. The thermodynamic representation of a system can in principle alternatively formulated in terms of different variable such as specific volume, entropy and salinity. Here we follow the TEOS-10 convention and use temperature, salinity and pressure.

3.2 The TEOS-10 Thermodynamics

The Gibbs potential function $\mathbf {G}$ in TEOS-10 is defined in terms of in-situ temperature, absolute salinity and pressure by

$$\begin{aligned} \begin{aligned} {\mathbf {G}}: ({T_{is}},\mathbf{S},p)&\mapsto {\mathbf {G}}({T_{is}},\mathbf{S},p), \end{aligned} \end{aligned}$$

(12)

where the Gibbs function ${\mathbf {G}}$ consists of two terms, one that describes the liquid pure water ${\mathbf {G}}^W$ and a second one ${\mathbf {G}}^S$ that takes into account salinity

$$\begin{aligned} \begin{aligned} {\mathbf {G}}({T_{is}},\mathbf{S},p)={\mathbf {G}}^W({T_{is}},p)+{\mathbf {G}}^S({T_{is}},\mathbf{S},p), \end{aligned}\end{aligned}$$

(13)

with ${\mathbf {G}}^s({T_{is}},0,p)=0$. TEOS-10 provides also a Gibbs potential function for ice but this is not pursued here. The two Gibbs functions ${\mathbf {G}}^W,{\mathbf {G}}^S$ are polynomials in the respective variables

$$\begin{aligned} \begin{aligned} {\mathbf {G}}^W({T_{is}},p)&=\sum _{j=0}^7\sum _{k=0}^{6}g_{j,k}{T_{is}}^j p^k,\\ {\mathbf {G}}^S({T_{is}},\mathbf{S},p),&=\sum _{i=0}^7\sum _{j=0}^6\sum _{k=0}^{5} \bigg \{g_{1jk}\mathbf{S}^2\ln \mathbf{S}+\sum _{i>1}g_{i,j,k}\mathbf{S}^i\bigg \} {T_{is}}^j p^k \end{aligned}\end{aligned}$$

(14)

From the Gibbs function fundamental thermodynamic quantities are derived in a thermodynamic consistent way as partial Frechet-derivatives. The density $\rho$ is defined via the specific volume $\upsilon$ and the specific volume is the partial Frechet-derivative of ${\mathbf {G}}$ with respect to the pressure variable

$$\begin{aligned} \begin{aligned} \rho ({T_{is}},\mathbf{S},p):=\frac{1}{\upsilon ({T_{is}},\mathbf{S},p)}\qquad \text {with}\qquad \upsilon ({T_{is}},\mathbf{S},p):=\frac{\partial \mathbf {G}}{\partial p}({T_{is}},\mathbf{S},p). \end{aligned}\end{aligned}$$

(15)

The ocean phenomenology enters the picture through the determinations of the coefficients $g_{ijk}$ of the Gibbs function (13). These are are determined by a least-square fit of partial derivatives of the Gibbs function to observations or laboratory measurements (for details we refer to [19]). This is the achievement of the thermodynamic formalism used in TEOS-10, it keeps the internal consistency of thermodynamics using a generic form of rational functions for the Gibbs potential and at at the same time integrates oceanic phenomenology by specifying the coefficients of the rational functions as an optimal fit to observations.

The TEOS-10 standard provides a description of the thermodynamics of seawater that is suitable for the whole spectrum of marine science. Numerical modelling of the large scale ocean dynamics, based on the primitive equations, comes with some simplifications such as the hydrostatic and the Boussinesq approximation. The Boussinesq approximation is thermodynamically relevant and will be discussed in Sect. 3.4. For actual ocean codes computational efficiency is of crucial importance and this places limitations on the computational resources spend for the evaluation of an equation of state. The computation of the density itself is costly as it involves the calculation of a high-order polynomial at each grid cell. Finally the Gibbs formalism of TEOS-10 is based on in-situ temperature and this requires a conversion in potential or conservative temperature.

Therefore TEOS-10 contains a formulation of the equation of state that is more suitable for ocean models, because it is formulated in terms of conservative temperature ${\varTheta }$ instead of the in-situ temperature ${T_{is}}$ and that is also more efficiently to compute than ${\mathbf {G}}$ in (13). Following TEOS-10 (see Appendix K in [19]) the specific volume $\upsilon$ is modelled in this case with a 75-term polynomial and the equation of state as rational function of the form

$$\begin{aligned} \begin{aligned} \rho _{rat}({\varTheta },\mathbf{S},p):=\frac{1}{\upsilon _{rat}({\varTheta },\mathbf{S},p)}\quad \text {and}\quad \upsilon _{rat}({\varTheta },\mathbf{S},p):=\sum _{i,j,k=0}^{N_i,N_j,N_k}c_{i,j,k}\mathbf{S}^i {\varTheta }^j p^k, \end{aligned}\end{aligned}$$

(16)

where $c_{i,j,k}\in \mathbb {R}$ are given coefficients (see table K1 in [19]) and $N_i=6, N_j=6, N_k=6$. The formulation of the equations of state (16) is based on [30]. In this work an additional performance gain was obtained by a modelling the equation of state directly with a 52-term polynomial

$$\begin{aligned} \begin{aligned} \rho _{pol}({\varTheta },\mathbf{S},p):=\sum _{i,j,k=0}^{M_i,M_j,M_k}d_{i,j,k}\mathbf{S}^i {\varTheta }^j p^k, \end{aligned}\end{aligned}$$

(17)

$d_{i,j,k}\in \mathbb {R}$ are given coefficients (see Appendix A in [30]) and $M_i=6, M_j=6, M_k=3$.

It is worth emphasizing that the oceanic equation of state does not describe a general fluid but a very special one, namely seawater. The equation of state in TEOS-10 are defined exclusively and only for seawater, aiming at capturing the unique relation between temperature,salinity, pressure and density that is measured in the world ocean. The thermodynamic formalism outlined above is designed for temperature, salinity and pressure values in a range that is observed in the real ocean. In the ocean the temperature varies between $\sim -2^\circ$ and $40^\circ$, salinity between 1 and 45 psu, the pressure up to 100 Mdba, the resulting density varies little around $1000\,\hbox {kg/m}^3$. The physically sensible range for temperature, salinity and pressure are approximately $I_\theta =[-2^\circ , 40^\circ ]$ and $I_S=[1\,\hbox {psu}, 45\,\hbox {psu}]$ and $I_p=[1,100\,\hbox {Mdb}]$. Relevant for us are not the particular numbers but that such a bounded interval for these fields exists and that the equation-of-state is constructed such that it provides corresponding density values that cover the range determined by ocean observations.

The ocean equation of state implicitly assumes that temperature, salinity and pressure values are within the observed range of the ocean. Outside this range the ocean density as well as the Gibbs function are not defined, because there are no observation available to close the defination of the Gibbs function. This motivates the following natural hypothesis on ocean thermodynamics:

Hypothesis on Ocean Thermodynamics:

• 1. There exist bounded intervals for temperature in the ocean $I_T:=[-K_1, K_2]$, for salinity $I_{S}:=[L_1,L_2]$, and for pressure $I_p:=[M_1,M_2]$ with $K_1,K_2$, $L_1,L_2,M_1,M_2>0$, such that for almost every $(x,y,z,t)\in \varOmega \times [0,\mathcal {T}]$ it holds

  $$\begin{aligned} \begin{aligned} T(x,y,z,t)\in I_T,\quad {S}(x,y,z,t)\in I_{S},\quad p(x,y,z,t)\in I_p. \end{aligned}\end{aligned}$$

  (18)

  We refer to the intervals $I_T, I_{S}, I_p$ as physically admissible ranges for the respective variable.

• 2. If $T\in I_T, {S}\in I_{S}$, $p\in I_p$, then the coefficients $c_{i,j,k}$ in (16) and $d_{i,j,k}$ in (17) are determined such that $\rho _{rat}, \rho _{pol}$ defined in (16), (17) satisfy

  $$\begin{aligned} \begin{aligned} \rho _{rat}, \rho _{pol}\in I_\rho :=[N_1,N_2], \quad N_1,N_2>0, \end{aligned}\end{aligned}$$

  (19)

  where $I_\rho$ is the physically admissible range for density.

The hypothesis above are valid for the different temperature and salinity variables defined in Sect. 3.1. In this work we focus on potential and conservative temperature, and practical and absolute salinity, because they are the most relevant variables for ocean general circulation models.

3.3 Different Nonlinear Equations of State

The focus of this work is the equation of state documented in the international standard TEOS-10. Since this standard is currently not fully implemented by all ocean general circulation models, we describe in this section also the UNESCO-80 equation of state [31], the predecessor of TEOS-10, and an intermediate equation of state due to McDougall et al. [26].

3.3.1 The UNESCO-80 Equation of State

The UNESCO-80 equations of state [31] is the oldest equations of state and the former international standard that is succeeded now by TEOS-10 (see [19] for shortcomings of the UNESCO-80 equation of state). The density is defined by

$$\begin{aligned} \begin{aligned} \rho _{\scriptscriptstyle UNO}({\theta _{pot}},{S_{prac}},p):=\frac{{\tilde{\rho }}({T_{is}},{S_{prac}})}{1-\frac{p}{K({T_{is}}, {S_{prac}},p)}}, \end{aligned}\end{aligned}$$

(20)

where ${\tilde{\rho }}$ and K are polynomials in their respective variables. For details we refer to [31] or Appendix 3 in [12]. For a description of the deficiencies of the UNESCO-80 equations of state we refer to [19].

3.3.2 The Equation of State of McDougall–Jackett–Wright–Feistel

in the literature one finds several suggestions to improve the accuracy of the UNESCO-80 equation of state that are used by codes of ocean general circulation models. One of these examples is due to McDougall et al. [26]. This equation of state is of the form

$$\begin{aligned} \begin{aligned} \rho _{\scriptscriptstyle MJWF}({\theta _{pot}},{S_{prac}},p):=\frac{P_1({\theta _{pot}},{S_{prac}},p)}{P_2({\theta _{pot}},{S_{prac}},p)}, \end{aligned}\end{aligned}$$

(21)

where $P_1,P_2$ are two polynomials with 12 and 13 terms respectively (see Appendix B in [26] for the actual coefficients). The equations of state (21) is formulated in terms of potential temperature rather in-situ temperature, a choice that facilitates its use in ocean models, where potential temperature is the most common temperature variable.

3.4 Thermodynamics for the Primitive Equations

The TEOS-10 equation of state describes the thermodynamics of seawater without any approximation while the ocean primitive equations rest on several fundamental simplifications. For thermodynamics the Boussinesq approximation is relevant. The Boussinesq approximation replaces the density in the dynamical equations in all terms except in the buoyancy term of the velocity equation and in the equation of state (see e.g. Section 2.4. in [32]). It is the second aspect that we are concerned with. The general equation of state (15) that is valid for non-Boussinesq oceanic models, depends on in-situ temperature, salinity and pressure. In the case of the primitive equations, i.e. of a hydrostatic Boussinesq fluid, there arises a subtlety which pressure p to be used in the equation-of state $\rho ({\varTheta },S,p)$. Using the hydrostatic pressure $p=p_h$, with $p_h$ from (3), then density would depend on hydrostatic pressure $\rho =\rho ({\varTheta },S,p_h)$ and the hydrostatic pressure $p_h=p_h(\rho )$ would depend on the density $\rho$, as a result we would have an implicit definition of the density through an nonlinear function of the form $\rho =\rho ({\varTheta },\mathbf{S},\rho )$. The solution of this dilemma comes from the Boussinesq approximation.

For the Boussinesq approximation the pressure in the equation of state has to be a function of the static pressure $g\rho _0 z$ with respect to a constant reference density $\rho _0>0$. As a consequence the pressure in the equation of state is a function of depth only. We refer to it as the static pressure and define

$$\begin{aligned} \begin{aligned} p_{stat}(z):=-g\rho _0z,\quad \text {for }z\le 0. \end{aligned}\end{aligned}$$

(22)

The reason for this choice is that for a general pressure, i.e. a pressure not only dependent on depth, no energy conservation law exists that takes into account that for Boussinesq fluids the concept of internal energy does not exist and that energy exchange is only possible between kinetic and potential energy. A different choice of the pressure in the equation of state would result in creating an internal energy reservoir. For more details we refer to Section 2.4.3. in [32], and to [33].

With the static pressure (22) the two equation of state $\rho _{rat}, \rho _{pol}$ in (16), (17) become now

$$\begin{aligned} \begin{aligned} \rho _{rat}({\varTheta },\mathbf{S},p)&\rightarrow \rho _{rat}({\varTheta }, \mathbf{S}, p_{stat}),\\ \rho _{pol}({\varTheta },\mathbf{S},p)&\rightarrow \rho _{pol}({\varTheta }, \mathbf{S}, p_{stat}), \end{aligned}\end{aligned}$$

(23)

We summarize this in the following definition.

Definition 2

Let the hypothesis on ocean thermodynamics (18) and (19) be satisfied.

(i):

Define $L^\infty _{I_\theta }(\varOmega ):=\{\theta \in L^\infty (\varOmega ): \theta \in I_\theta \}$ and $L^\infty _{I_S}(\varOmega ):=\{S\in L^\infty (\varOmega ): S\in I_S\}$, where $I_\theta , I_S$ denote the physically admissible range for potential temperature and salinity, introduced in (18). For $C\in \{\theta , S\}$ then $L^\infty _{I_C}(\varOmega )$ denotes either $L^\infty _{I_\theta }(\varOmega )$ or $L^\infty _{I_S}(\varOmega )$. The product space is defined by $L^\infty _{{I_\theta }\times I_{S}}(\varOmega ):=L^\infty _{I_\theta }(\varOmega )\times L^\infty _{I_S}(\varOmega )$.

(ii):

Let $\theta \in L^\infty _{I_\theta }(\varOmega ), S\in L^\infty _{I_S}(\varOmega )$. The density $\rho _{rat}$ is given by

$$\begin{aligned} \begin{aligned}&\rho _{rat}({\varTheta },\mathbf{S},p_h)=\upsilon _{rat}({\varTheta },\mathbf{S},p_{stat})^{-1},\\ \text {with }&\upsilon _{rat}({\varTheta },\mathbf{S}, p_{stat})=\sum _{i,j,k=0}^{N_i,N_j,N_k}c_{i,j,k}\mathbf{S}^i {\varTheta }^j p_{stat}^k, \end{aligned}\end{aligned}$$

(24)

where the static pressure $p_{stat}$ is given by (22).

(iii):

Let $\theta \in L^\infty _{I_\theta }(\varOmega ), S\in L^\infty _{I_S}(\varOmega )$. The density $\rho _{pol}$ is given by

$$\begin{aligned} \begin{aligned}&\rho _{pol}({\varTheta },\mathbf{S},p_h)=\sum _{i,j,k=0}^{M_i,M_j,M_k}d_{i,j,k}\mathbf{S}^i {\varTheta }^j p_{stat}^k, \end{aligned}\end{aligned}$$

(25)

where the thermodynamic pressure $p_{stat}$ is given by (22).

4 Well-Posedness of Primitive Equations with TEOS-10

Before we state the well-posedness theorem we provide two Lemma. The first lemma shows the fundamental property of temperature and salinity equations, namely that both quantities remain bounded within the range of there (bounded) initial values. This shows that if temperature and salinity values are initially within the physically sensible range, then the remain there at later times. This results implies the same property for density. Therefore it is sufficient to impose boundedness within a physically reasonable range to the initial temperature and salinity fields. This leads then to a density that is also in the physically admitted range, and is in particular bounded above and below and therefore integrable.

Lemma 1

(Boundedness of Temperature, Salinity and Density).

• (i) Denote by $C\in \{{\varTheta },\mathbf{S}\}$ the solution of the tracer equation for temperature (1d) or for salinity (1e). Suppose the initial conditions for temperature or salinity satisfy $C_0\in {V_{C}}\cap L^\infty _{I_C}(\varOmega )$. Then for almost every $(x,y,z,t)\in \varOmega \times [0,{T}]$

  $$\begin{aligned} C(x,y,z,t)\in I_C, \end{aligned}$$

  (26)

  where $I_C$ denotes the physically admissible range of values for temperature or salinity, i.e. $I_C=I_\theta$ if $C=\theta$ or $I_C=I_S$ if $C=S$, where $I_\theta , I_S$ were introduced in the fundamental assumptions on thermodynamics (see (18)).

• (ii) Let ${\varTheta }_0\in {V_{\theta }}\cap L^\infty _{I_\theta }(\varOmega )$, $\mathbf{S}_0\in {V_{S}}\cap L^\infty _{I_S}(\varOmega )$ be initial conditions for the equations for temperature (1d) and for salinity (1e). Then the densities $\rho _{rat}, \rho _{pol}$ are in the physically admissible range $\rho _{rat}, \rho _{pol}\in I_\rho$ for almost every $(x,y,z,t)\in \varOmega \times [0,{T}]$.

Proof

$Ad\, i)$ Define $\hat{C\!}\,(t):=C(t)-C_0$ and $\hat{C\!}_+:=\max \{0,\hat{C\!}\,\}$, then $\hat{C\!}\,$ satisfies

$$\begin{aligned} \begin{aligned} \partial _t \hat{C\!}\,+ (v\cdot \nabla ) \hat{C\!} +w\partial _z \hat{C\!} =\frac{1}{Rc_1}\triangle \hat{C\!}\,+\frac{1}{Rc_2}\partial _{zz}^2\hat{C\!}, \end{aligned}\end{aligned}$$

(27)

where $Rc_1\in \{Rt_1,Rs_1\}, Rc_2\in \{Rt_2,Rs_2\}$. Taking the scalar product in $L^2$ with $\hat{C\!}_+$ yields after integration by parts

$$\begin{aligned} \begin{aligned}&d_t ||\hat{C\!}_+||_{L^2} +\frac{1}{Rc_1}\int _\varOmega |\nabla \hat{C\!}_+|^2, dxdydz +\frac{1}{Rc_2}\int _\varOmega |\partial _z\hat{C\!}_+|^2, dxdydz\\&\quad +\,k_C||\hat{C\!}_+(z=0)||_{L^2} =0, \end{aligned}\end{aligned}$$

(28)

where $k_C=k_\theta$ in the case $C=\theta$ and where $k_C=0$ if $C=S$ (cf. boundary conditions (4). Neglecting the positive diffusion terms on left-hand side yields

$$\begin{aligned} \begin{aligned} d_t ||\hat{C\!}_+||_{L^2} +k_C||\hat{C\!}_+(z=0)||_{L^2} \le 0. \end{aligned} \end{aligned}$$

(29)

This implies $\hat{C\!}_+(t)\equiv 0$ for all t and the asserted upper bound in (26) follows. By replacing $\hat{C\!}_+$ by $\hat{C\!}_-:=\min \{0,\hat{C\!}\,\}$ one shows analogously the lower bound. $Ad\,$(ii) We derive from part (i) that ${\varTheta }\in I_\theta , \mathbf{S}\in I_S$ almost everywhere on $\varOmega \times [0,{T}]$. The assertion follows then with the fundamental assumption (19) on thermodynamics. $\square$

The next Lemma will be used in the proof of local existence and in the uniqueness part of Theorem 2.

Lemma 2

Let $({\varTheta }_1,\mathbf{S}_1)\in L^\infty _{{I_\theta }\times I_{S}}(\varOmega )$ and $({\varTheta }_2,\mathbf{S}_2)\in L^\infty _{{I_\theta }\times I_{S}}(\varOmega )$ be two pairs of admissible temperature and salinity fields.

• (i) Denote by $\rho _{rat}^{(1)}({\varTheta }_1,\mathbf{S}_1),\rho _{rat}^{(2)}({\varTheta }_2,\mathbf{S}_2)$ the two associated densities with respect to the equation of state (24). Then it holds almost everywhere in $\varOmega$

  $$\begin{aligned} \begin{aligned}&|\rho _{rat}^{(1)}-\rho _{rat}^{(2)}|\le K(|{\varTheta }_1-{\varTheta }_2|+|\mathbf{S}_1-\mathbf{S}_2|), \end{aligned}\end{aligned}$$

  (30)

  where $K=K(\theta _1,\theta _2,S_1,S_2,c_{i,j,k})$ is defined in (34).

• (ii) Denote by $\rho _{pol}^{(1)}({\varTheta }_1,\mathbf{S}_1),\rho _{pol}^{(2)}({\varTheta }_2,\mathbf{S}_2)$ the two associated densities with respect to the equation of state (25). Then it holds almost everywhere in $\varOmega$

  $$\begin{aligned} \begin{aligned}&|\rho _{pol}^{(1)}-\rho _{pol}^{(2)}|\le K(|{\varTheta }_1-{\varTheta }_2|+|\mathbf{S}_1-\mathbf{S}_2|), \end{aligned}\end{aligned}$$

  (31)

  where $K=K(\theta _1,\theta _2,S_1,S_2,d_{i,j,k})$ is defined in (34).

Proof

Adi) The difference of the two specific volumes $\upsilon _1-\upsilon _2$ reads as follows

$$\begin{aligned} \begin{aligned} \upsilon _1-\upsilon _2&=\sum _{i,j=0}^{N_i,N_j}c_{i,j}p_{stat}^k(\mathbf{S}^i_1 {\varTheta }^j_1-\mathbf{S}^i_2 {\varTheta }^j_2)\\&=\sum _{i,j=0}^{N_i,N_j}c_{i,j}(p_{stat}^k\big (\mathbf{S}^i_1( {\varTheta }^j_1-{\varTheta }^j_2) +{\varTheta }^j_2(\mathbf{S}^i_1-\mathbf{S}^i_2)\big ). \end{aligned}\end{aligned}$$

(32)

We write ${\varTheta }^j_1-{\varTheta }^j_2={\varTheta }_1({\varTheta }^{j-1}_1-{\varTheta }^{j-1}_2)+{\varTheta }^{j-1}_2({\varTheta }_1-{\varTheta }_2)$. Applying the same reformulation repeatedly to the first bracket with the $j-1$ terms creates a polynomial $\mathcal {P}_1^{(j-1)}({\varTheta }_1,{\varTheta }_2)$ in ${\varTheta }_1,{\varTheta }_2$ that factors out the difference $({\varTheta }_1-{\varTheta }_2)$ and that is of the form $\mathcal {P}_1^{(j-1)}({\varTheta }_1,{\varTheta }_2)=\sum _{l=1}^{j} {\varTheta }_1^{l-1}{\varTheta }2_2^{j-l}$. Finally we get

$$\begin{aligned} \begin{aligned} {\varTheta }^j_1-{\varTheta }^j_2=\mathcal {P}_1^{(j-1)}({\varTheta }_1,{\varTheta }_2)({\varTheta }_1-{\varTheta }_2). \end{aligned}\end{aligned}$$

(33)

After applying the same procedure to the salinity part in (32) and denoting by $\mathcal {P}_2^{(i-1)}(\mathbf{S}_1,\mathbf{S}_2)$ a polynomial of degree $i-1$ we obtain

$$\begin{aligned} \begin{aligned} |\upsilon _1-\upsilon _2|&=|\sum _{i,j=0}^{N_i,N_j}c_{i,j}p_{stat}^k\bigg (\mathbf{S}^i_1\mathcal {P}_1^{(j-1)}({\varTheta }_1,{\varTheta }_2)( {\varTheta }_1-{\varTheta }_2) +{\varTheta }^j_2\mathcal {P}_2^{(i-1)}(\mathbf{S}_1,\mathbf{S}_2)(\mathbf{S}_1-\mathbf{S}_2)\bigg )|\\&\le K(| {\varTheta }_1-{\varTheta }_2| +|\mathbf{S}_1-\mathbf{S}_2|),\\ \text {with }&K:=\sum _{i,j=0}^{N_i,N_j}\big \{|c_{i,j}|\, |p_{stat}^k|\big (\, |\mathbf{S}^i_1|\, |\mathcal {P}_1^{(j-1)}[{\varTheta }_1,{\varTheta }_2]| +|{\varTheta }^j_2|\, |\mathcal {P}_2^{(i-1)}[\mathbf{S}_1,\mathbf{S}_2]|\big )\big \}, \end{aligned}\end{aligned}$$

(34)

where ${\varTheta }_1,{\varTheta }_2$ and $\mathbf{S}_1,\mathbf{S}_2$ and the polynomials $\mathcal {P}_1^{(j-1)}({\varTheta }_1,{\varTheta }_2),\mathcal {P}_2^{(i-1)}(\mathbf{S}_1,\mathbf{S}_2)$ are bounded. From (34) follows the assertion of for density by using that the specific volumes are bounded below, according to the assumption (18). The assertion ii) follows analogously. $\square$

The next theorem gives the main result of this paper. In the case of a linear equation of state the original work is due to Cao and Titi [2] for the boundary conditions (5) and Kukavicas and Ziane [22] for the boundary conditions (6). The two theorems require different proofs but both rely on energy methods. It is reasonable to conjecture that the proof strategy via evolution equations methods as in [17] can also be extended to deal with a nonlinear equation of state.

Theorem 2

(Global Well-Posedness for Ocean Primitive Equations with Nonlinear Equation of State). Let a time interval $[0, \mathcal {T}]$, with $\mathcal {T} > 0$, be given. The boundary conditions for temperature and salinity are given by (4) and the boundary conditions for velocity are given either by (5) or by (6). Suppose the density $\rho$ is given either by the equation of state (24) such that $\rho =\rho _{rat}$ or by the equations of state (25) such that $\rho =\rho _{pol}$. The coefficients of the equation of state are determined such that (19) is satisfied. If the initial conditions for velocity, conservative temperature and absolute salinity satisfy $\mathbf{u}_0\in {V_{v}}$, ${\varTheta }_0\in {V_{\theta }}\cap L^\infty _{I_\theta }(\varOmega ),\mathbf{S}_0 \in {V_{S}}\cap L^\infty _{I_S}(\varOmega )$, then there exists on the time interval $[0, \mathcal {T} ]$ a unique strong solution $(\mathbf{u}, {\varTheta },\mathbf{S})$ of the system (1) in the sense of Definition 1. This solution depends continuously on the initial data with respect to the $C([0,\mathcal {T}], L^2(\varOmega ))$-topology.

Proof

We prove first the assertion for the equation of state given by the rational function (24) for both types of boundary conditions (5) and (6).

We observe that the first hypothesis on thermodynamics (18) is satisfied, because the initial conditions for temperature and salinity satisfy ${\varTheta }_0\in L^\infty _{I_\theta }(\varOmega ), \mathbf{S}_0 \in L^\infty _{I_S}(\varOmega )$, and Lemma 1 guarantees that ${\varTheta }(\cdot , t)\in L^\infty _{I_\theta }(\varOmega ), \mathbf{S}(\cdot ,t) \in L^\infty _{I_S}(\varOmega )$ for $t\in [0,\mathcal {T}]$. Since the second hypothesis on thermodynamics (19) is satisfied by assumption, it follows that the density remains in it physically admissible range, $\rho (\cdot , t)\in I_\rho$ for $t\in [0,\mathcal {T}]$.

We start with the velocity boundary conditions (5) considered by Cao and Titi [2]. Their proof contains an equation of state of the form $\rho =\alpha T$, with $\alpha =1$, where T denotes temperature. In this case we focus on the changes in the proof due to the nonlinear equation of state and stay otherwise close to the original proofs. Key estimates for velocity are recalled, for their derivation we refer to the original work [2].

Velocity boundary conditions (5): We use a Galerkin approximation of the system (1). For the velocity space ${V_{v}}$ we denote by $(\psi _m)_m$ the eigenfunctions of the Stokes operator, they constitute an orthonormal basis of of $\mathbf{L}^{2}(\varOmega )$. Define ${V_{v}}^m:=span\{\psi _k: k=1\ldots m\}$. The projection on ${V_{v}}^m$ is denoted by $\mathbf{P}_m$. The nature of the Stokes operator depends on the velocity boundary conditions. For boundary conditions (5) the Stokes operator is of Laplacian type.

Denote by $L_{\varTheta }:=\frac{1}{Rt_1}\triangle - \frac{1}{Rt_2}\partial _{zz}^2$ and $L_S:=-\frac{1}{Rs_1}\triangle - \frac{1}{Rs_2}\partial _{zz}^2$ the anisotropic Laplacian operators with domains $D(L_{\varTheta }), D(L_\mathbf{S})$, defined as the closure in ${{H^{1}}}$ of ${V_{\theta }}$ and ${V_{S}}$, respectively. The eigenfunctions of $L_{\varTheta }$ and $L_\mathbf{S}$ form orthonormal bases $(\phi _k^{\varTheta })_k$ and $(\phi _k^\mathbf{S})_k$ of $L^2(\varOmega )$ (cf. [1, 29]). Denote ${V_{\theta }}^m:=span\{\phi ^{\varTheta }_k: k=1\ldots m\}$ and ${V_{S}}^m:=span\{\phi ^\mathbf{S}_k: k=1\ldots m\}$. The projection on ${V_{\theta }}^m,{V_{S}}^m$ is denoted by $P_{{V_{\theta }}^m}, P_{{V_{S}}^m}$, respectively. We approximate velocity v, temperature and salinity ${\varTheta }, \mathbf{S}$ in terms of projection on the span of the respective basis functions

$$\begin{aligned} \begin{aligned}&{\mathbf{u}_m}(x,y,z,t):=\mathbf{P}_m \mathbf{u}(x,y,z,t):=\sum _{k=1}^m a_k(t)\psi _k(x,y,z),\\&{\varTheta }_m(x,y,z,t):=P_m{\varTheta }(x,y,z,t):=\sum _{k=1}^m d_k(t)\phi _k^{\varTheta }(x,y,z), \\&\mathbf{S}_m(x,y,z,t):=P_m\mathbf{S}(x,y,z,t):=\sum _{k=1}^m b_k(t)\phi _k^\mathbf{S}(x,y,z). \end{aligned} \end{aligned}$$

(35)

Step 1: Formulation of Approximative System. The Galerkin approximation for (1) reads as follows

$$\begin{aligned}&\partial _t {\mathbf{u}_m}+\mathbf{P}_m[({\mathbf{u}_m}\cdot \nabla ) {\mathbf{u}_m}] -\mathbf{P}_m\big [\big (\int _{-h}^z\nabla \cdot {\mathbf{u}_m}(\xi )\, d\xi \big )\partial _z {\mathbf{u}_m}\big ] +f\mathbf {k}\times {\mathbf{u}_m} \nonumber \\&\quad +\nabla \int _{-h}^zg\rho _m(x,y,z',t)\,dz' +\nabla (p_s)_m -\frac{1}{Re_1}{} \mathbf{P}_m\triangle \,{\mathbf{u}_m} - \frac{1}{Re_2}{} \mathbf{P}_m\partial ^2_{zz} {\mathbf{u}_m}=0, \end{aligned}$$

(36a)

$$\begin{aligned}&\partial _t {{\varTheta }_m}+ ({\mathbf{u}_m}\cdot \nabla ) {{\varTheta }_m} +P_m\big [\big (\int _{-h}^z\nabla \cdot {\mathbf{u}_m}(\xi )\, d\xi \big )\partial _z {{\varTheta }_m}\big ] +P_m\frac{1}{Rt_1}\triangle {{\varTheta }_m}\nonumber \\&\quad +P_m\frac{1}{Rt_2}\partial ^2_{zz}{{\varTheta }_m}=0, \end{aligned}$$

(36b)

$$\begin{aligned}&\partial _t {\mathbf{S}_m} + ({\mathbf{u}_m}\cdot \nabla ) {\mathbf{S}_m} +P_m\big [\big (\int _{-h}^z\nabla \cdot {\mathbf{u}_m}(\xi )\, d\xi \big )\partial _z {\mathbf{S}_m}\big ] +P_m\frac{1}{Rs_1}\triangle {\mathbf{S}_m}\nonumber \\&\quad +P_m\frac{1}{Rs_2}\partial ^2_{zz}\ {\mathbf{S}_m}=0, \end{aligned}$$

(36c)

here $\rho _m$ is defined in terms of the Galerkin approximation for ${{\varTheta }_m},{\mathbf{S}_m}$

$$\begin{aligned} \begin{aligned} \rho _m:=\rho ({{\varTheta }_m}, {\mathbf{S}_m}, p_{stat}), \end{aligned}\end{aligned}$$

(36d)

where the static pressure $p_{stat}$ is given by (22).

Step 2: Local existence. The Eq. (36) is an ordinary differential equation on the Hilbert space ${\mathcal {V}}^m:={V_{v}}^m\times {V_{\theta }}^m\times {V_{S}}^m$. Lipschitz continuity for the pressure gradient with density from the nonlinear equation of state follows from Lemma 2. For all other terms Lipschitz continuity is classical. Applying the Picard Theorem on Hilbert spaces it follows that for all $m\in \mathbb {N}$ the system (36) has a solution on a time interval $[0,T_m]$ whose length depends on m.

We continue with a priori estimates to establish global existence in time for the approximate system and then in the limit $m\rightarrow \infty$.

Step 3: $L^\infty ([0,\mathcal {T}], {{\mathcal {L}}^{2}}(\varOmega ))$- and $L^2([0,\mathcal {T}], {{\mathcal {H}}^{1}}(\varOmega ))$-Estimate for velocity and tracers. Let $C_m\in \{{{\varTheta }_m},{\mathbf{S}_m}\}$. Taking the $L^2$ scalar product of the tracer equation (1d) or (1e) with $C_m$ yields after integration-by-parts and with the tracer boundary condition (4)

$$\begin{aligned} \begin{aligned}&\frac{1}{2}d_t||C_m||^{2}_{{L^{2}}} + \frac{1}{Rc}\int _\varOmega |\nabla _3 C_m|^2\, dxdydz +k_C||C_m(z=0)||_{{L^{2}}}^2=0, \end{aligned}\end{aligned}$$

(37)

where $Rc=Rt$ if $C=\theta$ and $Rc=Rs$ if $C=S$. In (37) the integral with the tracer advection vanishes due to the incompressibility of $\mathbf{v}$. The boundary term vanishes in the case of salinity. From (37) follows

$$\begin{aligned} \begin{aligned}&d_t||C_m||_{{L^{2}}}\le 0. \end{aligned}\end{aligned}$$

(38)

With the Gronwall inequality and with $||C_m(t=0)||_{2}\le ||C(t=0)||_{2}$ for $m\in {\mathbb {N}}$, we conclude that

$$\begin{aligned} \begin{aligned} ||C_m(t)||_{{L^{2}}} \le ||C(t=0)||_{{L^{2}}} =:K_1 \end{aligned}\end{aligned}$$

(39)

where $K_1>0$ is independent of m. From (37) follows after integration with respect to time, with Lemma 1 and with $||C_m(t=0)||_{\infty }\le ||C(t=0)||_{\infty }$ that

$$\begin{aligned} \begin{aligned}&\int _0^t\bigg [\int _\varOmega |\nabla _3 C_m|^2\, dxdydz +k_C||C_m(z=0)||_{{L^{2}}}^{2}\bigg ]ds \le |\varOmega |||C(t=0)||_\infty t =:K_2(t), \end{aligned}\end{aligned}$$

(40)

with $K_2(t)$ bounded on $[0,\mathcal {T}]$ and independent on m. The estimate for velocity can be identically derived as in [2] (see pp 253-254). By taking the $L^2$-scalar product of (36) with ${\mathbf{u}_m}$ we obtain after integration by parts and by the fact that the Coriolis force acts perpendicular to the velocity

$$\begin{aligned} \begin{aligned} \frac{1}{2}\frac{d ||{\mathbf{u}_m}||_{{L^{2}}}^2}{d t}&+\int _\varOmega {\mathbf{u}_m}\cdot \bigg (\nabla \int _{-h}^z g\rho _m(x,y,\xi ,t)\, d\xi \bigg )dxdydz\\&+ \frac{1}{Re_1}||\nabla {\mathbf{u}_m}||^2_{{L^{2}}}+\frac{1}{Re_2}||\partial _z {\mathbf{u}_m}||^2_{{L^{2}}}=0. \end{aligned}\end{aligned}$$

(41)

The boundary term involving the surface pressure vanishes

$$\begin{aligned} \begin{aligned}&\int _\varOmega \nabla p_{s,m}\cdot {\mathbf{u}_m}\, dxdydz =\int _{\partial \varOmega } p_{s,m} \mathbf{n}\cdot {\mathbf{u}_m}\, ds -\int _\varOmega p_{s,m} \nabla \cdot {\mathbf{u}_m}\, dxdydz\\&\quad =\int _M p_{s,m}(\int ^0_{-h} \nabla \cdot {\mathbf{u}_m}\, dz)dxdy=\int _M p_{s,m}(w(0)-w(-h))\, dxdy=0. \end{aligned}\end{aligned}$$

(42)

With the inequality of Cauchy-Schwarz follows after integration by parts

$$\begin{aligned} \begin{aligned}&\frac{1}{2}\frac{d ||{\mathbf{u}_m}||_{{L^{2}}}^2}{d t} +\frac{1}{Re_1}||\nabla {\mathbf{u}_m}||_{{L^{2}}}^2 +\frac{1}{Re_2}||\partial _z {\mathbf{u}_m}||_{{L^{2}}}^2\\&\quad =-\int _\varOmega \nabla \cdot {\mathbf{u}_m}\bigg (\int _{-h}^z g\rho _m(x,y,\xi ,t)\, d\xi \bigg )dxdydz \le gh||\rho _m||_{{L^{2}}}||\nabla {\mathbf{u}_m}||_{{L^{2}}}, \end{aligned}\end{aligned}$$

(43)

and by the inequalities of Young and Poincaré

$$\begin{aligned} \begin{aligned} \frac{d ||{\mathbf{u}_m}||_{{L^{2}}}^2}{d t} + \frac{1}{Re_1C_{M}}|| {\mathbf{u}_m}||_{{L^{2}}}^2 \le (gh)^2Re_1||\rho _m||_{{L^{2}}}^2. \end{aligned}\end{aligned}$$

(44)

With Gronwall’s inequality we obtain

$$\begin{aligned} \begin{aligned} ||{\mathbf{u}_m}(t)||_{\mathbf{L}^{2}}^2&\le ||\mathbf{u}_m(t=0)||_{\mathbf{L}^{2}}^2e^{-\frac{1}{C_{M}Re_1}t}+(gh)^2C_{M}Re_1\int _0^t||\rho _m(s)||_{\mathbf{L}^{2}}^2ds\\&\le ||\mathbf{u}_m(t=0)||_{\mathbf{L}^{2}}^2 +(gh)^2C_{M}Re_1 K t, \end{aligned}\end{aligned}$$

(45)

where we have used that, due to the boundedness of density (cf. Lemma 1 (ii)), a constant K>0 exists such that $||\rho _m||_{L^\infty }^2\le K$. Integrating (45) with respect to time over [0,t] yields

$$\begin{aligned} \begin{aligned}&\int _0^t \frac{1}{Re_1}||\nabla {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2+\frac{1}{Re_2}||\partial _z {\mathbf{u}_m}||_{\mathbf{L}^{2}}^2ds \le ||{\mathbf{u}_m}(t=0)||_{\mathbf{L}^{2}}^2+ (gh)^2C_{M}Re_1Kt. \end{aligned}\end{aligned}$$

(46)

In summary we obtain from (39), (40) and from (45), (46)

$$\begin{aligned} \begin{aligned}&||{\mathbf{u}_m}(t)||_{\mathbf{L}^{2}}^2 +\int _0^t \frac{1}{Re_1}||\nabla {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2+\frac{1}{Re_2}||\partial _z {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2ds\\&\qquad + ||{{\varTheta }_m}(t)||_{\mathbf{L}^{2}}^2 +\mu \int _0^t\big [||\nabla _3{{\varTheta }_m}||_{\mathbf{L}^{2}}^2+k_\theta ||{{\varTheta }_m}(z=0)||_{2}^{2}\big ]ds\\&\qquad + ||{\mathbf{S}_m}(t)||_{{L^{2}}}^2 +\int _0^t\big [||\nabla _3 {\mathbf{S}_m}||_{{L^{2}}}^2 \big ]ds\\&\quad \le ||v(t=0)||_{\mathbf{L}^{2}}^2e^{-\frac{1 }{C_{M}Re_1}t} +(gh)^2C_{M}Re_1 Kt +h||v_0||_{\mathbf{L}^{2}}^2 =:K_3(t), \end{aligned}\end{aligned}$$

(47)

where $K_3(t)$ is bounded uniformly in m on $[0,\mathcal {T}]$. This shows that for all m the approximate solutions $({\mathbf{u}_m},{{\varTheta }_m},{\mathbf{S}_m})$ exist on the time interval $[0,\mathcal {T}]$.

Step 4: $L^\infty ([0,\mathcal {T}], {{\mathcal {H}}^{1}}(\varOmega ))$-Estimate for velocity and tracers. The proof of the $L^\infty ((0,T),H^1(\varOmega ))$-estimates for velocity follows [2] (cf. Sect. 3.3.1–3.3.3). This is the crucial part of the proof. It requires little changes due to the nonlinear equation of state. We sketch the steps only and refer the reader for more details to the original paper [2].

To derive the $H^1$-estimate for velocity, the velocity equation is vertically averaged, the vertical average is referred to as “barotropic” part and it is subtracted from the original velocity equation to obtain the “baroclinic” part. Both systems of equations contain as coupling terms the respective other velocity field (see Eqs. (36)–(37) and (38)–(39) in [2]). A major advantage of this decomposition is that the three dimensional baroclinic equations do not contain the surface pressure any more, the surface pressure is part of the barotropic equations, but these are two-dimensional. The barotropic system can be handled similarly as the 2D Navier–Stokes equations. The hydrostatic pressure that is related to the nonlinear equation of state is part of the barotropic as well as of the baroclinic equations.

One derives first a $L^\infty ([0,\mathcal {T}], L^6(\varOmega ))$-bound for the baroclinic velocity field ${\tilde{\mathbf{u}}}$ (see sect. 3.2. in [2]). Here the nonlinear equation of state within the hydrostatic pressure does not pose difficulties. The upper bound depends on the vertically averaged density ${\bar{\rho }}$ and its norms $||{\bar{\rho }}||_{{L^{2}}}, ||\nabla {\bar{\rho }}||_{{L^{2}}}$. The density itself is bounded. The estimation of the density gradient in terms of temperature and salinity gradients is shown below in (49). Then the $L^2$-scalar product of the barotropic equation with $\varDelta {\bar{\mathbf{u}}}$ is taken. As a consequence of the boundary conditions the integral with surface pressure vanishes. Due to the incompressibility of the barotropic velocity the term with the hydrostatic pressure gradient vanishes also and the presence of nonlinear thermodynamics has no consequences for this estimate (see sect. 3.3.1 in [2]). The upper bound depends on the previous estimate of the $L^6$-norm. The next step in Sect. 3.3.2 of in [2] is an a priori estimate on the vertical derivative $\partial _z\mathbf{u}$ in the $L^\infty ([0,\mathcal {T}], {L^{2}}(\varOmega ))$ norm. Here the hydrostatic pressure can be treated as in the case of a linear equation of state.

In the final step (sect. 3.3.3. of in [2]) the $L^2$-scalar product of the original velocity equation with the horizontal Laplace operator $\triangle v$ is taken. By using the important Proposition 2.2 in [1] one arrives at the estimate below (see p. 261 in [2])

$$\begin{aligned} \begin{aligned}&d_t||\nabla {\mathbf{u}_m}||_{\mathbf{L}^{2}}^2+\frac{1}{Re_1}||\triangle {\mathbf{u}_m}||_{\mathbf{L}^{2}}^2+\frac{1}{Re_2}||\nabla \partial _z {\mathbf{u}_m}||_{L^2}^2\\&\quad \le C\big (K_6(t)+K_z(t)||\nabla {\mathbf{u}_m} ||_{\mathbf{L}^{2}}^2\big )||\nabla {\mathbf{u}_m}||_{\mathbf{L}^{2}}^2 + C||\nabla \rho _m||_{L^2}^2, \end{aligned}\end{aligned}$$

(48)

where $K_6(t),K_z(t)$ are bounded functions on $[0,\mathcal {T}]$ (see Eqs. (67) and (76) in [2]).

As a consequence of the nonlinear equation of state the pressure gradient on the right-hand side has to be estimated. This is done with (24) as follows

$$\begin{aligned} \begin{aligned}&||\nabla \rho _m||_{L^2}^2 =\int _\varOmega |\frac{\nabla \upsilon _m}{|\upsilon _m|^2}|^2dxdydz \le c \int _\varOmega |\nabla \upsilon _m|^2dxdydz\\ {}&\quad =c\sum _{i,j,k=0}^{N_i,N_j,N_k}c_{i,j,k}\int _\varOmega |\nabla \big (\mathbf {S}^i_m {\varTheta }^j_m p_{stat}^k\big ) |^2dxdydz\\ {}&\quad \le C\sum _{i,j=0}^{N_i,N_j}\int _\varOmega |\mathbf {S}^i_m\nabla {\varTheta }^j_m|^2 dxdydz +C\sum _{i,j=0}^{N_i,N_j}\int _\varOmega |{\varTheta }^j_m\nabla \mathbf {S}^i _m|^2 dxdydz\\ {}&\quad \le C\sum _{i,j=0}^{N_i,N_j}\int _\varOmega |\mathbf {S}^i_m{\varTheta }^{j-1}_m |^2\, |\nabla {\varTheta }_m|^2 dxdydz +C\sum _{i,j=0}^{N_i,N_j}\int _\varOmega |{\varTheta }^j_m \mathbf {S}^{i-1}_m |^2\, |\nabla \mathbf {S}_m |^2 dxdydz\\ {}&\quad \le K (||\nabla {\varTheta }_m||_{L^2}^2+||\nabla \mathbf {S}_m||_{L^2}^2). \end{aligned}\end{aligned}$$

(49)

The $K>0$ on the right-hand side of (49) depends on the coefficients $c_{i,j,k}$ of the equation of state and the powers of ${\varTheta }_m,\mathbf{S}_m$. Using the estimate $|{\varTheta }^j_m\mathbf{S}^i_M|\le ||{\varTheta }||_{L^\infty }^j||\mathbf{S}||_{L^\infty }^i\le |I_{\varTheta }|^j\, |I_\mathbf{S}|^i$, where $|I_{\varTheta }|,\, |I_\mathbf{S}|$ denote the size of the physically admissible range of ${\varTheta }$ and $\mathbf{S}$, introduced in Definition 2, it follows that K is not only finite but also independent of m. Since ${\varTheta }_m, \mathbf{S}_m\in L^2([0,\mathcal {T}], {H^{1}}(\varOmega ))$ we conclude from (49) that $\rho _m\in L^2([0,\mathcal {T}], {H^{1}}(\varOmega )))$. Using this estimate we can complete the $L^\infty ([0,\mathcal {T}],\mathbf{H}^{1}(\varOmega ))$-estimate for velocity as in sect. 3.3.3 of [2] and obtain with the Gronwall inequality

$$\begin{aligned} \begin{aligned}&||\nabla {\mathbf{u}_m}(t)||_{\mathbf{L}^{2}}^2+||\partial _z {\mathbf{u}_m}(t)||_{\mathbf{L}^{2}}^2 +\int _0^t \frac{1}{Re_1}||\triangle {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2 +\frac{1}{Re_2}||\partial _{zz}^2 {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2\\&\quad +\frac{1}{Re_2}||\nabla \partial _z {\mathbf{u}_m}(s)||_{\mathbf{L}^{2}}^2\, ds \le K_4(t), \end{aligned}\end{aligned}$$

(50)

where $K_4(t)$ is bounded on $[0,\mathcal {T}]$, uniform in m. The upper bound $K_4(t)$ contains the $L^2([0,\mathcal {T}], {H^{1}}(\varOmega )$-norms of ${\mathbf{u}_m},{{\varTheta }_m},{\mathbf{S}_m}$ but these are, according to (47), uniformly bounded with respect to m.

The $L^\infty ([0,\mathcal {T}],{H^{1}}(\varOmega )$-bound for temperature and salinity that is proven in [2] does not change, with $C_m\in \{ {\varTheta }_m, \mathbf{S}_m\}$ it reads as follows (see (79) in [2])

$$\begin{aligned} \begin{aligned}&||\nabla C_m(t)||_{{L^{2}}}^2+||\partial _z C_m(t)||_{{L^{2}}}^2 +k_C||\nabla C_m(z=0)||_{{L^{2}}}^2 +\int _0^t\bigg [ \frac{1}{Rc_1}||\triangle C_m(s)||_{\mathbf{L}^{2}}^2\\&\quad +\big (\frac{1}{Rc_1}+\frac{1}{Rc_2})||\nabla \partial _z C_m(s)||_{\mathbf{L}^{2}}^2+k_C||\nabla C_m(z=0)||_{{L^{2}}}^2\big )\frac{1}{Rc_2}||\partial _{zz}^2 C_m(s)||_{\mathbf{L}^{2}}^2\bigg ] \, ds\\&\quad \le K_C(t), \end{aligned}\end{aligned}$$

(51)

where the boundary term on the left-hand side is absent in case $C_m=\mathbf{S}_m$ and where $K_C$ on the right-hand side is bounded on $[0,\mathcal {T}]$.

Step 5: Estimate on the time derivative. We have to show that $\partial _t {\mathbf{u}_m}$ is bounded uniformly in the $L^2([0,\mathcal {T}], \mathbf{L}^{2}(\varOmega ))$-norm. For this purpose one takes the $L^2$.scalar product of the velocity equation with $\partial _t {\mathbf{u}_m}$. For the pressure gradient follows with (49)

$$\begin{aligned} \begin{aligned} \int _\varOmega \partial _t {\mathbf{u}_m}\nabla \int _z^0\rho _m(x,y,z')dz'\, dxdydz&\le c||\partial _t {\mathbf{u}_m}||_{\mathbf{L}^{2}}||\nabla \rho _m||_{{L^{2}}}\\&\le c||\partial _t {\mathbf{u}_m}||_{\mathbf{L}^{2}}(||\nabla {{\varTheta }_m}||_{{L^{2}}}+||\nabla {\mathbf{S}_m}||_{{L^{2}}})\\&\le ||\partial _t {\mathbf{u}_m}||_{\mathbf{L}^{2}}(K_{\varTheta }(t)+K_\mathbf{S}(t)), \end{aligned}\end{aligned}$$

(52)

where $K_{\varTheta }, K_\mathbf{S}$ are the bounded functions from (51) for temperature and salinity respectively. With the Young inequality the time derivative on the right-hand side of (52) can be compensated with the left hand side. Together with the corresponding bounds of the scalar products between $\partial _t {\mathbf{u}_m}$ and the remaining terms in the velocity equation it follows that $(\partial _t {\mathbf{u}_m})_m$ is uniformly bounded in $L^2([0,\mathcal {T}], \mathbf{L}^{2}(\varOmega ))$. The same results applies to $(\partial _t{{\varTheta }_m})_m,(\partial _t {\mathbf{S}_m})_m$. Since $({\mathbf{u}_m})_m,({{\varTheta }_m})_m,({\mathbf{S}_m})_m$ are bounded in $L^2([0,\mathcal {T}], \mathbf{L}^{2}(\varOmega ))$ and $L^2([0,\mathcal {T}], {L^{2}}(\varOmega ))$ respectively, it follows from the Lions–Aubin compactness lemma [9] that a subsequence converges in $L^2([0,\mathcal {T}], {{\mathcal {H}}^{1}}(\varOmega ))$ to a limit $(\mathbf{u},{\varTheta },\mathbf{S})$.

Step 6: Global Existence. This is proven by contradiction. Assuming that there exists a maximal time interval of existence that is finite, $\mathcal {T}^*<\infty$. This implies for the $H^1$-norms $\lim \sup _{t\rightarrow \mathcal {T}^*_-}||\mathbf{u}(t)||_{\mathbf{H}^{1}}+||{\varTheta }(t)||_{{H^{1}}}+||\mathbf{S}(t)||_{{H^{1}}}=+\infty$. But this contradicts the a priori estimate above. Therefore $\mathcal {T}^*=\infty$ and the strong solution exists globally in time.

Step 7: Uniqueness For the uniqueness consider two strong solutions $(\mathbf{u}_1,{\varTheta }_1,\mathbf{S}_1)$ and $(\mathbf{u}_2,{\varTheta }_2,\mathbf{S}_2)$, and define the differences

$$\begin{aligned} \begin{aligned}&\hat{\mathbf {v}}:=\mathbf {u}_1-\mathbf {u}_2,\ {\hat{{\varTheta }}}:={\varTheta }_1-{\varTheta }_2,\ \hat{S}:=\mathbf {S}_1-\mathbf {S}_2,\\ \text{ and }&\hat{\rho }:=\rho ({\varTheta }_1,\mathbf {S}_1,p_{stat})-\rho ({\varTheta }_2,\mathbf {S}_2,p_{stat}). \end{aligned}\end{aligned}$$

(53)

Now we take the $L^2$ scalar product of $(\hat{\mathbf {v}},{\hat{{\varTheta }}}, \hat{S})$ with the equations that govern the evolution of this differences (see (81), 82) in [2]). For the velocity difference equation we have

$$\begin{aligned}\begin{aligned}&\frac{1}{2}\frac{d}{d t}|| \hat{\mathbf {v}}||_{\mathbf {L}^{2}}^2+\frac{1}{Re_1}||\nabla \hat{ \mathbf {v}}||_{\mathbf {L}^{2}}^2 +\frac{1}{Re_2}||\partial _z \hat{\mathbf {v}}||_{\mathbf {L}^{2}}^2\\ {}&\quad = -\int _\varOmega \bigg [(\mathbf {u}_1\cdot \nabla ) \hat{\mathbf {v}}+(\hat{\mathbf {v}}\cdot \nabla ) \mathbf {u}_2 -\big (\int _{-h}^z\nabla \cdot \mathbf {u}_1(x,y,\xi ,t)\, d\xi \big )\partial _z \hat{\mathbf {v}}\bigg ]\cdot \hat{\mathbf {v}}\, dxdydz\\ {}&\qquad -\bigg [\int _{-h}^z\nabla \cdot \hat{\mathbf {v}}(x,y,\xi ,t)\, d\xi \big )\partial _z \mathbf {u}_2\cdot \hat{\mathbf {v}}\, dxdydz - \int _\varOmega \nabla \big (\int _{-h}^z \hat{\rho }(x,y,\xi ,t)d\xi \bigg ]\cdot \hat{ \mathbf {v}}\, dxdydz =0. \end{aligned}\end{aligned}$$

The only term that is different from the uniqueness proof of [2] is the term involving the difference of the pressure gradients. Using the boundedness of density and Lemma 2 it can be estimated as follows

$$\begin{aligned}\begin{aligned}&|\int _\varOmega \big (\nabla \int _{-h}^z\hat{\rho }(x,y,\xi )\, d\xi \big )\hat{ \mathbf {v}}\, dxdydz| =|\int _\varOmega \big (\int _{-h}^z\hat{\rho }(x,y,\xi )\, d\xi \big )\nabla \cdot \hat{\mathbf {v}}\, dxdydz|\\ {}&\quad \le c||\hat{\rho }||_{{L^{2}}}||\nabla \hat{\mathbf {v}}||_{\mathbf {L}^{2}}\le c||\upsilon _1-\upsilon _2||_{{L^{2}}}||\nabla \hat{\mathbf {v}}||_{\mathbf {L}^{2}} \le c(||\hat{\theta }||_{{L^{2}}}+||\hat{S}||_{{L^{2}}})||\nabla \hat{\mathbf {v}}||_{\mathbf {L}^{2}}, \end{aligned}\end{aligned}$$

where $\upsilon _1-\upsilon _2$ denotes the difference of the specific volumes. By carrying out the other estimates in the velocity difference equation exactly as in Sect. 4 of [2] and by applying the estimates for the temperature difference equation in [2] to temperature and salinity the we arrive at the following inequality

$$\begin{aligned} \begin{aligned}&\frac{1}{2}\frac{d}{d t}\big (|| \hat{\mathbf {v}}||_{\mathbf {L}^{2}}^2+||\hat{{\varTheta }}||_{{L^{2}}}^2+||\hat{\mathbf {S}}||_{{L^{2}}}^2\big )\\ {}&\quad \le C\big (||\nabla \mathbf {u}_2||_{\mathbf {L}^{2}}^4+||\nabla {\varTheta }_2||_{{L^{2}}}^2+||\nabla \mathbf {S}_2||_{{L^{2}}}^2 +||\partial _z \mathbf {u}_2||_{\mathbf {L}^{2}}^2||\nabla \partial _z \mathbf {u}_2||_{\mathbf {L}^{2}}^2 \\ {}&\qquad +||\partial _z{\varTheta }_2||_{{L^{2}}}^2||\nabla \partial _z{\varTheta }_2||_{{L^{2}}}^2 +||\partial _z\mathbf {S}_2||_{{L^{2}}}^2||\nabla \partial _z\mathbf {S}_2||_{{L^{2}}}^2\big ) \big (|| \hat{\mathbf {v}}||_{\mathbf {L}^{2}}^2+||\hat{{\varTheta }}||_{{L^{2}}}^2+||\hat{\mathbf {S}}||_{{L^{2}}}^2\big ). \end{aligned}\end{aligned}$$

(54)

From the Gronwall inequality follows

$$\begin{aligned} \begin{aligned}&|| \hat{\mathbf {v}}(t)||_{\mathbf {L}^{2}}^2+||\hat{{\varTheta }}(t)||_{{L^{2}}}^2+||\hat{\mathbf {S}}(t)||_{{L^{2}}}^2\\ {}&\quad \le (|| \hat{\mathbf {v}}(t=0)||_{\mathbf {L}^{2}}^2+||\hat{{\varTheta }}(t=0)||_{{L^{2}}}^2+||\hat{\mathbf {S}}(t=0)||_{{L^{2}}}^2\big ) exp K(t), \end{aligned}\end{aligned}$$

(55)

where K(t) is given by the integral with respect to time over the interval from initial time $t=0$ to time t of the first bracket on the right-hand side of (54). Since $(\mathbf{u}_2,{\varTheta }_2,\mathbf{S}_2)$ is by assumption a strong solution, the function K(t) is bounded. This proves the continuous dependency of the strong solution on the initial condition. It follows also from (55) that if the two strong solutions $(\mathbf{u}_1,{\varTheta }_1,\mathbf{S}_1), (\mathbf{u}_2,{\varTheta }_2,\mathbf{S}_2)$ have the same initial conditions then they coincide for all times, i.e. strong solutions are unique.

Velocity boundary conditions (6): This case that was investigated by Kukavica and Ziane [22] and we follow their proof idea. The original proof in [22] proceeds in two steps. First, the primitive equations are considered in the following form

$$\begin{aligned}&\partial _t \mathbf{u}+(\mathbf{u}\cdot \nabla ) \mathbf{u}+w\partial _z \mathbf{u}+\nabla p_s -\frac{1}{Re_1}\triangle \,\mathbf{u}- \frac{1}{Re_2}\partial ^2_{zz} \mathbf{u}=F, \end{aligned}$$

(56)

$$\begin{aligned}&\nabla \cdot \mathbf{u}+\partial _zw=0, \end{aligned}$$

(57)

The hydrostatic pressure $p_h$ and consequentially the whole thermodynamics are absent in this equation (cf. (3)) and only the surface pressure is present. In other words, the authors assume $\rho \equiv 0$. Also the Coriolis force is not considered in (56). For the boundary conditions (6) the Stokes operator is related to a Stokes problem. This requires a special treatment of the surface pressure by using results from the regularity theory of the Stokes problem. For details we refer to [22]. For a square-integrable right-hand side $F\in L^2([0,\mathcal {T}], {L^{2}}(\varOmega ))$ the authors show the existence and uniqueness of solution to (56), (57) for initial conditions $\mathbf{u}_0\in {H^{1}}$.

The second step of the proof includes linear thermodynamics by specifying the right-hand side F to include the hydrostatic pressure and the Coriolis force

$$\begin{aligned} \begin{aligned} F(x,y,z,t):&= g\int ^0_z\nabla \rho (x,y,\xi ,t)\, d\xi -f\mathbf {k}\times \mathbf{u}\\ \text {with }&\nabla \rho :=\alpha \nabla {\varTheta }-\beta \nabla \mathbf{S}. \end{aligned}\end{aligned}$$

(58)

The assertion in [22] follows now by showing that F in (58) satisfies $F\in L^2([0,\mathcal {T}], {L^{2}}(\varOmega )$. This follows from $({\varTheta },\mathbf{S})\in L^2([0,\mathcal {T}], {H^{1}}(\varOmega )$ for solutions of (1d), (1e) with initial conditions ${\varTheta }_0,\mathbf{S}_0\in {H^{1}}$ (see above, (40)).

For the case of a nonlinear equation of state considered here we adopt this strategy. We have to show that $F\in L^2([0,\mathcal {T}], {L^{2}}(\varOmega )$ with F as in (58) but now with density given by (24) in Definition 2. From (49) follows that $\nabla \rho \in L^2([0,\mathcal {T}], {L^{2}}(\varOmega ))$, because ${\varTheta },\mathbf{S}\in L^2([0,\mathcal {T}], {H^{1}}(\varOmega ))$. This concludes the proof for the case $\rho =\rho _{rat}$.

Let the equation of state be given by the polynomial (25) such that $\rho =\rho _{pol}$. The proof is analogous to the case $\rho =\rho _{rat}$, the only difference is the estimation of the density gradient (49), for which holds

$$\begin{aligned}\begin{aligned}&||\nabla \rho _m||_{L^2}^2 =\int _\varOmega |\nabla \big ( \sum _{i,j,k=0}^{M_i,M_j,M_k}d_{i,j,k}\mathbf{S}^i_m {\varTheta }^j_m p_{stat}^k\big ) |^2dxdydz\\&\quad \le K\sum _{i,=0}^{M_i,M_j}\int _\varOmega |\nabla \big (\mathbf{S}^i_m {\varTheta }^j_m\big ) |^2dxdydz \le K (||\nabla {\varTheta }_m||_{L^2}^2+||\nabla \mathbf{S}_m||_{L^2}^2). \end{aligned}\end{aligned}$$

The intermediate steps of this estimate are the same as in (49). The proof can now be concluded for both types of velocity boundary conditions in the same way as for the equation of state given by the rational function $\rho =\rho _{rat}$. $\square$

We state a generalization of Theorem 2 in which we use equations of state that are rational functions and where we replace the hypothesis on thermodynamics (19) by a condition on the denominator of the rational function, while hypothesis (18) remains valid. More specifically, for $T\in L^\infty _{I_\theta }(\varOmega ),S\in L^\infty _{I_S}(\varOmega )$ and static pressure $p_{stat}$ given by (22), let $\rho =\rho (T,S,p_{stat})$ be an equation of state that is of the form

$$\begin{aligned} \begin{aligned} \rho (T,S,p_{stat})=\frac{Q_1(T,S,p_{stat})}{Q_2(T,S,p_{stat})}, \end{aligned}\end{aligned}$$

(59)

where $Q_1,Q_2$ are polynomials in $(T,S,p_{stat})$ of arbitrary order. Assume that for $T\in L^\infty _{I_\theta }(\varOmega ),S\in L^\infty _{I_S}(\varOmega )$ the function $Q_2(x,y,z,t):=Q_2(T(x,y,z,t),S(x,y,z,t), p_{stat})$ is bounded below, i.e. $|Q_2(x,y,z,t)|\ge c>0$ on $\varOmega \times [0,\mathcal {T}]$. Note that under these conditions the density $\rho$ is bounded, but not necessarily within the physically admissible range. This would require an additional assumption on $\rho$ or $Q_1$ but for the well-posedness theorem the boundedness of density is sufficient.

Corollary 1

(Well-Posedness for Alternative Equations of State) Let T, S be generic temperature and salinity variable that satisfy advection–diffusion equations (1d), (1e) with boundary conditions (4). Let the equation of state given by (59). If the initial conditions for velocity, temperature and salinity satisfy $\mathbf{u}_0\in {V_{v}}$, $T_0\in {V_{\theta }}\cap L^\infty _{I_\theta }(\varOmega )$, $S_0 \in {V_{S}}\cap L^\infty _{I_S}(\varOmega )$, then the statement of Theorem 2 remains valid.

Proof

The assumptions of the advection–diffusion equations for the temperature and salinity variables T and S imply with Lemma 1 that $T(t)\in L^\infty _{I_\theta }(\varOmega ), S(t)\in L^\infty _{I_S}(\varOmega )$ for $t\in [0,\mathcal {T}]$. The assumed properties of the equation of state (59) guarantees that the density satisfies $\rho \in L^\infty (\varOmega )$. The estimation of the density gradient (cf. (49) follows analogously by using the lower bound on $Q_2$. The proof can be carried out along the lines of the proof of Theorem 2. $\square$