Development and validation of a global 1∕32° surface-wave–tide–circulation coupled ocean model: FIO-COM32

Model resolution and the included physical processes are two of the most important factors that determine the realism or performance of ocean model simulations. In this study, a new global surface-wave–tide–circulation coupled ocean model FIO-COM32 with a resolution of $\mathrm{1}/\mathrm{32}$^∘ × $\mathrm{1}/\mathrm{32}$^∘ is developed and validated. Promotion of the horizontal resolution from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘ leads to significant improvements in the simulations of surface eddy kinetic energy (EKE), the main paths of the Kuroshio and Gulf Stream, and the global tides. We propose the integrated circulation route error (ICRE) as a quantitative criteria to evaluate the simulated main paths of the Kuroshio and Gulf Stream. The non-breaking surface-wave-induced mixing (B_V) is proven to still be an important contributor that improves the agreement of the simulated summer mixed-layer depth (MLD) against the Argo observations even with a very high horizontal resolution of $\mathrm{1}/\mathrm{32}$^∘. The mean error in the simulated mid-latitude summer MLD is reduced from −4.8 m in the numerical experiment without B_V to −0.6 m in the experiment with B_V. By including the global tide, the global distributions of internal tide can be explicitly simulated in this new model and are comparable to the satellite observations. Based on Jason-3 along-track sea surface height (SSH), wavenumber spectral slopes of mesoscale ranges and wavenumber frequency analysis show that the unbalanced motions, mainly internal tides and inertia-gravity waves, induced SSH undulation and are a key factor for the substantially improved agreement between model and satellite observations in the low latitudes and low-EKE regions. For the ocean model community, surface waves, tidal currents and ocean general circulations have been separated into different streams for more than half a century. This paper demonstrates that it is time to merge these three streams for a new generation of ocean model development.

1 Introduction

As available computing power has been increasing rapidly, by approximately 1 order of magnitude every 5 years, the state-of-the-art computing ability of modern global ocean numerical models is becoming enormously high, leading to recent achievements in high-resolution global ocean models. The term of “high resolution” as regards present global ocean models may refer to those with horizontal resolutions ranging from 1 to 5 km, which are well beyond the mesoscale resolving threshold in most open oceans (Hallberg, 2013). Further improved resolution has a significant impact on the simulated eddy activities (Thoppil et al., 2011; Sasaki and Klein, 2012; Biri et al., 2016; Ajayi et al., 2020), the vertical mass and buoyancy fluxes (Capet et al., 2016; Su et al., 2018; Dong et al., 2020a), and even the representation of large-scale circulations (Lévy et al., 2010; Chassignet and Xu, 2017). These high-resolution models are not only helpful in improving our scientific understanding of mesoscale and sub-mesoscale eddies and mixed-layer dynamics but are also very useful when evaluating the satellite products of both present (Amores et al., 2018) and future generations, such as SWOT projects (Uchida et al., 2022). Table 1 summarises recent developments in high-resolution global ocean models. Not only is the model horizontal resolution significantly increased, which now ranges from $\mathrm{1}/\mathrm{20}$ to $\mathrm{1}/\mathrm{48}$^∘, but the ocean model physical processes in ocean models are also making notable progresses. For example, HYCOM25 and LLC4320 are global tide–circulation coupled models, and FIO-COM32 is a global surface-wave–tide–circulation coupled ocean model.

Table 1 Recent developments in global high-resolution ocean models.

Improving the representation of physical processes in ocean models has been the most fundamental aspect for ocean general circulation model (OGCM) development. Since the establishment of the first OGCM (Bryan and Cox, 1967), surface wave models, tide models, ocean internal wave models and OGCMs have been separated into different streams (Mellor and Blumberg, 2004; Qiao et al., 2004). The most uncertain term in all OGCMs is ocean turbulence. As a result, the vertical structures of ocean temperature and salinity cannot be accurately simulated or predicted. For example, the simulated mixed-layer depth (MLD) in the upper ocean is too shallow, and the sea surface temperature (SST) overheats during summer in nearly all OGCMs. Teixeira and Belcher (2002) investigated the interaction of background turbulence with a monochromatic irrotational surface wave and found that the Stokes drift associated with the surface wave tilts vertical vorticity into the horizontal direction, subsequently stretching it into elongated streamwise vortices, and generates vertical mixing. Qiao et al. (2004, 2016) proposed an upper-ocean mixing scheme of B_V induced by surface waves, which plays a dominant role in turbulent mixing of the upper ocean, and analytically expressed B_V as a function of the wavenumber spectrum that can be directly calculated from a wave model (Qiao et al., 2004, 2017). Ocean surface wave–circulation coupling has become one of the most important directions of ocean model development. B_V has been widely adopted in a series of different OGCMs and climate models, and all models showed dramatic improvements (Qiao et al., 2004; Shu et al., 2011; Song et al., 2012; Fan and Griffies, 2014; Wang et al., 2019). B_V is particularly effective in remedying the shallow biases of the simulated MLD and overestimated SST in summer. B_V brought the following two surprising results into OGCMs. First, after closing the traditional vertical turbulence schemes, only B_V can properly simulates the global ocean, which indicates that B_V plays dominant role in vertical mixing in the upper ocean (Qiao and Huang, 2012). Second, B_V can cause the MLD to become shallower in winter in climate models, which is also an improvement for climate models (Chen et al., 2018). Although the effect of B_V has been comprehensively verified in many coarse-resolution OGCMs, most of these numerical experiments were based on climatological simulations and validated against climatological observations. In order to investigate the impacts of horizontal resolution on the model MLD biases, taking the significantly different behaviours of coarse- and high-resolution models into consideration, it is still necessary to validate the effect of B_V in the high-resolution OGCM. For the first time, this paper investigate the B_V effects in a high-resolution OGCM of $\mathrm{1}/\mathrm{32}$^∘ validated against in situ observations.

Ocean tides have been recognised as a fundamental aspect regulating the hydrodynamic environments in shallow regions (e.g. Simpson and Hunter, 1974; Garrett et al., 1981; Holt and Umlauf, 2008; Lin et al., 2020); thus, ocean tides are often included in regional ocean models. The tidal mixing controls the formation of thermal fronts in coastal regions and generates upwelling (Lü et al., 2006, 2008, 2010). The regional surface-wave–tide–circulation coupled model of the sea around China has shown excellent performance in revealing the three-dimensional circulation structure of Yellow Sea cold water mass (Xia et al., 2006) and has been applied in operational ocean forecasting systems (Wang et al., 2016). In recent years, high-resolution global OGCMs forced by atmospheric forcings begins to include ocean tide explicitly (Arbic et al., 2010, 2012, 2018; Rocha et al., 2016). One of the benefits of including ocean tide in eddy-resolving global OGCMs is that the global distribution of internal tide and mesoscale eddies can be resolved concurrently (e.g. Arbic et al., 2010, 2012, 2018; Buijsman et al., 2015; Shriver et al., 2012; Ansong et al., 2018; Timko et al., 2019). In addition, inclusion of ocean tide in a global $\mathrm{1}/\mathrm{48}$^∘ ocean model (MITgcm llc4320) leads to more realistic representation of the unbalanced and ageostrophic motions (Rocha et al., 2016).

To examine whether high-resolution OGCMs can faithfully reproduce the ocean environment, the comparison of wavenumber and frequency spectra against observations is widely adopted (e.g. Sasaki and Klein, 2012; Richman et al., 2012; Chassignet and Xu, 2017, hereafter CX17; Savage et al., 2017a, b; Biri et al., 2016). Therefore, the wavenumber spectral slope in the 70–250 km mesoscale range becomes an important criteria for OGCMs' validation. Sasaki and Klein (2012, hereafter SK12) noted that as the horizontal resolution improved to $\mathrm{1}/\mathrm{30}$^∘, the wavenumber spectral slope in the 70–250 km mesoscale range of the model shows quite inconsistent patterns with those of satellite observations. The satellite observations show strong latitudinal variability (Xu and Fu, 2012). In high-latitude and high-eddy-kinetic-energy (EKE) regions, the slopes range between −5 and −4, which is generally consistent with the theory of mesoscale turbulence. However, the slopes are much flatter in tropical and low-EKE regions, which is substantially distinct from the theory of mesoscale turbulence. In high-resolution models such as SK12 and CX17, the slopes range between −5 and −4 in both low- and high-latitude regions. More importantly, Chassignet and Xu (2021, hereafter CX21) showed that as the tide is included in their $\mathrm{1}/\mathrm{50}$^∘ Atlantic regional ocean model, the resolved internal tide in sea surface height (SSH) signals is the main explanation for the inconsistency between model and satellite observations. Although ocean tides are included in regional OGCMs, there have still been few attempts to include tides in global OGCMs.

This paper aims to answer the following three key questions through establishing a global $\mathrm{1}/\mathrm{32}$^∘ surface-wave–tide–circulation coupled ocean model (FIO-COM32). What are the main effects of increasing horizontal resolution from $\mathrm{1}/\mathrm{10}$^∘ to $\mathrm{1}/\mathrm{32}$^∘ in a global ocean model? What are the main effects of B_V in high-resolution OGCMs? Inspired by the results of an Atlantic regional high-resolution ocean model of CX21, what is the tidal effect in a global $\mathrm{1}/\mathrm{32}$^∘ high-resolution OGCM of simulating the wavenumber spectral slopes?

This paper is organised as follows. Section 2 describes the model configurations and design of the numerical experiments. In Sect. 3.1, we present the results of basic aspects of the new global $\mathrm{1}/\mathrm{32}$^∘ surface-wave–tide–circulation coupled ocean model and illustrate the effects of increased horizontal resolution by comparing it with a previous global $\mathrm{1}/\mathrm{10}$^∘ model. Section 3.2 shows the effects of ocean surface wave and tide coupling in the global $\mathrm{1}/\mathrm{32}$^∘ ocean model. Section 4 provides the summary and discussion.

2 Model description and numerical experimental design

2.1 Model description

FIO-COM consists of a global ocean circulation model, a sea ice model and a global ocean wave model. The ocean circulation component is based on the Modular Ocean Model 5 (Griffies, 2012), the sea ice component is based on the Sea Ice Simulator (Winton, 2000) and the surface wave component is based on the MASNUM (laboratory of MArine Sciences and NUmerical Modelling of the Ministry of Natural Resources of the People's Republic of China) ocean wave model (Yang et al., 2005; Qiao et al., 2017). Currently, the data exchanged between the ocean wave and the ocean circulation model are limited to B_V values. The effects of ocean surface currents on surface waves are closed pro tempore to save computer resources and are tested near the end of this paper. B_V is calculated in the ocean wave model following Qiao et al. (2004), as in Eq. (1) below, where S(k) represents the wavenumber spectrum and k is the wavenumber. As B_V is incorporated into the vertical mixing schemes of the ocean circulation model, α is a tuneable parameter for practical purposes, and we set it as 0.3 in this paper (Wang et al., 2010).

$$\begin{array}{}\text{(1)}& B_V=\mathit{\alpha }{\int }_{k}S\left(k\right)e^{-\mathrm{2}\left|k\right|z}\mathrm{d}k\cdot {\displaystyle \frac{\partial }{\partial z}}{\left[{\int }_k{\mathit{\omega }}^{\mathrm{2}}S\left(k\right)e^{-\mathrm{2}\left|k\right|z}\mathrm{d}k\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

The variable exchanged between the ocean wave and the sea ice models is sea ice concentration (SIC, from the sea ice model to the ocean wave model). The SIC is used to calculate time-varying masks of the ocean wave model, the model domain with SIC exceeding 15 % is masked out as land during the ocean wave model integration.

According to different research purposes, there are two different coupling strategies for the wave model component. Firstly, in the lower-resolution numerical experiments, where the computational costs are not expensive, the wave model component is coupled online with the ocean circulation and sea ice models. The real-time online data exchanges are achieved based on the subroutine version of the MASNUM wave model. In this way, the wave model component could be coupled with ocean circulation and ice components through direct calling of the MASNUM wave model as subroutines in the ocean circulation model. The ocean wave, ocean circulation and sea ice components share the same model grids. Secondly, in the numerical experiments with high computational costs and research focusing on ocean dynamics, the wave component can be turned off to save computational resources (the computational cost of the coupled ocean–ice model is approximately half that of the surface-wave–ocean–ice model). In this configuration, the surface-wave-induced mixing coefficients previously saved in data files are then read into the OGCM; the drawback is that the wave field cannot interact with the ocean circulation fields as they change.

The horizontal grid of FIO-COM32 is a tripolar grid (Murray, 1996) with a horizontal resolution of $\mathrm{1}/\mathrm{32}$^∘ × $\mathrm{1}/\mathrm{32}$^∘. The model covers the entire global ocean, with a latitudinal coverage of 82^∘ S to 90^∘ N and a longitudinal coverage of 280^∘ W to 80^∘ E. Vertically, the z^* coordinate (Adcroft and Campin, 2004) is adopted, and the thickness is 2 m at the surface and gradually increases to 367 m at the bottom. The vertical grid has 54 or 57 levels depending on whether ocean tide is introduced. In the numerical experiments without ocean tide, the maximum depth is 5500 m. While the ocean tide is explicitly included, the maximum depth extends to 7000 m by including additional three bottom levels in OGCM. The maximum total grid size is 11 520 × 5504 × 57. Model topography is from GEBCO (IHO-IOC, 2018) and is smoothed by applying a radial filter (Arbic et al., 2004). Model bottom cells are set as partial cells (Adcroft et al., 1997). The horizontal mixing scheme is a bi-harmonic operator with diffusive velocities of 1.96 cm s⁻¹ for momentum and 0.65 cm s⁻¹ for tracers. Both the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models use the same model settings. The viscosity is calculated from diffusive velocity times the cube of the grid spacing, which means that the viscosity is more than 30 times smaller in the $\mathrm{1}/\mathrm{32}$^∘ model than in the $\mathrm{1}/\mathrm{10}$^∘ model (FIO-COM10). The vertical mixing scheme uses the K-profile parameterization (KPP) scheme (Large et al., 1994) and B_V values. The background vertical viscosity and diffusivity are set as 1.0 × 10⁻⁴ and 3.0 × 10⁻⁵ m² s⁻¹, respectively. B_V is calculated from a global $\mathrm{1}/\mathrm{4}$^∘ online coupled FIO-COM model and incorporated into both vertical diffusivity and viscosity. Neither sea surface temperature nor salinity is restored to observational data. The bottom drag is quadratic with a coefficient of 2.5 × 10⁻³.

Ocean tide is explicitly included through introducing eight main tidal generating potentials, i.e. M₂, S₂, N₂, K₂, K₁, O₁, P₁ and Q₁, in the momentum equations following Schiller and Fiedler (2007). The treatment of self-attraction and loading (SAL) is a simple scalar approximation following Arbic et al. (2010), and the scalar alpha is set to 0.93. The SAL is applied only on the tidal elevations, thus a 25 h running average time filter is adopted to get the slowly varying sea surface height related to large-scale ocean circulation. A topographic drag scheme (Jayne and St. Laurent, 2001) is introduced in the barotropic momentum equation, and this scheme is enabled in the barotropic numerical experiments but closed in the baroclinic experiments because a large portion of the barotropic-to-baroclinic tidal energy conversion can be explicitly resolved in the baroclinic experiments.

The initial conditions, including sea temperature, salinity, velocity and sea surface height (SSH), are interpolated from outcomes of a global operational ocean forecasting system based on FIO-COM10 with a horizontal resolution of $\mathrm{1}/\mathrm{10}$^∘ (Qiao et al., 2019; Shi et al., 2018; Sun et al., 2020). As a result, the new $\mathrm{1}/\mathrm{32}$^∘ model enters a steady state after a short spin-up period of about 1 year (Fig. A1). The atmospheric forcing is from the Global Forecast System (GFS) of the National Centers for Environmental Prediction (NCEP) (https://www.nco.ncep.noaa.gov/pmb/products/gfs/, last access: 17 March 2023) and has a horizontal resolution of $\mathrm{1}/\mathrm{4}$^∘. The ocean surface heat and momentum fluxes are calculated via the bulk formula (Large and Yeager, 2004), and wind stress is calculated using relative wind.

2.2 Numerical experimental design

To investigate to what extent surface waves and tides contribute to the newly established FIO-COM32 model and taking the extremely expensive computational cost into consideration, two numerical experiments are designed. In EXP1, only the OGCM is active, neither the ocean surface waves nor the tides are included, and the simulated period is from 1 June 2016 to 31 December 2019. In EXP2, wave–tide–circulation coupling is enabled, which means both B_V and tidal currents are activated. Since the online coupled surface wave model will almost double the computational cost compared with the ocean-only counterpart, in this paper the $\mathrm{1}/\mathrm{32}$^∘ numerical experiments uses the B_V data calculated from a global $\mathrm{1}/\mathrm{4}$^∘ online coupled FIO-COM for the exactly same time period. Here we mainly focus on the effect of B_V in the new $\mathrm{1}/\mathrm{32}$^∘ ocean model, and due to the large-scale characteristics of the surface wave simulations and the extremely expensive computational costs, we turned off the global $\mathrm{1}/\mathrm{32}$^∘ online surface wave model directly and tested the modulation of ocean currents on surface waves near the end of this paper (Fig. 16).

EXP1 starts from outcomes of the global FIO-COM10 forecast system on 1 June 2016. EXP2 branches off from EXP1 on 1 July 2017. The data analysis focuses on the period from 1 January 2018 to 31 December 2019. The model output frequency of three-dimensional variables is daily, and the output frequency of SSH and steric SSH is hourly. An additional two experiments with a horizontal resolution of $\mathrm{1}/\mathrm{10}$^∘ are also conducted to investigate the effects of model resolution. The settings of the two $\mathrm{1}/\mathrm{10}$^∘ models are kept identical with those of EXP1 and EXP2, respectively; hence, they are named EXP1Low and EXP2Low, respectively. As has been mentioned in the model description, neither sea surface temperature nor salinity is restored to observational data in all of the numerical experiments, and the model drift during the experiment period is examined to ensure that the model maintains a reasonable evolution (Fig. A2). Since the time span of the numerical experiments is not long (3.5 years), the model drift is generally quite small.

Another two numerical experiments, OGCM + B_V and OGCM + tide, should be conducted to clearly identify the exact effects of surface waves and tides, respectively. However, considering the high computational costs, we only perform EXP1 and EXP2.

In order to understand how the model resolution affects the simulated barotropic tide, simulation results of the two global barotropic tide models with horizontal resolutions of $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ are compared. Ocean temperature and salinity in the global barotropic tide models are set as uniform values and are kept as constants during model integrations. The parameterised topography drag scheme is needed in the barotropic tide models since the energy transfer mechanism from barotropic to baroclinic tide is missing. As has already been noted, a topography drag scheme (Jayne and St. Laurent, 2001) is incorporated into the momentum equation following previous work (Jayne and St. Laurent 2001; Egbert et al., 2004; Arbic et al., 2004), and the drag coefficient is best tuned individually for each numerical experiment. The model topography in this paper is different from our previous work (Xiao et al., 2016) in which the treatment of ice shelf topography of Antarctica is more realistic due to the use of the water column thickness of BEDMAP2 (Fretwell et al., 2013), the use of which would yield a smaller root-mean-square error (RMSE).

3 Results

3.1 Effects of horizontal resolution

In this section, the surface EKE and SSH values simulated by the $\mathrm{1}/\mathrm{32}$^∘ models are compared with those of $\mathrm{1}/\mathrm{10}$^∘ models to investigate the sensitivity of different horizontal resolution.

One of the most noteworthy elements of these simulations is that the EKE increases significantly as the model horizontal resolution increases. As CX17 pointed out, this increase is due to smaller effective horizontal viscosity and more active mesoscale and sub-mesoscale motions of higher-resolution models. The simulated surface EKE and satellite observations are compared here to investigate the effects of improved horizontal resolution in FIO-COM32. The surface EKE is calculated as follows:

$$\begin{array}{}\text{(2)}& \mathrm{EKE}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle {u^{\prime }}^{\mathrm{2}}+{v^{\prime }}^{\mathrm{2}}\rangle ,\end{array}$$

where u^′ and v^′ are anomalous zonal and meridional components of surface geostrophic velocity, respectively, which are calculated from sea level anomalies referenced to the mean sea level of 2018 and 2019. The bracket indicates the time averaging of the model years of 2018 and 2019 or the spatial averaging of the global ocean. Satellite-derived surface EKE is calculated using absolute dynamic topography data of Ssalto/Duacs altimeter products produced by the Copernicus Marine and Environment Monitoring service (CMEMS, http://marine.copernicus.eu, last access: 17 March 2023). Generally, increasing the model resolution from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘ significantly improves the agreement of the simulated EKE and satellite observations (Fig. 1). Compared with the $\mathrm{1}/\mathrm{10}$^∘ model, the $\mathrm{1}/\mathrm{32}$^∘ model shows greatly enhanced EKE almost everywhere in the global ocean, especially in the regions of subtropics, mid-latitudes and western boundary current system regions, where surface EKE in the $\mathrm{1}/\mathrm{10}$^∘ model is too weak against satellite data. This indicates the horizontal resolution of $\mathrm{1}/\mathrm{10}$^∘ is insufficient to properly resolve mesoscale and sub-mesoscale motions, while the simulations of the $\mathrm{1}/\mathrm{32}$^∘ model are much improved. For the Kuroshio and Gulf Stream extensions, the EKE amplitudes and spatial distribution structures, such as the separation point and the extension range, are enhanced and agree with the satellite observations much better in the $\mathrm{1}/\mathrm{32}$^∘ model. The calculated global mean EKE time series (Fig. 2a) also show improved agreement for the $\mathrm{1}/\mathrm{32}$^∘ model against the satellite observations, while the $\mathrm{1}/\mathrm{32}$^∘ model has higher global mean EKE values. In order to understand how does the small-scale motions resolved by the $\mathrm{1}/\mathrm{32}$^∘ model contribute to the enhanced EKE level, a spatial filter with radius of ∼ 50 km is applied to separate the large-scale (mesoscale or larger) and small-scale (some sub-mesoscale) motions. As a result, the output of SSH of the $\mathrm{1}/\mathrm{32}$^∘ model is low-pass filtered, the corresponding low-pass EKE (Fig. 2b) is calculated and the residuals of total EKE are regarded as small-scale contributions (Fig. 2c). Both the low-pass model EKE and CMEMS are on the same $\mathrm{1}/\mathrm{4}$^∘ grid, which makes them more comparable. The time series of the low-pass EKE of the $\mathrm{1}/\mathrm{32}$^∘ model (the thick solid red line in Fig. 2a) is significantly lower than that of original $\mathrm{1}/\mathrm{32}$^∘ model and also lower than that of CMEMS, which indicates the sub-mesoscale motions (such as small-scale eddies and fronts) have important contributions to the enhanced EKE level of the $\mathrm{1}/\mathrm{32}$^∘ model (Fig. 2c). On the other hand, the low-pass EKE of the $\mathrm{1}/\mathrm{32}$^∘ model is still much more energetic than that of $\mathrm{1}/\mathrm{10}$^∘ model. This indicates that the mesoscale eddy fields are more energetic in the $\mathrm{1}/\mathrm{32}$^∘ model than in the $\mathrm{1}/\mathrm{10}$^∘ model. The calculated RMSE of the 2-year-averaged EKE between model simulations and satellite observations decreases from 63.5 cm² s⁻² in the $\mathrm{1}/\mathrm{10}$^∘ model to 31.7 cm² s⁻² in the $\mathrm{1}/\mathrm{32}$^∘ model. The RMSE is calculated for regions that have water depths exceeding 1000 m and that are located between 65^∘ S and 65^∘ N.

Figure 1 Eddy kinetic energy (EKE) of CMEMS from all merged and gridded satellite data (a) and the FIO-COM10 (EXP1Low) (b) and FIO-COM32 (EXP1) (c) model results.

Figure 2 Time series of globally averaged EKE (a), EKE of the $\mathrm{1}/\mathrm{32}$^∘ model filtered onto a $\mathrm{1}/\mathrm{4}$^∘ grid by a low-pass spatial filter with a radius of ∼ 50 km (b) and the EKE residual as a high-pass component (c).

Both the mean SSH and SSH standard deviation (SD) of the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models are compared in Fig. 3. The large-scale circulation patterns are reasonably simulated in both the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models (Fig. 3a, c, e). For the detailed simulation structures, such as pathways and separation points of the west boundary currents, the $\mathrm{1}/\mathrm{32}$^∘ model shows significant improvements over those of the $\mathrm{1}/\mathrm{10}$^∘ model. As the SSH SD comparisons show, the $\mathrm{1}/\mathrm{32}$^∘ model is more energetic than the $\mathrm{1}/\mathrm{10}$^∘ model in the well-known mesoscale active regions (Fig. 3b, d, f), which has already been discussed in the EKE analysis. In the Kuroshio region, the $\mathrm{1}/\mathrm{10}$^∘ model shows over-concentrated energy near the separating point, and the Kuroshio Extension is too weak to penetrate the interior of the Pacific Ocean, but these characteristics are much improved in the $\mathrm{1}/\mathrm{32}$^∘ model.

Figure 3 Mean SSH (a, c, e) and SSH SD (b, d, f) based on altimeter observations of CMEMS (a, b), FIO-COM10 (c, d) and FIO-COM32 (e, f) simulations.

To better understand the effect of model resolution on the large-scale circulations, we investigated the simulated main paths of the Kuroshio and Gulf Stream. The annual mean SSH of $\mathrm{1}/\mathrm{10}$ (EXP2Low) and $\mathrm{1}/\mathrm{32}$^∘ (EXP2) models of 2019 are compared with the corresponding altimeter observations of CMEMS. In this paper, we propose the integrated circulation route error (ICRE) as a quantitative criteria to assess the quality of the simulated paths of the Kuroshio and Gulf Stream. The ICRE is calculated as follows. First, the path of the western boundary current is defined as the contour edge of SSH at level of 0.16 m. Note that the regionally averaged value is extracted from the SSH before the ICRE analysis is conducted to ensure the SSH of the model and observation are comparable. Second, the ICRE is calculated as the total integrated misfit area between the contour edge of the model and CMEMS, and the regions of the ICRE diagnostics are defined as those in Figs. 4 and 5.

Figure 4 Mean SSH of the Kuroshio region. The route of the Kuroshio is represented by the contour line at level of 0.16 m. CMEMS (a), the $\mathrm{1}/\mathrm{10}$^∘ model (EXP2Low) (b) and its circulation route errors (d), and the $\mathrm{1}/\mathrm{32}$^∘ model (EXP2) (c) and its circulation route errors (e) are shown. Note that the regionally averaged value is extracted from the SSH before the ICRE analysis is conducted to ensure the SSH of the model and observation are comparable.

Figure 5 The same as Fig. 4 but for the Gulf Stream region. Note that the ICRE of the Gulf stream recirculation is masked out to focus on the main path of Gulf Stream.

As the resolution increases from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘, the simulated Kuroshio path is significantly improved (Fig. 4), the ICRE is decreased from 3.01 × 10¹¹ to 1.73 × 10¹¹ m². The most notable improvement is that the $\mathrm{1}/\mathrm{32}$^∘ model is able to reproduce the Kuroshio's large meanders, while the $\mathrm{1}/\mathrm{10}$^∘ model cannot. The ICRE of the Gulf Stream is decreased from 4.85 × 10¹¹ m² in the $\mathrm{1}/\mathrm{10}$^∘ model to 3.27 × 10¹¹ m² in the $\mathrm{1}/\mathrm{32}$^∘ model. We should note that for the comparisons of the Gulf Stream region, the recirculation part of the ICRE is masked out to focus on the simulation of main path of the Gulf Stream. Compared to the $\mathrm{1}/\mathrm{32}$^∘ model, the $\mathrm{1}/\mathrm{10}$^∘ model is not able to reproduce the deep penetration of the Gulf Stream into the Atlantic Ocean (Fig. 5).

Figures 6 and 7 show the simulated relative vorticity of sea surface current in the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models. They clearly show that as the horizontal resolution increases, the model can resolve many more sub-mesoscale and mesoscale motions. The $\mathrm{1}/\mathrm{32}$^∘ model also clearly shows more significant seasonal variations over the $\mathrm{1}/\mathrm{10}$^∘ model. This is because the mixed-layer instabilities are substantially enhanced in the $\mathrm{1}/\mathrm{32}$^∘ model, especially during winter; this phenomenon is consistent with previous publications (e.g. Sasaki et al., 2017; CX17; Dong et al., 2020b). In the higher-resolution model, the small-scale eddy activities are greatly enhanced during winter compared to those of summer, while this is not as obvious in $\mathrm{1}/\mathrm{10}$^∘ model. Stokes drift is the main mechanism that “selectively amplifies the vertical vorticity component as streamwise vorticity” according to Teixeira and Belcher (2002). Since the horizontal scales of the resolved fronts become much smaller in $\mathrm{1}/\mathrm{32}$^∘ model, the potential effects of explicit Stokes shear force can be estimated and compared among the models with different horizontal resolutions. The estimation is based on the parameter ϵ proposed by McWilliams and Fox-Kemper (2013):

$$\mathit{ϵ}={\displaystyle \frac{V^{\mathrm{s}}}{f{H}^{\mathrm{s}}}}{\displaystyle \frac{H}{L}},$$

where V^s is the Stokes drift characteristic velocity; f is the Coriolis parameter; H and L are the depth and width of the front, respectively; and H^s is the Stokes drift decay depth. The estimation is made for a typical section in the Kuroshio Extension region, which is marked as a black line in Fig. 6a. Figures A3 and A4 show the evolution of relative vorticity and temperature along a typical section in the Kuroshio Extension region. The typical horizontal scale of resolved fronts in the $\mathrm{1}/\mathrm{10}$^∘ model is about 60 km, while the horizontal scale in the $\mathrm{1}/\mathrm{32}$^∘ model is about 3 times smaller (∼ 20 km). The vertical scales of the resolved fronts of both the $\mathrm{1}/\mathrm{10}$ model and the $\mathrm{1}/\mathrm{32}$^∘ model are quite similar, the summer (winter) has a vertical scale of about 30 m (300 m). The dimensionless parameter representing the wave parameters in ϵ is $\frac{V^{\mathrm{s}}}{f{H}^{\mathrm{s}}}$, which is typically O(10–100) (Suzuki et al., 2016). The major difference in the estimations between the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models is the frontal aspect ratio in ϵ, which is $\frac{H}{L}$. During summer, the frontal aspect ratio is estimated to be ∼ $\mathrm{1}/\mathrm{1000}$ and ∼ $\mathrm{3}/\mathrm{1000}$ in the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models, respectively, which yield an estimated ϵ of ∼ $\mathrm{1}/\mathrm{10}$ and ∼ $\mathrm{3}/\mathrm{10}$, respectively. During winter the frontal aspect ratio is estimated to be ∼ $\mathrm{1}/\mathrm{200}$ and ∼ $\mathrm{3}/\mathrm{200}$ in the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models, respectively, yielding an estimated ϵ of ∼ $\mathrm{1}/\mathrm{2}$ and ∼ $\mathrm{3}/\mathrm{2}$, respectively. Hence, it can be concluded that as the horizontal resolution increased from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘, the critical resolution may be reached that the explicit Stokes shear force is not to be neglected (at least during winter). In this paper, the B_V value is included and could be treated as a “bulk” mixing term accounting for these wave–turbulence interactions. The effect of the explicit implementation of surface-wave-induced forces in the $\mathrm{1}/\mathrm{32}$^∘ model remains to be explored in the future.

Figure 6 Relative vorticity ($\frac{\mathit{\zeta }}{f}$) of the sea surface current in the $\mathrm{1}/\mathrm{10}$ (a, c) and $\mathrm{1}/\mathrm{32}$^∘ (b, d) model results. Snapshots of summer (1 September 2019) of the Kuroshio region (a, b) and the Gulf Stream region (c, d).

Figure 7 The same as Fig. 6 but for winter (1 March 2019).

Reasonable global tide simulation is an important prerequisite for tides and OGCMs to be able to be coupled, since tuning global barotropic tide can also yield better model topography settings, especially as some numerically unstable topography features can be fixed through this practice. The accuracy of simulated global barotropic tide is quite sensitive to the model resolution (Egbert et al., 2004). Here we show the effects of improving horizontal resolutions from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘ in the global barotropic tide simulations of FIO-COM. We also present the simulated global tide in EXP2, which is a global baroclinic tide model.

For the $\mathrm{1}/\mathrm{10}$^∘ barotropic model, the tidal amplitude is smaller than TPXO9 in some regions, such as west of Panama and in the northeastern Pacific, while the tidal amplitude is stronger than TPXO9 in the Labrador Sea, Sea of Okhotsk and Andaman Sea (Fig. 8b, c). These biases are obviously reduced in the $\mathrm{1}/\mathrm{32}$^∘ barotropic model (Fig. 8d), which indicates that the tidal dissipation in the regions of complex terrain is resolved better as the horizontal resolution is increased. The RMSE distributions show that the tidal accuracy of the $\mathrm{1}/\mathrm{32}$^∘ barotropic tide model is improved almost everywhere over that of the $\mathrm{1}/\mathrm{10}$^∘ barotropic model (Fig. 8e). The statistics show that for the M₂ constituent, the globally averaged RMSE of tidal elevation is decreased from 9.65 to 8.06 cm, i.e. reduced by 16.4 %. Here, the tidal elevation RMSE is calculated for regions with water depth exceeding 1000 m and located between 65^∘ S and 65^∘ N following previous publications (Arbic et al., 2004; Egbert et al., 2004).

Figure 8 M₂ co-tidal charts of TPXO9 (a), barotropic tide models of $\mathrm{1}/\mathrm{10}$^∘ (b) and $\mathrm{1}/\mathrm{32}$^∘ (d) resolution with best-tuned topographic drag parameterisation, a baroclinic tide model of $\mathrm{1}/\mathrm{32}$^∘ resolution (EXP2) (f) without topographic drag parameterisation, and the RMSEs of the three models (c, e and g, respectively).

The simulated global tide of EXP2, which is a global baroclinic tide model, is shown in Fig. 8f. The overall pattern of the M₂ global tide agrees well with the TPXO9, and most of the amphidromic points in the open basins are reproduced reasonably. The topography drag scheme is turned off in EXP2. As a result, the simulated tidal amplitude is significantly larger than that of the barotropic experiments (Fig. 8b, d) with topography drag optimally tuned, especially in the Atlantic Ocean and the eastern Pacific Ocean. The globally averaged tidal elevation RMSE is 16.1 cm (almost double the values in optimally tuned experiments). While the tide in the western and central Pacific Ocean and most of Indian Ocean agrees well with that of TPXO9, the amplitude bias and RMSE are much smaller than those of Atlantic Ocean. This result may indicate that a considerable portion of tidal energy conversion in the western Pacific region may be explicitly resolved in EXP2. As a result, the tidal dissipation parameterisation for EXP2 must be able to distinguish the unresolved and resolved tidal energy dissipation and conversion. To our knowledge, this kind of tidal dissipation scheme still does not exist in the present literature and is a daunting challenge for the ocean model scientific community. Development of such a tidal dissipation scheme is beyond the scope of this paper. Hence, we choose to disable the topographic drag scheme used in the barotropic runs to not have to tune it, which is aligned with the configuration of LLC4320 of MITgcm (Arbic et al., 2018). Although this leads to reasonably anticipated larger errors in simulated global tide of EXP2, we believe that in the current stage, the representation of tide–circulation coupled processes will benefit from this strategy. On the other hand, the barotropic tide is filtered out in the wavenumber spectrum analysis to avoid the direct impacts of the barotropic tide.

In order to investigate whether the tide of the $\mathrm{1}/\mathrm{32}$^∘ model could interact more strongly with the circulation than the $\mathrm{1}/\mathrm{10}$^∘ model, the incoherent tide is calculated and compared between the $\mathrm{1}/\mathrm{10}$ and $\mathrm{1}/\mathrm{32}$^∘ models (EXP2) in Fig. 9. First, the coherent tidal constants of eight main tidal constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁ and Q₁) are obtained by applying harmonic analysis to the hourly SSH model outputs. Second, the incoherent signals are obtained by extracting the predicted coherent signals of the eight main tidal constituents. Finally, a Butterworth 10th-order band-pass filter with semi-diurnal (1.73–2.13 cycles per day) bands is adopted to calculate the incoherent tide time series. The incoherent tide amplitude is defined as the standard deviation of the incoherent tide time series. The incoherent tide amplitude of the $\mathrm{1}/\mathrm{32}$^∘ model is obviously stronger than that of $\mathrm{1}/\mathrm{10}$^∘ model. The increased incoherent tide amplitude in the $\mathrm{1}/\mathrm{32}$^∘ model should be attributed to the increased eddy activities as indicated in the EKE comparisons (Fig. 1).

Figure 9 Incoherent tide amplitude of semi-diurnal tidal band of the $\mathrm{1}/\mathrm{10}$ (a) and $\mathrm{1}/\mathrm{32}$^∘ (b) models.

3.2 Effects of surface-wave–tide–circulation coupling

Prior to examining the effects of surface-wave-induced mixing in the newly established FIO-COM32 model, the accuracy of the ocean wave model is evaluated. Figure 10 shows the comparisons of satellite data from Jason-3 and model-simulated significant wave height (SWH). From the comparison, the model results are interpolated onto the satellite ground tracks. Both the distributions and seasonal cycle can be reproduced well by the global $\mathrm{1}/\mathrm{4}$^∘ resolution surface wave model. The RMSE of the simulated SWH is calculated against the Jason-3 along-track geophysical data records. The Jason-3 data are produced and distributed by AVISO (Archiving, Validation, and Interpolation of Satellite Oceanographic data, http://www.aviso.altimetry.fr/en/home.html, last access: 17 March 2023). We use the low-pass-filtered Jason-3 data for the comparison, and the global mean RMSE of the simulated SWH in 2018 and 2019 is 0.57 m.

Figure 10 Along-track seasonal mean significant wave height of Jason-3 (a, b) and EXP2 (c, d) data. Note that the surface wave model of EXP2 is offline and has a horizontal resolution of $\mathrm{1}/\mathrm{4}$^∘.

The ocean upper mixed layer is a crucial layer located between the ocean interior and atmosphere. The MLD determines the heat and momentum content of the upper-ocean boundary layer and is a key factor that reflects the ability of an ocean model. Figure 11 shows the derived summer MLD based on Argo observations (Fig. 11a) and EXP1 (Fig. 11b) and EXP2 (Fig. 11c) simulations. The summer is the mean of July, August and September (JAS) for the Northern Hemisphere and the mean of January, February and March (JFM) for the Southern Hemisphere. The MLD is defined as the depth at which potential density increases by 0.125 kg m⁻³ from the corresponding surface values. The model results are interpolated onto the Argo profiles in both space and time.

Figure 11 Summer (JAS and JFM for the Northern Hemisphere and Southern Hemisphere, respectively) MLD based on Argo observations (a) and EXP1 (b) and EXP2 (c) simulations. The difference in MLD between EXP1 and Argo and between EXP2 and EXP1 and the zonal mean MLD are shown in (d), (e), and (f), respectively.

Both EXP1 and EXP2 are able to reproduce the general patterns of summer MLD distributions. However, the MLD of EXP1 without B_V is shallower than that of Argo observations in most of the extratropical regions (Fig. 11d). These shallow biases are significantly alleviated in EXP2 (Fig. 11e), which is due to the inclusion of B_V (as previous work has shown). The comparison of the zonal mean MLD shows great improvement in EXP2 over EXP1. EXP2 fits quite well with the observations, except in the Antarctic Circumpolar Current regions. The mean error in MLD at mid-latitudes (10–40^∘ N/S) decreases from −4.8 (EXP1) to −0.6 m (EXP2). Although the resolution is increased to $\mathrm{1}/\mathrm{32}$^∘, the B_V is still an important contributor that improves the summer MLD simulations. The impact of B_V in the $\mathrm{1}/\mathrm{10}$^∘ model is quite similar to that of the $\mathrm{1}/\mathrm{32}$^∘ model (not shown), which means increasing horizontal resolution alone does not solve the long-standing MLD shallow biases in summer. As B_V deepens the simulated summer MLD, we further investigate whether it has a detectable effect on the length scale of the surface mixed-layer instability (MLI) (L_sml). We follow Dong et al. (2020b) to calculate the L_sml, which is defined as

$$L_{\mathrm{sml}}=\mathrm{6}{\displaystyle \frac{N_{\mathrm{ML}}{H}_{\mathrm{ML}}}{f}},$$

where the N_ML is the depth-averaged buoyancy frequency of the mixed layer and H_ML is the MLD determined by the relative variance methods (Huang et al., 2018). As shown in Fig. 12, compared with the Argo data, the L_sml of EXP1 with the KPP vertical mixing scheme has a negative bias in most regions similar to that of Dong et al. (2020b), and the L_sml of EXP2 with B_V is improved and conforms better with that of the Argo data (Fig. 12d).

Figure 12 Length scale of summer (JAS and JFM for the Northern Hemisphere and Southern Hemisphere, respectively) surface mixed-layer instability (L_sml) of Argo observations (a). The differences in L_sml between EXP1 and Argo and between EXP2 and EXP1 and the zonal mean L_sml are shown in (b), (c) and (d), respectively. The mean error in L_sml decreases from −0.96 km in EXP1 to 0.05 km in EXP2.

As the model resolution is increased and global tide is explicitly included in EXP2, the model is able to simulate global internal tide (IT), which is activated when tidal currents flow over rough topography. Since a large portion of barotropic-to-baroclinic energy transfer can be explicitly resolved in EXP2, we disable the parameterised topography drag that has been used in the barotropic model. Figure 13 compares the global internal tide fields of satellite MIOST (Multivariate Inversion of Ocean Surface Topography) (Ubelmann et al., 2022) and model simulations. We apply the radial filter with a radius of 4^∘ as a high-pass filter to remove the barotropic tide signals in the open ocean in the model SSH output. Then, harmonic analysis of 1-year high-pass SSH time series is performed to obtain the simulated M₂ internal tide amplitude. As shown in Fig. 13, for both the M₂ and K₁ constituents, the MIOST and EXP2 agree well with each other in terms of the positions of the IT “hot spot” generation sites and the long-range propagation patterns. For example, the well-known generation site of M₂ IT near Hawaii and world's strongest K₁ IT near Luzon Strait are reproduced well in the model. The model simulations tend to be more energetic than the MIOST data, and whether this bias is due to the mesoscale contamination of the satellite data in the MIOST or insufficient IT-related energy dissipation in present model remains to be explored in the future. The globally averaged bias of amplitude of M₂ internal tide at the sea surface is 0.7 cm. Some published global baroclinic tide models such as HYCOM (Arbic et al., 2012; Shriver et al., 2012) show a more conformed internal tide amplitude when compared to the satellite observations than this paper. This may due to the fact that the topographic drag used in their model damps both the amplitude of barotropic and internal tide, while we do not use this kind of dissipation scheme in EXP2.

Figure 13 M₂ internal tide amplitude at the surface from MIOST (a) and EXP2 (b); (c) and (d) are the same as (a) and (b) but for K₁.

Previous studies (SK12, CX17) reported that the high-resolution model-simulated SSH wavenumber spectral slope in the 70–250 km mesoscale domain is quite different from that of satellite observations. In the previous models, the slopes range between −5 and −4 in both low- and high-latitude regions. While the satellite observations show strong latitudinal variability, with a much flatter slope of about −1 in tropical ocean that increases to about −5 and −4 in high-latitude and high-EKE regions. CX21 showed that internal tide-induced SSH signals are a critical factor influencing the wavenumber spectral slopes. With tide included, the wavenumber spectral slopes in the 70–250 km mesoscale range of their $\mathrm{1}/\mathrm{50}$^∘ Atlantic regional ocean model show much better agreement with those of satellite observations than the model results without tide. Figure 13 shows that the internal tide-induced sea surface undulations with spatial scale of tens to hundreds kilometres can reach several centimetres, and thus there is a need to explore the different responses of wavenumber spectral slopes in EXP1 and EXP2.

The SSH wavenumber spectral slope in the 70–250 km range is calculated as follows. First, the model SSH data are interpolated onto Jason-3 ground tracks. Second, the spectral slopes are calculated on a 1^∘ × 1^∘ grid for the global ocean. For each grid point, tracks falling into a 10^∘ × 10^∘ super-domain centred at the grid point are taken for spectral calculation. In addition, it also requires sub-tracks that are longer than 500 km and that have 90 % or more good-quality data. Prior to spectral analysis, the trend of each sub-track is removed, and a Hanning window is applied. The spectra is obtained by averaging all the SSH spectra obtained via a fast Fourier transform (FFT). Finally, the spectral slope is calculated in the 70–250 km mesoscale range, which is in line with Xu and Fu (2012).

Figure 14 SSH wavenumber spectrum slope of 70–250 km in the Jason-3 along-track data (a), EXP1 (b) and EXP2 (c). SSH wavenumber spectrum at three sites in regions of the Kuroshio (d), tropical Pacific Ocean (e) and Southern Ocean (f); the black lines denote the wavenumber at 70 and 250 km. Snapshots of SSH from EXP1 (g) and EXP2 (h) on 1 December 2018 in the region shown with the rectangle in (a). Note that the sign of the spectrum slope is reversed to make it positive in (a)–(c).

Figure 14a, b and c show SSH wavenumber spectral slopes of 70–250 km range of Jason-3 along-track filtered data for EXP1 and EXP2. Figure 14a resembles the results of Xu and Fu (2012), which shows the strong latitudinal variability in wavenumber spectral slopes in the satellite observations. Figure 14b shows the wavenumber spectral slopes of EXP1 have much weaker latitudinal variability, which resemble that of SK12 and CX17. The wavenumber spectral slopes of EXP2 show substantially improved agreement with those of satellite observations shown in Fig. 14a. The flattened wavenumber spectral slopes of EXP2 in low-EKE regions result from internal tide-induced SSH undulations that can be clearly observed in SSH snapshots (see Figs. 14g and h and the Hovmöller diagram in Fig. A5). These internal tide signals have spatial scales of tens to hundreds of kilometres and become nontrivial where the background eddy and circulation field is relatively weak. Figure 14d, e and f show wavenumber spectra at three representative sites. Three locations are chosen for representing typical conditions: (i) high EKE and active internal tides, (ii) typical equatorial regions, and (iii) high EKE but inactive internal tides. In the high-latitude and high-EKE regions, the spectra are less affected by the internal tide signals. However, in the low-latitude and low-EKE regions, the presence of internal tide substantially alters the spectra and yields much more realistic simulations. The textures of SSH in EXP1 and EXP2 shown in Fig. 14g and h manifest many differences that mainly result from the internal tide signals, and EXP2 would be much closer to reality. Therefore, for high-resolution global ocean models, tide–circulation coupling is important for more detailed and realistic SSH simulations. Although the slopes of EXP2 (Fig. 14c) are more in agreement with the observations than those of EXP1 (Fig. 14b), there are still some discrepancies, especially in the low-latitude regions of the Atlantic. In this region, the slope of EXP2 is apparently more flat than the observations, which may indicate that the internal tide is stronger here than normal. A more realistic internal tidal field may further improve the agreement between the model and observations.

Figure 15 SSH variance spectrum amplitude [m² (cpkm × cph)⁻¹] in horizontal wavenumber–frequency space estimated from 10^∘ × 10^∘ boxes centred at 26^∘ N, 138^∘ E (a, d, g), 0^∘ N/S, 175^∘ E (b, e, h), and 55^∘ S, 175^∘ E (c, f, i). The upper and middle panels are based on hourly outputs of EXP2 and EXP1, respectively, and the lower panels are based on the $\mathrm{1}/\mathrm{10}$^∘ EXP2Low outputs. Dashed white lines denote the inertial frequencies, and dashed black lines denote the tidal constituents of S₂, M₂, K₁, and O₁. Solid black curves denote the dispersion relations for inertia-gravity waves of the 1st (upper) and 10th (lower) vertical modes. Note that the colour scale is exponential.

To further explore how the tide–circulation coupling and the increased resolution affect the simulations of different processes (including the rotational processes of sub-mesoscale turbulence and divergent processes of inertial gravity waves, IGWs), the wavenumber–frequency spectra of SSH variance are calculated from the hourly outputs of EXP1, EXP2 and EXP2Low simulations. Focusing on the high-frequency sub-mesoscale processes, the wavenumber–frequency spectra are calculated every 30 d and average the results to obtain a spectral estimate. Before calculating the discrete Fourier transform (DFT), we remove linear trends and multiply the data using a three-dimensional Hanning window. Meanwhile, the dispersion curves of IGWs are estimated from the World Ocean Atlas 2013 climatological temperature–salinity profiles (Qiu et al., 2018).

With respect to sub-mesoscale IGWs, these are shown by clear discrete beams above the 10th vertical normal mode curve in EXP2 (Fig. 15a–c), especially in the subtropical gyre of the western Pacific, which is invisible in EXP1 (Fig. 15d–f). By comparing the SSH wavenumber–frequency spectrum of EXP1 and EXP2, it shows that SSH wavenumber–frequency spectra are more energetic after including tidal forcing in the near-inertial, diurnal and semi-diurnal tidal frequency bands and along the dispersion curves of discrete vertical modes of inertial-gravity waves. This can further explain how the substantially enhanced unbalanced motions (according to Chereskin et al., 2019, “unbalanced motions” refers to the ageostrophic motions, but here the term mainly refers to internal tides and inertia-gravity waves) are responsible for the better match of EXP2 to Jason-3 than EXP1 in Fig. 14. As shown in Fig. 15a–b, in the regions with active IGWs, the energy distribution extends continuously from the sub-mesoscale regime to the IGW regime, which may indicate a strong non-linear interaction between them. However, in the low-resolution experiments shown in Fig. 15g–h there is either a gap between the sub-mesoscale and IGW regimes or insufficient sub-mesoscale activities, which may limit the non-linear interaction between them. It is obvious that SSH variance density is more energetic in the sub-mesoscale regime as resolution increases (<50 km and below the dispersion relation curve associated with the 10th vertical normal mode).

4 Summary and discussions

In this paper, for the first time a global $\mathrm{1}/\mathrm{32}$^∘ surface-wave–tide–circulation coupled model of FIO-COM32 is developed and validated. Increasing the model horizontal resolution from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘ leads to significant improvements in the simulation of surface EKE, the main paths of the Kuroshio and Gulf Stream, and global tides. We proposed ICRE as a quantitative criteria to evaluate the simulated main paths of the Kuroshio and Gulf Stream. As the model horizontal resolution increases from $\mathrm{1}/\mathrm{10}$ to $\mathrm{1}/\mathrm{32}$^∘, the ICRE of the simulated Kuroshio (Gulf Stream) path is decreased from 3.01 × 10¹¹ (4.85 × 10¹¹) to 1.73 × 10¹¹ (3.27 × 10¹¹) m². The RMSE of the EKE between the model simulations and satellite observations decreased from 63.5 cm² s⁻² in the $\mathrm{1}/\mathrm{10}$^∘ model to 31.7 cm² s⁻² in the $\mathrm{1}/\mathrm{32}$^∘ model. The global barotropic M₂ RMSE decreased from 9.65 cm in the $\mathrm{1}/\mathrm{10}$^∘ model to 8.06 cm in the $\mathrm{1}/\mathrm{32}$^∘ model. Although the resolution has been increased to $\mathrm{1}/\mathrm{32}$^∘, the non-breaking wave-induced mixing (B_V) is still a key factor in improving the summer MLD simulation against the Argo observations, with the mean error in the mid-latitude summer MLD reduced from −4.8 m in the numerical experiment without B_V to −0.6 m with B_V. Internal tides can be explicitly simulated in this new model, and global comparisons of along-track SSH wavenumber spectral slopes in 70–250 km domain show that the internal tide-induced SSH undulations are critical to substantially improved agreement between the global $\mathrm{1}/\mathrm{32}$^∘ model and satellite data, indicating that tide is a generator of sub-mesoscale process.

In future work, the effects of sea surface currents (SSCs) on surface wave simulation need to be taken into consideration for a fully coupled wave–circulation model, as this phenomenon has been shown to have impact on surface wave simulation (Ardhuin et al., 2017). As has been stated in this paper, the online coupled surface-wave–ocean circulation version of the FIO-COM32 is turned off pro tempore and is ready to be implemented when the demands for greater computational resources are met in the near future, and at that time the effects of ocean currents on the surface waves will be explored in depth. Figure 16 shows an example of the effect of SSCs on the simulated significant wave height. The hourly SSC output from FIO-COM32 in the western Pacific is fed into the MASNUM wave model to test SSC effects on surface waves. It is clear that SSCs have a notable influence on the simulated significant wave height, especially around the Kuroshio. In MASNUM wave model, SSCs influence the surface wave simulations through advection, refraction and the energy source function due to interaction between surface waves and SSC shear (Yang et al., 2005).

Figure 16 (a) Sea surface current (SSC) snapshot of 8 October 2018. (b) The differences in significant wave height between MASNUM wave models with and without SSC effects.

Since the internal tide can be resolved in the FIO-COM32 model, finding a more robust representation of energy dissipation of internal tide propagation is still an open task that needs to be addressed in the future. A proper dissipation scheme of the internal tide and its adaption to the traditional viscosity schemes for OGCM is still a daunting challenge. We have conducted numerical experiments (not shown) with an inordinately large background vertical viscosity of 1.0 m² s⁻¹, and in these experiments the amplitude of the simulated internal tide is brought down to the comparable level of the MIOST observations. These preliminary explorations indicate that a proper scheme of vertical mixing considering the internal gravity wave dissipation in the tide–circulation coupled ocean model needs to be developed. At the same time, proper representation of the unresolved barotropic-to-baroclinic tidal conversion also needs to be considered. More realistic treatments of SAL will improve the accuracy of global barotropic tide models (Arbic et al., 2004).

A better prediction of ocean is the inexhaustible goal of the ocean model community. As the ocean, especially the upper ocean, plays a dominant role in the climate system, ocean model improvement can shed light on new generations of climate model development. Surface waves, tides and ocean circulations are often separately simulated by different numerical models. In this paper, we show that surface wave–tide–circulation coupling can dramatically improve simulations. Therefore, it is time for us to regard the ocean as a fully coupled dynamic system through the key conjunction of turbulence and develop fully coupled surface-wave–internal-wave–tide–circulation models for the next generation of ocean model development.

Appendix A

Since the time span of the numerical experiments is not long (3.5 years), the model drift is generally quite small. The model drift in temperature shows a strong seasonal cycle in the upper ocean and a warming trend in the sub-surface layer, and the maximum drift in temperature at the end is about 0.2^∘ (Fig. A2a). The upper 150 m becomes slightly salty, and the sub-surface layer shows freshening drift, while the salinity drift value is generally small, i.e. less than 0.02 PSU (Fig. A2b). The zonal mean temperature and salinity and their corresponding drifts are displayed in Figs. A2c–f, and the zonal mean drifts are moderate.

Figure A1 Time series of simulated total kinetic energy of EXP1. Figure A1 shows the total kinetic energy of EXP1 since the beginning of the numerical experiment. The $\mathrm{1}/\mathrm{32}$^∘ model starts from a $\mathrm{1}/\mathrm{10}$^∘ model state and enters a new steady state after a swift adjustment period of about 1 year. After the short spin-up period, the total kinetic energy reaches quasi-steady state with small fluctuations.

Figure A2 Model drift of temperature (a) and salinity (b), zonally averaged temperature (c), salinity (d), and the corresponding drift (e, f) at the end of EXP1. As the last day of model integration is 31 December 2019, the drifts in (e) and (f) are calculated using the simulation difference between December 2016 and December 2019.

Figure A3 Hovmöller diagram of surface relative vorticity along a section in the Kuroshio Extension region, marked as black line shown in Fig. 6a.

Figure A4 Temperature along the section in the Kuroshio Extension region of $\mathrm{1}/\mathrm{10}$ (a, c) and $\mathrm{1}/\mathrm{32}$^∘ (b, d) models during summer (a, b) and winter (c, d), respectively. Figures A3 and A4 show the evolution of relative vorticity and temperature along a typical section in the Kuroshio Extension region, which is shown in Fig. 6a. The typical horizontal scale of resolved fronts in the $\mathrm{1}/\mathrm{10}$^∘ model is about 60 km, and the horizontal scale in the $\mathrm{1}/\mathrm{32}$^∘ model is about 3 times smaller (∼ 20 km). The vertical scales of the resolved fronts of both the $\mathrm{1}/\mathrm{10}$ model and $\mathrm{1}/\mathrm{32}$^∘ model are quite similar, and the summer (winter) has a vertical scale of about 30 m (300 m).

Figure A5 Hovmöller diagram of steric SSH along the section of 135^∘ E. EXP2 with tide shows obvious signals of the propagation of internal tide. This figure, together with Fig. 14g–h, explains how the internal tide signature affects the simulated SSH.