1 Introduction

The oceans have larger heat capacity than the atmosphere and are thought to play a central role in stabilizing the Earth’s climate system. It has recently been determined that the deep ocean (> 2000 m depth) warmed significantly from the 1990s to the 2000s, accounting for ~ 30% of the Earth’s energy change inventory (IPCC 2014). This highlights the importance of understanding the mechanisms that cause temporal variations in ocean temperature through the entire water column, including the deep ocean, when assessing the present climate state. However, it is difficult to investigate such mechanisms using only the limited observational data available, so ocean general circulation models (OGCMs) must be utilized. Turbulent vertical mixing, which needs to be parameterized in OGCMs, is a key element in simulating the deep-ocean state, as it is one of the most important physical processes controlling the global meridional overturning circulation (MOC) (Munk and Wunsch 1998).

Deep-water formation does not occur in the Pacific Ocean, which thus corresponds to the endpoint of the global thermohaline circulation. This suggests that vertical mixing plays an important role in the balance determining the deep-water properties in this region. In addition, modern, high-quality, full-depth, hydrographic section observations repeated at decadal intervals have revealed modest, but basin-wide warming of the deep water in the Pacific Ocean (Fukasawa et al. 2004; Kawano et al. 2006). When integrated over each basin, the increase in heat content in the Pacific is larger than in other basins, with the exception of the Southern Ocean (Kouketsu et al. 2011). The particularly large increase in the Southern Ocean is probably related to the deep-water formation in that region.

Vertical mixing below the surface mixed layer in the open ocean, which we refer to as “interior mixing”, is driven by a wide range of processes with distinct governing physics in the ocean, and is largely related to the breaking of internal waves, the main sources of which are tides and winds (MacKinnon et al. 2013). Barotropic tidal currents generate internal waves when flowing over topographic obstacles. These waves are called “internal tides”. The global energy conversion rate from barotropic to internal tides was estimated at ~ 0.9 TW on the basis of altimetric observations (Egbert and Ray 2000) and numerical experiments (Niwa and Hibiya 2001; Simmons et al. 2004a). Two types of mixing related to these waves are recognized: near-field mixing (NFM) caused by waves breaking near their generation site owing to direct breaking and enhanced rates of interaction with other internal waves, and far-field mixing (FFM) associated with low-mode waves propagating over greater distances and causing vertical mixing through various mechanisms, such as a cascade to smaller scales through wave–wave interactions. Wind stress excites inertial motions in the surface mixed layer, and a part of their energy can radiate away downward and equatorward in the form of propagating near-inertial internal waves (Alford et al. 2016). The global wind energy input is estimated to be of the same order as the input by tides, but the energy propagating downward to the deep ocean is limited (e.g. Furuichi et al. 2008). Thus, we focus on tidally induced mixing in this study, while the effect of the wind is implicitly included in the background mixing.

Various schemes have been proposed for parameterizing vertical mixing on the basis of knowledge about mixing derived from observations, theories, and numerical simulations. Enhanced vertical mixing over rough topography is confirmed by numerous direct observations of turbulence and indirect estimates using fine-scale observations (e.g. Polzin et al. 1997). In addition, the distribution of the rate of conversion from barotropic to internal tidal energy can be estimated using tidal models and satellite observations of sea surface height (e.g. Jayne and St. Laurent 2001); this has enabled the development of energetically feasible NFM parameterization schemes, which were first proposed by St. Laurent et al. (2002) and have since been modified in various ways. These types of NFM schemes are now widely used in OGCMs and climate models (e.g. Jayne 2009). However, OGCMs using the NFM schemes do not necessarily simulate the MOC realistically, particularly over the Pacific Ocean (Simmons et al. 2004b; Hasumi et al. 2010).

It is still unclear where and how low-mode waves cause vertical mixing (MacKinnon et al. 2017). Therefore, FFM has not been explicitly parameterized in OGCMs, but is instead generally accounted for through an ad hoc background mixing. Oka and Niwa (2013) incorporated an FFM scheme in an OGCM for the first time by extending the approach for the NFM scheme. They demonstrated that the simulated MOC and distribution of age tracers in the Pacific Ocean are improved by incorporating FFM in addition to NFM. Their FFM scheme is based on dissipation rates of baroclinic tidal energy that are diagnosed in a three-dimensional tide model, and assumes that the energy dissipation is distributed uniformly in the vertical. As they suggest, however, this vertical profile is merely a crude expression of the wave dissipation processes, and should be improved. In addition, there is room for improvement in the horizontal distribution. Hibiya et al. (2006) showed that the supplied semidiurnal internal tide energy is not efficiently cascaded down to dissipation scales at latitudes over 35°, on the basis of the diapycnal diffusivity in the open ocean as estimated from observed fine-scale vertical shear. This distribution suggests that the cascade of the supplied semidiurnal internal tide energy to dissipation scales is dominated by wave–wave interaction known as parametric subharmonic instability (Staquet and Sommeria 2002), which numerical models predict to occur only equatorward of 30° latitude (Hibiya et al. 1996, 1998, 2002; Furuichi et al. 2005).

Other attempts to estimate vertical mixing have involved inverse estimate approaches that do not use any information about vertical mixing. Using an one-dimensional advection–diffusion model and a vertical temperature profile, Munk (1966) estimated a vertical diffusivity of ~ \(1 \times {10}^{-4}\) m2 s−1 in the Pacific Ocean. This value is roughly consistent with estimates based on observed broad-scale hydrographic characteristics derived using inverse box models (e.g. Ganachaud 2003; Lumpkin and Speer 2007), but is an order of magnitude higher than the direct observation of interior vertical mixing in the open ocean (e.g. Ledwell et al. 1998). It is commonly accepted that at least part of this discrepancy may be explained by the neglect of physical processes in the simplified models (e.g. adiabatic upwelling of North Atlantic deep water in the Southern Ocean; Toggweiler and Samuels 1998; Webb and Suginohara 2001). This illustrates the possible confusion introduced by the assumptions built into the simple model, and motivates the use of a more realistic model.

The inverse approach using OGCMs and data assimilation techniques is a promising means of estimating mixing parameters that brings OGCMs into consistency with the observed hydrography (e.g. Stammer 2005). Recent data synthesis experiments have estimated vertical and other mixing coefficients on each model grid cell through the four-dimensional variational (4D-VAR) assimilation approach, and show that adjusting mixing coefficients is effective in reducing a global misfit between a model simulation and ocean observations over recent decades (Liu et al. 2012; Forget et al. 2015). In those experiments, however, the adjustment of the vertical mixing coefficient is small, especially in the deep ocean. The distribution obtained in previous studies that lacks the bottom-intensified structure is inconsistent, at least with current understanding of NFM. This must be partly related to the integration periods of the models, which are far shorter than the time scale of the MOC and/or the residence time of the deep water.

When the number of unknown parameters is small, the Green’s function method, which is a simpler data-assimilation technique, is also an effective tool for correcting model biases (Menemenlis et al. 2005). Toyoda et al. (2015) reduced the number of unknown parameters by representing the interior mixing in an OGCM by a linear combination of three semi-empirical schemes, and estimated optimal values for the weighting factor for each scheme and other physical parameters by conducting a linear optimization experiment using the Green’s function method. They successfully reproduced the climatological state of temperature and salinity in the deep Pacific Ocean in an OGCM. The three mixing schemes are (a) a horizontally uniform background profile (Tsujino et al. 2000), (b) a parameterization depending on the local static stability (Gargett, 1984), and (c) a parameterization for bottom intensified mixing depending on the bottom topography (Hasumi and Suginohara 1999). The optimized vertical mixing was largely represented by the sum of (a) and (c), suggesting that broadly distributed mixing is needed in addition to the bottom intensified mixing above the rough topography. This is roughly consistent with Oka and Niwa (2013), who implemented a broadly distributed FFM in addition to the NFM above rough topography. Note that the two semi-empirical schemes are not fully consistent with the current knowledge of tidally induced mixing, such as the energetic constraints estimated by the tide models, thus it is not clear whether the tidally induced mixing schemes work as effectively as the semi-empirical schemes. However, this is an important example suggesting that we could synthesize knowledge of mixing and modeling of the large-scale ocean state through refinement of the existing vertical mixing schemes based on hydrographic observations within the framework of OGCMs. This motivates us to follow the approach of Toyoda et al. and incorporate the current knowledge of tidally induced mixing.

In this study, we revise the Toyoda et al. (2015) model by replacing the interior mixing schemes from the semi-empirical schemes with two tidally induced schemes for NFM and FFM that will be described in detail in Sect. 2.2, and conduct an optimization experiment to assess the usefulness of the combination of tidally induced schemes and discuss the roles of NFM and FFM in maintaining the deep ocean state in the Pacific.

2 Method

2.1 OGCM

As a base model, we used Version 3 of the Modular Ocean Model (Pacanowski and Griffies 2000), modified by Toyoda et al. (2015). It involves established and effective parameterizations, including a turbulent closure scheme for the surface mixed layer (Noh 2004), an isopycnal mixing scheme (Gent and McWilliams 1990), and parameterizations for double-diffusion due to salt fingering (Schmitt 1988) and diffusive convection (Marmorino and Caldwell 1976). To represent convection, vertical diffusivity was set to 0.1 m2 s−1 where stratification is unstable. In addition to implementing the two vertical diffusivity schemes for interior mixing (Sect. 2.2), we added geothermal heat flux based on the global distribution map of Davies (2013).

The quasi-global model domain covers the region between 75° S and 80° N. The horizontal resolution is 1° for both longitude and latitude, and there are 45 levels in the vertical with thicknesses of 10–400 m, together with an additional layer representing the bottom boundary layer. The initial condition is the climatological state in the optimized model of Toyoda et al. (2015). The model is forced by 10-day climatological air–sea fluxes (heat, freshwater, and momentum) derived from the 6-hourly National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) ‘Reanalysis 1’, 1957–2014. We also applied a commonly used flux-correction method by relaxing the sea surface temperature (SST) and salinity (SSS) to monthly climatologies on a timescale of 30 days. The SST climatology was derived from a historical 10-day-mean SST dataset for 1957–2014 compiled from the Optimally Interpolated SST (OISST) dataset (Reynolds et al. 2002), extended back using the reconstructed Reynolds SST dataset. The SSS climatology is from the World Ocean Atlas 1998 (WOA98; see Antonov et al. 1998). Potential temperature (T) and salinity (S) were restored to annual mean values in WOA98 through the whole water column, on the same timescale, along the northern and southern open boundaries and in regions where the model cannot simulate important processes, such as deep convection and inflow from the unresolved marginal seas; these include the northern polar regions north of 74° N, the continental shelf and slope shallower than 2000 m around Antarctica south of 69° S, the Sea of Okhotsk, the area off the coast of Greenland, and the entrances to the Mediterranean Sea, Red Sea, and Persian Gulf. The model was integrated for 2900 years using the acceleration method (Bryan 1984), and thereafter for 1160 years without acceleration (total 4060 years) to obtain the current climatological state.

2.2 Vertical diffusivity schemes

Two tidally induced vertical diffusivity schemes were implemented, rather than the three semi-empirical schemes that Toyoda et al. (2015) used for interior mixing: the NFM scheme of St. Laurent et al. (2002) and the FFM scheme of Hibiya et al. (2006), as described below.

Near-field mixing

As described in the introduction, the NFM scheme first proposed by St. Laurent et al. (2002) is widely used in OGCMs. Following their scheme, we represent the turbulent dissipation rate, \(\varepsilon\), as:

$$\begin{array}{c}\varepsilon =\frac{q}{\rho }{E}_{\rm g}\left(x,y\right)F\left(z\right), \end{array}$$

where q is the local dissipation efficiency, \(\rho\) is the reference density of seawater (kg m−3), \({E}_{\rm g}\) (W m−2) is the rate of conversion of barotropic tidal energy into internal waves (related to the generation of internal waves), and \(F\left(z\right)\) is a vertical distribution function that assumes \(\varepsilon\) decays exponentially away from the ocean bottom as:

$$\begin{array}{c}F\left(z\right)=\frac{{\rm exp}\left(-\frac{H+z}{\zeta }\right)}{\zeta \left(1-{\rm exp}\left(-\frac{H}{\zeta }\right)\right)}, \end{array}$$

where z is the vertical coordinate (positive upward), \(H\) is the bottom depth, and \(\zeta\) is the vertical decay scale. Vertical diffusivity, \({\kappa }_{{\rm NEAR}}\), is calculated following the Osborn (1980) relationship for the mechanical energy budget of turbulence as \({\kappa }_{{\rm NEAR}}=\Gamma \varepsilon {N}^{-2}\), where \(\Gamma\) is the mixing efficiency of turbulence (assumed to be 0.2) and \(N\) is the buoyancy frequency.

We use \({E}_{\rm g}\) values of the major semidiurnal (M2 and S2) and diurnal (K1 and O1) tidal constituents estimated from the three-dimensional global numerical simulation of Niwa and Hibiya (2011). While St. Laurent et al. (2002) used \(q=0.3\), more recent studies using a numerical model have indicated that q can take larger values where tidal frequency is subinertial (lower than the inertial frequency; Tanaka et al. 2010). We assume that q takes different values for subinertial and superinertial tidal frequencies, referred to as \({q}_{{\rm sub}}\) and \({q}_{{\rm sup}}\), respectively, and use \({q}_{{\rm sub}}\) only for diurnal constituents poleward of 30°. We use \({q}_{{\rm sub}}\), \({q}_{{\rm sup}}\), and \(\zeta\) as the control variables in the optimization, and set their first-guess values as \({q}_{{\rm sub}}^{0}=1\), \({q}_{{\rm sup}}^{0}=0.3\), and \({\zeta }^{0}=500\) m.

Far-field mixing

For the FFM scheme that models the mixing induced by the cascade of internal waves to dissipation scales, we adopted Hibiya et al. (2006) instead of Oka and Niwa (2013), with an emphasis on consistency with observations. The Hibiya et al. (2006) scheme is based on the empirical relationship between diapycnal diffusivity in the thermocline (estimated from observed fine-scale vertical shear) and the numerically predicted available energy density of the semidiurnal internal tide (the major energy source for diapycnal diffusivity) at each location. The vertical diffusivity, \({\kappa }_{{\rm FAR}}\), at longitude (\(\theta\)) and latitude (\(\phi\)) is calculated as:

$${\kappa }_{{\rm FAR}}=\left({F}_{{\rm FAR}}\left(\theta ,\phi \right)+{b}_{{\rm FAR}}\right)\times \frac{\phi }{20}\, {\rm for} \,0\le |\phi |<20$$
$${\kappa }_{{\rm FAR}}={F}_{{\rm FAR}}\left(\theta ,\phi \right)+{b}_{{\rm FAR}}\, {\rm for} \, 20\le |\phi |<30$$
$${\kappa }_{{\rm FAR}}={F}_{{\rm FAR}}\left(\theta ,\phi \right)\times \frac{35-\phi }{5}+{b}_{{\rm FAR}}\, {\rm for} \, 30\le |\phi |<35$$
$${\kappa }_{{\rm FAR}}={b}_{{\rm FAR}}\, {\rm for} \,35\le \left|\phi \right|$$

with

$$\begin{aligned}{F}_{{\rm FAR}}\left(\theta ,\phi \right)=&\,{r}_{{\rm FAR}}\times {{\rm log}}_{10}\left(\frac{{E}_{{\rm internal}}\left(\theta ,\phi \right)}{0.1}\right) \,\\ &\quad {\rm for} \,{E}_{{\rm internal}}\left(\theta ,\phi \right)\ge 0.1 \; {\rm J }{{\rm m}}^{-3}\end{aligned}$$
$${F}_{{\rm FAR}}\left(\theta ,\phi \right)=0\, {\rm for} \,{E}_{{\rm internal}}\left(\theta ,\phi \right)<0.1 \; {\rm J}{{\rm m}}^{-3},$$

where \({b}_{{\rm FAR}}\) and \({r}_{{\rm FAR}}\) are unknown parameters and \({E}_{{\rm internal}}\) is the available energy density of the internal tide. The dependency on latitude is determined under the assumption that the cascade of internal waves is dominated by parametric subharmonic instability, and is consistent with observations. We used \({E}_{{\rm internal}}\) estimated by Niwa and Hibiya (2011) for the four major tidal constituents. We assumed that vertical diffusivity, \({\kappa }_{{\rm FAR}}\), is vertically constant, on the basis of the consensus from observations of the relationship between dissipation and stratification in the open ocean away from rough topography (Waterhouse et al. 2014). We used \({b}_{{\rm FAR}}\) and \({r}_{{\rm FAR}}\) (representing background diffusivity and the ratio between diffusivity and \({{{\rm log}}_{10}(E}_{{\rm internal}})\), respectively) as control variables in the optimization, and set their first-guess values as \({b}_{{\rm FAR}}^{0}=0.2 \times {10}^{-4}\) m2 s−1 and \({r}_{{\rm FAR}}^{0}=0.67 \times {10}^{-4}\) m2 s−1.

Vertical diffusivity, κ, in our model is represented as the sum of \({\kappa }_{{\rm NEAR}}\), \({\kappa }_{{\rm FAR}}\), and the diffusivity from the other schemes for the surface mixed layer, double diffusion, and convection. Figure 1 shows κ in the control (CTL) run using the first-guess values for all parameters (Table 1). The overall pattern is very similar over a wide range of depths (e.g. 1000 m and 3000 m; Fig. 1a, b), with this pattern representing \({\kappa }_{{\rm FAR}}\), which depends on latitude and \({E}_{{\rm internal}}\), but not depth. Vertical diffusivity is apparently greater at 4000 m than at 1000 m (Fig. 1c). This enhancement of κ in the deeper ocean is associated with \({\kappa }_{{\rm NEAR}}\), which has a bottom-confined structure; κ is large around rough topography like the Central Ridge in the Atlantic Ocean and the Hawaiian Ridge, because \({E}_{{\rm g}}\) is large and the seabed is relatively shallow, rising to near this depth. Large κ values in closed narrow basins, such as those off the west coast of South America, are associated with weak stratification. This could be exaggerated in a coarse-resolution model, but does not violate energy constraints, as the scheme is based on the energy dissipation rate. Similar enhancement of κ around rough topography also occurs at 3000 m depth. At 1000 m depth, κ values are relatively high around shallow sills, such as the Kuril Straits and the Aleutian Passes, as associated with \({\kappa }_{{\rm NEAR}}\). Values of κ exceeding 100 × 10−4 m2 s−1 at 1000 m depth off Antarctica and Greenland reflect unstable stratification associated with dense water formation.

Fig. 1
figure 1

Vertical diffusivity (m2 s−1) in the control (CTL) case at a 1000 m, b 3000 m, and c 4000 m

Full size image
Table 1 First-guess, perturbed, and optimized values of the control variables, with their prescribed standard deviations, \({\sigma }_{j}\)
Full size table

2.3 Optimization

A Green’s function approach (Menemenlis et al. 2005) was applied in searching for an optimal set of model parameter values, by minimizing the following cost function:

$$J={J}_{{\rm bg}}+{J}_{{\rm obs}}=\frac{1}{2}[{\varvec{a}}-{{\varvec{a}}}_{0}{]}^{{\rm T}}{\mathbf{B}}^{-1}\left[{\varvec{a}}-{{\varvec{a}}}_{0}\right]+\frac{1}{2}[H\left({\varvec{a}}\right)-{\varvec{y}}{]}^{{\rm T}}{\mathbf{R}}^{-1}\left[H\left({\varvec{a}}\right)-{\varvec{y}}\right],$$

where \({J}_{{\rm bg}}\) and \({J}_{{\rm obs}}\) are background and observation costs, respectively, \({\varvec{a}}\) and \({{\varvec{a}}}_{0}\) represent the vector of control variables and their first-guess values, respectively, \(\mathbf{B}\) and \(\mathbf{R}\) are the background and observation error covariance matrices, respectively, and H is a nonlinear operator that transforms \({\varvec{a}}\) to a model equivalent of the observation vector \({\varvec{y}}\). We assume that control variables are independent, and therefore set the background matrix as \(\mathbf{B}={\rm diag}[{\sigma }_{1}^{2},{\sigma }_{2}^{2},\dots ,{\sigma }_{m}^{2}],\) where \({\sigma }_{j}^{2}\) is the variance of the jth control variable, and \(m\) is the number of control variables. \(\mathbf{R}\) is assumed to be diagonal, whose diagonal term is determined based on the squared difference between model and WOA observation, as well as on the grid volume corresponding to each term. For observational data we used WOA98; monthly climatological data were used for the upper layer shallower than 2000 m, and annual climatological data for the remaining layer deeper than 2000 m.

The vector \({\varvec{a}}\) comprises 12 control variables that are globally constant: 5 are parameters used in the two tidally induced vertical diffusivity schemes; 5 are empirical parameters associated with diffusivity for interior mixing as selected by Toyoda et al. (2015), including double-diffusion parameterizations due to salt fingering and diffusive convection, \({f}_{{\rm sf}}\) and \({f}_{{\rm dc}}\), respectively; isopycnal diffusivity, \({A}_{{\rm isp}}\) (Redi 1982); thickness diffusivity of eddy-induced transport parameterization, \({A}_{{\rm thk}}\) (Gent and McWilliams 1990); diffusivity between the bottom boundary layer and the layer immediately above, \({K}_{{\rm BBL}}\). The two remaining variables are horizontal diffusivity, \({A}_{{\rm h}}\), and a parameter for the tapering function for isopycnal and thickness diffusion, \({S}^{{\rm dm}}\) (Danabasoglu and McWilliams 1995). The tapering function is given by:

$${\rm taper}=0.5\left[1-{\rm tanh}\left(\frac{\left|\mathbf{S}\right|-{\delta }^{{\rm dm}}}{{S}^{{\rm dm}}}\right)\right],$$

where \(\mathbf{S}\) is the isoneutral slope vector defined as:

$$\mathbf{S}=-\left(\frac{{\nabla }_{{\rm h}}\rho }{{\rho }_{{\rm z}}}\right).$$

\({S}^{{\rm dm}}\) is the half-width over which the function changes from ~ 1 to ~ 0, and \({\delta }^{{\rm dm}}\) is the slope at which the function takes a value of 0.5. We assumed that \({\delta }^{{\rm dm}}=4{S}^{{\rm dm}}\).

If a linear approximation is used for the observation operator, the cost function has the following stationary solution:

$$\widehat{{\varvec{a}}}={{\varvec{a}}}_{0}-[{{\mathbf{B}}^{-1}+\mathbf{D}\mathbf{H}}^{{\varvec{T}}}{\mathbf{R}}^{-1}\mathbf{D}\mathbf{H}{]}^{-1}{\mathbf{D}\mathbf{H}}^{{\varvec{T}}}{\mathbf{R}}^{-1}\left[\mathbf{H}\left({{\varvec{a}}}_{0}\right)-{\varvec{y}}\right],$$

whose uncertainty is represented with the analysis error covariance:

$$\mathbf{P}=[{{\mathbf{B}}^{-1}+\mathbf{D}\mathbf{H}}^{T}{\mathbf{R}}^{-1}\mathbf{D}\mathbf{H}{]}^{-1}.$$

The derivative \(\mathbf{D}\mathbf{H}\) is approximated by a Green’s function (Menemenlis et al. 2005):

$${(\mathbf{D}\mathbf{H})}_{{\varvec{i}},{\varvec{j}}}\cong \frac{{\rm H}_{{\varvec{i}}}\left({{\varvec{a}}}_{{\varvec{o}}}+{{\varvec{E}}}_{{\varvec{j}}}\right)-{{\rm H}}_{{\varvec{i}}}\left({{\varvec{a}}}_{{\varvec{o}}}\right)}{{{\varvec{e}}}_{{\varvec{j}}}}$$

where \({\rm H}_{{\varvec{i}}}\) is the ith component of the observation operator, and \({{\varvec{E}}}_{{\varvec{j}}}\) is the perturbation vector that has all zero components except the jth, \({{\varvec{e}}}_{{\varvec{j}}}\). In addition to CTL that uses the first-guess values for all of the control variables, \({{\varvec{a}}}_{{\varvec{o}}}\), we conducted 12 perturbed model integrations, each of which perturbed only one control variable from \({{\varvec{a}}}_{{\varvec{o}}}\) (Table 1) to obtain the derivative, \(\mathbf{D}\mathbf{H}\). We then computed the optimal solution \(\widehat{{\varvec{a}}}\) (Table 1) by applying the above solution formula. We also integrated the model using \(\widehat{{\varvec{a}}}\), and refer to the climatological state, \({\rm H}\left(\widehat{{\varvec{a}}}\right)\), as the optimized (OPT) case. To evaluate the error derived from the linear assumption, we also constructed a climatological state predicted using the Green’s function as:

$${H}^{{\rm GF}}=H\left({{\varvec{a}}}_{{\varvec{o}}}\right)+\mathbf{D}\mathbf{H}\left(\widehat{{\varvec{a}}}-{{\varvec{a}}}_{{\varvec{o}}}\right),$$

which we refer to as the Green’s function (GF) case.

3 Results and discussion

3.1 Optimized vertical diffusivity

As a result of the optimization, \({r}_{{\rm FAR}}\) decreased, and \({b}_{FAR}\) increased, with the \({\widehat{r}}_{{\rm FAR}}\) and \({\widehat{b}}_{{\rm FAR}}\) values being acceptable when compared with microstructure observations. Hibiya et al. (2007) estimated κ values based on microstructure observations around the North Pacific Ocean, and investigated the relation between κ and \({E}_{{\rm internal}}\) as in Hibiya et al. (2006), who used κ estimated from fine-scale parameterization. The Hibiya et al. (2007) results agree qualitatively with Hibiya et al. (2006), but suggest that the observed mixing hotspot from 20° N to 30°N, which is related to \({r}_{{\rm FAR}}\), is 2–3 times less intense. The \({\widehat{r}}_{{\rm FAR}}\) of this study is consistent with their result. Hibiya et al. (2007) also showed that in the region poleward of 30° N and equatorward of 20° N, where κ is not related to \({E}_{{\rm internal}}\), but is determined by \({b}_{{\rm FAR}}\) in the scheme of Hibiya et al. (2006), the maximum κ value is 0.3 × 10−4 m2 s−1, and most values are < 0.1 × 10−4 m2 s−1. The \({\widehat{b}}_{{\rm FAR}}\) value in this study is below this observed upper limit. Waterhouse et al. (2014) compiled 5200 microstructure observations from all over the world, and calculated mean profiles of κ for each of three different topographic regimes. \({\widehat{b}}_{{\rm FAR}}\) lies within the range of their mean profiles, when we consider the value well above the topography where the influence from the bottom intensified mixing is small. Considering the nature of turbulence and its variability, κ estimated on the basis of limited microstructure observations has a large uncertainty. Therefore, it remains to be tested whether these \({\widehat{r}}_{{\rm FAR}}\) and \({\widehat{b}}_{{\rm FAR}}\) values are realistic by using more observational data and/or more sophisticated statistical approaches based on a generative model for the observed profiles (Sugiura et al. 2019).

The \({\widehat{q}}_{{\rm sup}}\) value is ~ 50% higher than \({q}_{{\rm sup}}^{0}\), but is consistent with the global mean q for M2 estimated by Vic et al. (2019), who suggest that small-scale internal tides, hitherto overlooked, account for > 50% of global internal tide generation, breaking, and mixing. A \({\widehat{q}}_{{\rm sub}}\) value larger than unity means that the energy constraints are violated, if NFM is caused only by the direct breaking of internal tides. However, NFM is also related to the wave breaking through the interaction with the background internal waves, which may supply additional energy. In addition, this excess is small enough to be neglected relative to the differing estimates of \({E}_{{\rm g}}\) among different methods, tidal models, and satellite data (e.g. Egbert and Ray 2003).

Note that \(\widehat{\zeta }\) is significantly larger than \({\zeta }^{0}\). This may be related to uncertainties in \(\zeta\) and/or F(z). Schmittner and Egbert (2014) suggested that the underestimation of mixing at shallow depths, as a result of unrealistic representation of narrow subgrid-scale topographic features, may be a reason for MOC being weak in coarse-resolution ocean models. This is consistent with \(\widehat{\zeta }\) being larger than \({\zeta }^{0}.\) Although we assume \(\zeta\) is constant globally, it depends on the wave number of the bottom topography and the tidal flow amplitude (Hibiya et al. 2017), and this could be another reason for the underestimation of mixing at shallow depths. Although we used an exponentially decaying profile as in St. Laurent et al. (2002), Polzin (2004) proposed another profile that decays algebraically on the basis of an one-dimensional radiation balance equation. Melet et al. (2013) showed that the use of this latter profile results in stronger mixing at shallow depths and significantly modifies stratification and the MOC.

An error range was added to each optimal value in Table 1. Although quantitative and objective assessments are difficult to make, the posterior error decreased compared with the prior error (\({\sigma }_{j}\)). Among five parameters for tidally induced mixing schemes, the decrease is relatively small for \({q}_{{\rm sub}}\) and \(\zeta\). This may suggest that the ocean state is less sensitive to—and/or the assumption of a global constant is less reasonable for—those parameters.

Although κ in OPT seems very similar to κ in CTL (Fig. 1) with the chosen color scale (not shown), the deviation of κ from CTL to OPT shows a characteristic pattern (Fig. 2). The modification of \({\kappa }_{{\rm FAR}}\) is clearly seen in the deviation of κ at 1000 m depth (Fig. 2a). The decrease of κ in low-latitude regions and its moderate increase in most remaining regions are consistent with the decrease of \({r}_{{\rm FAR}}\) and increase of \({b}_{{\rm FAR}}\), respectively. This pattern of deviation of κ at 1000 m depth is completely different from that at 3000 m depth (Fig. 2b), which shows increases around rough topography. This is consistent with the increase of both \({q}_{{\rm sub}}\) and \({q}_{{\rm sup}}\). Note that the increase of \(\zeta\) leads not only to the increase of ε at shallow depths, but also to the decrease of ε near the bottom. The decrease of κ near the bottom topography at 4000 m depth must be at least partly related to this effect (Fig. 2c).

Fig. 2
figure 2

Deviation of vertical diffusivity (m2 s−1) in the optimized (OPT) case from the control (CTL) case at a 1000 m, b 3000 m, and c 4000 m

Full size image

3.2 Potential temperature and salinity

The observation cost decreased by a total of 11% (Table 2). By component, the cost for T decreased by 11%, and that for S decreased by 10%, meaning that the climatological state is improved in both T and S. A significant decrease in cost occurs in the Pacific Ocean, especially for T, while a relatively large cost increase occurs for T in the North Atlantic Ocean and the eastern Indian Ocean (Fig. 3). This spatial distribution of positive and negative signs is well predicted by the GF case. The cost increase in GF is weaker than in OPT, and is insignificant in comparison with the cost decrease seen broadly across the Pacific Ocean. This suggests that the number of control variables may be too small to reduce the cost globally, but the linear optimization performed well, especially in the Pacific Ocean.

Table 2 Observation cost values for temperature (T), salinity (S), and total, in the control (CTL), optimized (OPT), and Green’s function (GF) cases
Full size table
Fig. 3
figure 3

Vertically integrated cost difference over the 0–5500 m depth range for temperature (T), salinity (S), and T + S, between the optimized (OPT) case and the control (CTL) case (top), the Green’s function (GF) case and CTL (middle), and OPT and GF (bottom)

Full size image

Regions with relatively large cost increases in OPT correspond to regions where the cost deviation from the GF to OPT is large and thus the linear assumption causes large errors. An iterative approach as well as additional control variables may be effective in order to reduce costs in those regions.

We note that the change in cost is generally small in regions where T and S are restored to WOA98 to simulate the effects of the processes that cannot be reproduced by the model. Since restoration always reduces differences between models, the sensitivity of T and S to control variables must be inherently weak in those restoring regions in our model. This can explain the small cost change in those restoring regions, and may explain why the predicted cost change in the GF case is weak but positive in and around some restoring regions (e.g. that for T in the North Atlantic Ocean and in the Southern Ocean). Considering that deep water formation is related to the restoration in our model and is thus less sensitive to the control variables, it is possible that the difference between OPT and CTL for deep water is caused mainly by the change of circulation and/or the modification of water properties by mixing and/or geothermal heating along the path of advection, rather than that of water mass formation.

From here on, we focus on the Pacific Ocean, where the cost is decreased considerably and the linear assumption likely worked well. When compared with the observation-based climatology of WOA98, deep water (4000–5500 m) has lower temperature and higher salinity in the North Pacific in CTL (Fig. 4). This cold bias is apparently reduced in OPT. This temperature distribution is improved even when compared with the optimized state of the former version of this model (Toyoda et al. 2015): water colder than 0.7 °C from the Southern Ocean intrudes farther north, and the meridional gradient of temperature along the western boundary strengthens; intrusion of water colder than 1 °C to the eastern North Pacific weakens; and the warm bias south of Australia is reduced. The high salinity bias in the North Pacific is reduced in OPT, although the meridional salinity gradient is still weaker than in WOA98. When we compare with the optimized state of Toyoda et al. (2015), this gradient is weaker and thus less realistic, but there are some improvements in OPT: in Toyoda et al. (2015) the location of the 34.7 psu contour (Fig. 4) is farther south, and the salinity is significantly lower in the Philippine Sea in their optimized state than in WOA98, but these are more realistic in OPT.

Fig. 4
figure 4

Climatological temperature and salinity in the deep-ocean layer (4000–5500 m) in observations (OBS), the control (CTL) case, and the optimized (OPT) case

Full size image

To evaluate the contributions of the two tidally induced vertical diffusivity schemes to this improvement of the deep-ocean state, we conducted three additional experiments: (1) a near-field mixing (NFM) run using \(\widehat{{\varvec{a}}}\) only for \({q}_{{\rm sub}}\), \({q}_{{\rm sup}}\), and \(\zeta ;\) a far-field mixing (FFM) run using \(\widehat{{\varvec{a}}}\) only for \({b}_{{\rm FAR}}\) and \({r}_{{\rm FAR}}\); and an “other parameters” (OTH) run using \(\widehat{{\varvec{a}}}\) only for the remaining seven parameters. We refer to a deviation of experimental results from CTL by the term \(\Delta\) with a subscript (e.g. \(\Delta {T}_{{\rm OPT}}\)).

The \(\Delta T\) and \(\Delta S\) for the deep ocean, 4000–5500 m, are shown in Figs. 5 and 6, respectively. \(\Delta {T}_{{\rm OPT}}\) shows that the deep water warms throughout the Pacific Ocean as a result of the optimization. This warming increases northward, intensifying the positive meridional gradient of T. \(\Delta {T}_{{\rm OPT}}\) is well reproduced by the sum of the three components, \(\Delta {T}_{{\rm NFM}}+\Delta {T}_{{\rm FFM}}+\Delta {T}_{{\rm OTH}}\), supporting the suggestion from Fig. 3 that the linear assumption works well. \(\Delta {T}_{{\rm NFM}}\) and \(\Delta {T}_{{\rm OTH}}\) contribute to the improvement of T, the warming and the intensification of the meridional gradient, at about the same level. \(\Delta {T}_{{\rm FFM}}\) has a slight cooling effect.

Fig. 5
figure 5

Differences of temperature in the deep layer (4000–5500 m) with respect to the control (CTL) run in the optimized (OPT), near-field mixing (NFM), far-field mixing (FFM), and other parameters (OTH) runs. \(\Delta {T}_{{\rm GT}}\) is the difference in CTL compared with the no geothermal run, showing the effect of geothermal heating (GT)

Full size image
Fig. 6
figure 6

Same as in Fig. 5, but for salinity in the deep layer (4000–5500 m)

Full size image

\(\Delta {S}_{{\rm OPT}}\) indicates that the deep water is freshened throughout the Pacific Ocean as a result of the optimization, increasing northward. This contrast intensifies and thus improves the negative meridional gradient of S. \(\Delta {S}_{{\rm OPT}}\) is well reproduced by the sum of the three components, \(\Delta {S}_{{\rm NFM}}+\Delta {S}_{{\rm FFM}}+\Delta {S}_{{\rm OTH}}\), meaning that the linear assumption works well also for S. All of the three components are negative throughout most of the Pacific Ocean, while \(\Delta {S}_{{\rm FFM}}\) is very weak. Both \(\Delta {S}_{{\rm NFM}}\) and \(\Delta {S}_{{\rm OTH}}\) intensify the meridional gradient of S to about the same extent.

These results clearly show that modification of near-field mixing, \(\Delta {\kappa }_{{\rm NFM}}\), effectively improves the climatological state of the Pacific deep ocean, at least in our OGCM. Note that this effect is as large as the effect of modification of other parameters including isopycnal and thickness diffusivities. This contrasts with the results of the data synthesis experiment using 4D-VAR by Liu et al. (2012), who showed that the effect of modification of vertical mixing is one order of magnitude weaker than those of isopycnal and thickness diffusivity. Since the magnitude of the effect of modifying each parameter depends on their first guess values, it is difficult to compare the relative importance. However, this supports the importance of NFM.

We conducted another optimization experiment in which the two parameters for FFM, \({b}_{{\rm FAR}}\) and \({r}_{{\rm FAR}}\), were excluded from the control variables. The experiment yielded an 11% decrease in observation cost, as compared to 13% in GF case (Table 2). It is not surprising that the magnitude of the decrease in cost is reduced, given the reduction in the number of control variables. However, in the additional experiment the adjustment of the three parameters for NFM is weaker, and the reproducibility of the predicted deep-water properties is lower (not shown) in comparison with GF case. This suggests that the adjustment of parameters in FFM indirectly contributes to the improved reproducibility of the deep-ocean water properties.

Intensification of the basin-scale positive (negative) meridional gradient of T (S) means that the change of deep water in the North Pacific Ocean cannot be explained by the change of inflowing source water. Both \(\Delta {T}_{{\rm OPT}}\) and \(\Delta {S}_{{\rm OPT}}\) intensify the meridional density gradient. This suggests that the pressure distribution and thus the deep ocean circulation were also modified as a result of the optimization. We now discuss the impact on the Pacific meridional overturning circulation (PMOC), which is not directly constrained in the optimization.

3.3 Pacific meridional overturning circulation

Using the annual mean meridional velocity, we calculated the meridional overturning stream function, \(\psi\), in the Pacific Ocean for each experiment (Fig. 7). The lower cell, representing the inflow of cold, saline Antarctic bottom water (AABW) from the Southern Ocean, is ~ 8 Sv at 10° S in OPT. This value is slightly lower than that for the northward transport (10.6 ± 1.7 Sv) estimated using observations from hydrographic surveys and current meter moorings (Roemmich et al. 1996). However, this difference is acceptable considering the large temporal variability in values indicated by mooring observations in the Samoan Passage: Rudnick (1997) reported a mean value of 6.0 Sv (standard deviation of 1.5 Sv, a minimum of 1.1 Sv, and a maximum of 10.7 Sv) for 1992–1994, while Voet et al. (2016) reported a mean value of 5.4 Sv for 2012–2013.

Fig. 7
figure 7

Same as in Fig. 5, but for the meridional overturning stream function, \(\psi\), in the Pacific Ocean. Color shading indicates the difference in \(\psi\) from that in CTL and contour lines indicate \(\psi\) in the optimized (OPT) case

Full size image

\(\Delta {\psi }_{{\rm OPT}}\) indicates that as a result of the optimization, the lower cell connecting the Southern Ocean and North Pacific is intensified to about 30° N, while both the upper and lower cells are weakened north of 30° N. This intensification is ~ 0.5 Sv at 20° N. This suggests that \(\Delta {\psi }_{{\rm OPT}}\) reduces the discrepancy in the mean transport between CTL and the estimate based on observations, although further investigation using more observational data is required to determine if this is an improvement. As was the case for deep T and S, this \(\Delta {\psi }_{{\rm OPT}}\) is well reproduced by the sum of the three components, \(\Delta {\psi }_{{\rm NFM}}+\Delta {\psi }_{{\rm FFM}}+\Delta {\psi }_{{\rm OTH}}\). The intensification of the lower cell is explained only by \(\Delta {\psi }_{{\rm NFM}}\), while the weakening of the lower cell is explained by \(\Delta {\psi }_{{\rm NFM}}\) and \(\Delta {\psi }_{{\rm OTH}}.\) This suggests that NFM plays an important role in controlling the inflow of AABW.

The intensification of the lower cell in \(\Delta {\psi }_{{\rm NFM}}\) is reasonable, since an increase in both \({q}_{{\rm sub}}\) and \({q}_{{\rm sup}}\) strengthens upwelling, which must be strong near the bottom reflecting the bottom intensified distribution of \(\varepsilon .\) The weakening of the lower cell in both \(\Delta {\psi }_{{\rm NFM}}\) and \(\Delta {\psi }_{{\rm OTH}}\) consists of a narrow strong downwelling near 45° N, and a broad weak upwelling. This downwelling is associated largely with a downwelling at the entrance of the Sea of Okhotsk, where T and S are restored to WOA98, and strong upwelling occurs in CTL. Considering that T and S are improved in both NFM and OTH runs, it is suggested that the weakening of the lower cell is associated with the weakening of the artificial upwelling induced by the restoration.

Although beyond the scope of this study, the contribution of \(\Delta {\psi }_{{\rm FFM}}\) is significant in the upper ocean, indicating the importance of FFM in reproducing the circulation and thus water properties in the upper ocean. This \(\Delta {\psi }_{{\rm FFM}}\) value can be explained by \(\Delta {\kappa }_{{\rm FFM}}\): a decrease in \({r}_{{\rm FAR}}\) weakens the upwelling in low-latitude regions, and an increase in \({b}_{{\rm FAR}}\) strengthens upwelling in high-latitude regions.

3.4 Effects of geothermal heating

While this study focuses on the effects of tidally induced mixing, the inclusion of the geothermal heat flux provides an important improvement for our OGCM. To demonstrate its impact, we conducted another experiment using \({{\varvec{a}}}_{{\varvec{o}}}\) as CTL, but without the geothermal heat flux, referred to here as a “no Geothermal” run (noGT). We investigated the impact on deep-water T, S, and the PMOC by calculating the deviation of CTL from noGT, as represented by \(\Delta\) with subscript “GT” (\(\Delta {T}_{{\rm GT}}={T}_{{\rm CTL}}-{T}_{{\rm noGT}}\)).

The \(\Delta {T}_{{\rm GT}}\) value is positive, meaning that geothermal heating warms the deep ocean (Fig. 5). The \(\Delta {T}_{{\rm GT}}\) value is comparable with \(\Delta {T}_{{\rm OPT}}\), and larger than \(\Delta {T}_{{\rm NFM}}\). As with \(\Delta {T}_{{\rm OPT}}\), \(\Delta {T}_{{\rm GT}}\) increases northward and intensifies the meridional gradient of T. This gradient of \(\Delta {T}_{{\rm GT}}\) is, however, steeper than that of \(\Delta {T}_{{\rm OPT}}\). In the zonal pattern of the North Pacific Ocean, \(\Delta {T}_{{\rm GT}}\) has an effect opposite to that of \(\Delta {T}_{{\rm OPT}}\): \(\Delta {T}_{{\rm GT}}\) increases eastward, while \(\Delta {T}_{{\rm OPT}}\) increases westward. This eastward increase of \(\Delta {T}_{{\rm GT}}\) may be related at least partly to the distribution of the geothermal heat flux, which has strong zonal contrast with high values along the East Pacific Rise and the west coast of North America. Note that WOA98 indicates that deep water is warmer in the eastern region in the North Pacific (Fig. 4). This zonal contrast is not reproduced in OPT even with the geothermal heating, but is improved in comparison with Toyoda et al. (2015). This suggests that geothermal heating plays an important role in determining the spatial distribution of deep-ocean T.

The \(\Delta {S}_{{\rm GT}}\) value is positive everywhere except in the Southern Ocean (Fig. 6). The positive deviation is greatest around 30° N, decreasing the negative meridional gradient of S south of 30° N. This means that geothermal heating counteracts the improvement by \(\Delta {S}_{{\rm OPT}}\). This could explain why the deep-ocean S is reproduced less well in OPT when compared with Toyoda et al. (2015), although it is improved in comparison with CTL. The absolute value of \(\Delta {S}_{{\rm GT}}\) is comparable to \(\Delta {S}_{{\rm NFM}}\) and \(\Delta {S}_{{\rm OTH}}\), suggesting that geothermal effects cannot be neglected when estimating ocean states, even for S and the global budget of freshwater. The zonal structure has a meridional maximum that cannot be explained by the change of source water alone, suggesting that \(\Delta {S}_{{\rm GT}}\) is influenced by changes in circulation.

\(\Delta {\psi }_{{\rm GT}}\) indicates that the lower cell is intensified to about 30° N, but weakened north of 30° N. This overall pattern is similar to \(\Delta {\psi }_{{\rm NFM}}\) (Fig. 7). However, the intensification (weakening) of the lower cell within the North Pacific south (north) of 30° N is weaker (stronger) than \(\Delta {\psi }_{{\rm NFM}}\), meaning that \(\Delta {\psi }_{{\rm GT}}-\Delta {\psi }_{{\rm NFM}}\) shows a weakening of the lower cell throughout the North Pacific. This suggests that geothermal heating is less effective in increasing northward transport of bottom water to the northern North Pacific than \(\Delta {\kappa }_{{\rm NFM}}\).

It appears that the meridional maximum of \(\Delta {S}_{{\rm GT}}\) corresponds to the region where \(\Delta {\psi }_{{\rm GT}}\) indicates an upward circulation deviation. This may support the idea that \(\Delta {S}_{{\rm GT}}\) is influenced by the change of circulation as we have suggested.

The impact of adding the geothermal heat flux is thus comparable with that of adjusting the mixing parameters. This suggests that the geothermal heating cannot be neglected in order to provide reliable simulations that yield a realistic diffusivity representativeness, while the influence of uncertainty in the geothermal heat flux may be relatively small. It is worthy of note that the relation between \(\Delta T \: {\rm and} \: \Delta S\) differs for the impacts of geothermal heating and optimization. This implies that adjusting geothermal heat flux may be effective, at least qualitatively, in improving model simulation of the deep ocean.

4 Summary and future work

Using a Green’s function approach, we have conducted an optimization experiment for an OGCM using tidally induced vertical diffusivity schemes for NFM and FFM for interior mixing below a surface mixed layer rather than the semi-empirical schemes used in the earlier version of our model (Toyoda et al. 2015). The optimization improves both T and S in the deep Pacific Ocean, and the deep ocean state is represented in the optimized model reasonably well. This optimized deep ocean state is at the same level as the optimized state in the earlier version of our model: T is apparently improved, while an improvement in S is less clear, as its reproducibility seems to be reduced when considering the meridional gradient, although some improvements are apparent (e.g., the location of the 34.7 psu contour). This represents a meaningful advance, considering also that our new model provides a realistic diffusivity representativeness. As a result of the optimization, the lower cell of PMOC is intensified. We have shown that the adjustment of parameters in the NFM scheme works effectively to improve the deep-ocean water properties and to intensify PMOC. An additional optimization experiment, in which the two parameters for FFM were excluded from the control variables, suggests that the adjustment of parameters in FFM indirectly contributes to the improved reproducibility of the deep-ocean water properties. These results indicate that the combination of two tidally induced diffusivity schemes, NFM and FFM, works well, with NFM being especially important for the deep-ocean state of the Pacific Ocean.

We have also shown that the impact of adding the geothermal heat flux is comparable with that of adjusting the mixing parameters, thereby demonstrating the importance of geothermal heating in performing reliable simulations with a realistic diffusivity representativeness.

It is likely that further improvement of the mixing scheme will provide more reliable simulations of deep-ocean states that are more realistic and consistent with current understanding of mixing. The vertical structure of mixing is one of the most important factors in improving tidally induced mixing schemes for both NFM and FFM. We used the FFM scheme related to the cascade to smaller scales through wave–wave interactions in the open ocean, and it may be important to account for FFM in other mechanisms, such as the evolution of internal waves on continental slopes and shelves (e.g. Nash et al. 2004; Klymak et al. 2011). It is also important to assess the validity of the optimized parameters based on observations. While we focused only on tidally induced mixing, mixing that is related to processes involving winds (e.g. Alford et al. 2016) and geostrophic currents (e.g. Bell 1975) should also be considered.

In the context of ocean state estimation, the three-dimensional distributions of vertical and other diffusivities have recently been adjusted through a more sophisticated data assimilation technique known as the four-dimensional variational adjoint method, with the adjusted diffusivities being effective in reducing the difference between OGCMs and temperature and salinity observations (e.g. Liu et al. 2012; Forget et al. 2015). Nevertheless, a variational method usually requires consistent first-guess control variables, such as diffusivity parameters. Our approach, of adjusting parameters in mixing schemes on the basis of the recent knowledge, provides a basis for ocean state estimation that is energetically more reliable. Such estimations may contribute to refinements of the assessment of the present climate state and understanding of the role of the ocean in the Earth’s climate system.