Mesoscale Eddy-Induced Ocean Dynamic and Thermodynamic Anomalies in the North Pacific

Zhou, Jiang; Zhou, Guidi; Liu, Hailong; Li, Zhuhua; Cheng, Xuhua

doi:10.3389/fmars.2021.756918

ORIGINAL RESEARCH article

Front. Mar. Sci., 12 November 2021

Sec. Physical Oceanography

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.756918

Mesoscale Eddy-Induced Ocean Dynamic and Thermodynamic Anomalies in the North Pacific

$\r\nJiang Zhou,$ Jiang Zhou^1,2

Guidi Zhou^1,3*

Hailong Liu⁴

Zhuhua Li^1,3

Xuhua Cheng^1,3,5

¹College of Oceanography, Hohai University, Nanjing, China
²Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China
³Key Laboratory of Marine Hazards Forecasting, Ministry of Natural Resources (MNR), Hohai University, Nanjing, China
⁴State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China
⁵Southern Laboratory of Ocean Science and Engineering Guangdong Laboratory, Zhuhai, China

Oceanic mesoscale eddies are associated with large thermodynamic anomalies, yet so far they are most commonly studied in terms of surface temperature and in the sense of composite mean. Here we employ an objective eddy identification and tracking algorithm together with a novel matching and filling procedure to more thoroughly examine eddy-induced thermodynamic anomalies in the North Pacific, their relationship with eddy amplitude (SSH), and the percentage of variability they explain on various timescales from submonthly to interannual. The thermodynamic anomalies are investigated in terms of sea surface temperature (SST), isothermal layer depth (ITD), and upper ocean heat content (HCT). Most eddies are weak in amplitude and are associated with small thermodynamic anomalies. In the sense of composite mean, anticyclonic eddies are generally warm eddies with deeper isothermal layer and larger heat content, and the reverse is true for cyclonic eddies. A small fraction of eddies, most probably subsurface eddies, exhibits the opposite polarities. Linear relationships with eddy amplitude are found for each of the thermodynamic parameters but with different level of scatter and seasonality. HCT-amplitude relation scatters the least and has the smallest seasonal difference, ITD-amplitude relation has the largest scatter and seasonality, while SST-amplitude relation is in between. For the Kuroshio and Oyashio Extension region, the most eddy-rich region in the North Pacific, eddies are responsible for over 50% of the total SSH variability up to the intra-seasonal scale, and ITD and HCT variability up to interannual. Eddy-induced SST variability is the highest along the Oyashio Extension Front on the order of 40–60% on submonthly scales. These results highlight the role of mesoscale eddies in ocean thermodynamic variability and in air-sea interaction.

Introduction

Mesoscale eddies are ubiquitous in the global oceans. They are coherent rotating vortices on horizontal scales of O(100 km) and temporal scales of O(100 days). Surface eddies are characterized by sea level rise or fall that are associated with significant geostrophic velocity anomalies and therefore can be detected and tracked from sea level height maps (Morrow and Le Traon, 2012; Chelton, 2013; Petersen et al., 2013). Thanks to the advancement of satellite altimeter observations and objective eddy detection and tracking algorithms, we now have some knowledge about the spatial, statistical, and temporal distribution of some of the key dynamic and kinematic properties of eddies (Hwang et al., 2004; Itoh and Yasuda, 2010; Cheng et al., 2014). From the perspective of energetics, kinetic energy of eddies accounts for a considerable fraction of the total kinetic energy (e.g., Wyrtki et al., 1976; Richardson, 1983; Zhou and Cheng, 2021b), dynamically adjusted by interaction with the time-mean flow (Ferrari and Wunsch, 2009; Qiu and Chen, 2010; Waterman et al., 2011; Yang and Liang, 2016; Sérazin et al., 2018; Yang et al., 2018). In the North Pacific, the Kuroshio Extension and the Subtropical Counter Current are two of the most eddy-energetic regions with dynamic interactions between eddies and the time-mean flow. Other regions, such as along the eastern boundary and in the Gulf of Alaska, are also abundant with eddies, although there eddies are generally weaker (Cheng et al., 2014).

Besides dynamic anomalies like sea surface height and velocity, eddies are also associated with thermodynamic and biogeochemical perturbations at the surface and beneath, and in their atmospheric imprints. Satellite-based studies showed that on average, mesoscale eddies can cause temperature anomalies at the surface on the order of 0.5°C (Frenger et al., 2013; Faghmous et al., 2015; Chen et al., 2017). Dong et al. (2014) used Argo profiles to show that eddy-induced temperature anomalies extends to about 1,000 m below surface, maximized at the thermocline, having huge impacts on the upper ocean heat content. Remarkable eddy-signature in salinity and chlorophyll is also found (Chelton et al., 2011; Gaube et al., 2013; Wang et al., 2018; Zhao et al., 2021). It is also established that ocean mesoscale eddies can modulate the air temperature, wind velocity, and precipitation in the atmospheric boundary layer directly above (e.g., O’Neill, 2012; Frenger et al., 2013; Putrasahan et al., 2017; Sugimoto et al., 2017), and may even be able to change the free atmospheric circulation as well as atmospheric transient eddy activity (e.g., Ma et al., 2015, 2017; Zhou et al., 2015; Zhang et al., 2019, 2020; Liu X. et al., 2021).

However, our current understanding on eddy-induced dynamic and thermodynamic anomalies is still poor. First, most studies only examined temperature perturbations on the sea surface, while eddies’ impact on subsurface temperature in the upper ocean is mostly untouched. Lack of high-resolution subsurface observations is the main reason. Second, eddy-induced temperature anomalies are commonly studied in the sense of composite mean over a large number of eddies. Although some studies presented the statistical distribution (or at least a standard deviation besides the mean) of height anomalies, such distributions of temperature anomalies are rarely investigated. The relationship between dynamic and thermodynamic anomalies, moreover, is still to be explored. Third, aside from the absolute values, it remains unknown how important the eddy-induced perturbations are relative to the total variability of dynamic and thermodynamic quantities on various timescales. This problem is of particular importance in, e.g., midlatitude air-sea interaction. It is historically under debate as to what extent is the midlatitude surface temperature variability on interannual and shorter timescales driven by atmospheric forcing, and to what extent by ocean internal processes that can potentially force atmospheric responses. Recent studies started to highlight the ocean control on surface temperature (Smirnov et al., 2014; Bishop et al., 2017), based on autoregressive models linearized using either observed or model-simulated temperature and heat flux. These studies extend the conventional stochastic climate model proposed by Hasselmann (1976) and Frankignoul and Hasselmann (1977) by enabling ocean noise forcing, yet they still rely on statistical fitting, and the shortest timescale investigated is limited to monthly due to data resolution. Direct description of eddy-induced anomalies and multi-scale variability down to submonthly timescales is therefore in keen interest.

This paper is an attempt to tackle the above issues using direct eddy identification and matching approaches. It answers the following key questions: (a) how are eddy morphologic, kinematic, and thermodynamic parameters statistically related to eddy amplitude; (b) how much dynamic and thermodynamic variability do the eddy-induced anomalies explain. These questions are important to ocean dynamics and air-sea interaction, but have not been thoroughly investigated. Here by thermodynamic parameters, we consider both surface temperature and subsurface heat anomalies in terms of isothermal layer depth and upper ocean heat content. The rest of the paper is organized as follows: In section “Data and Methods” we introduce the data used in this work, and the eddy identification and tracking method. Section “Results” presents the results in three subsections: Statistical relationship of eddy morphologic and kinematic parameters with eddy amplitude (section “Detected Eddy Properties”); eddy-induced thermodynamic anomalies and their statistical relationship with the amplitude (section “Eddy-Induced Anomalies”); and the relative importance of eddy-induced variability on different timescales (section “Eddy-Induced Variabilities”). Section “Conclusion and Discussion” includes a summary of the results and a brief discussion.

Data and Methods

Data

In this study, we make use of daily sea surface height (SSH) and sea surface temperature (SST) data provided by the Estimating the Circulation and Climate of the Ocean, phase II, project (Menemenlis et al., 2008) on a spatial grid of 1/4°× 1/4°. The study domain is the North Pacific (120°E–90°W, 10°–60°N). Ocean sub-surface potential temperature and salinity from the same dataset with 3-day interval are used to estimate the isothermal layer depth and upper ocean heat content, which are then interpolated into daily records. There are 50 levels in the vertical from 5 m to 6150 m with the layer thickness ranging from 10 m below the surface to 456 m near the bottom. The time range of all data is 2008–2017 (10 years).

The reasons for using the ECCO2 reanalysis product are twofold. First, there lacks long-term Eulerian observations of ocean sub-surface temperature necessary for estimating isothermal layer depth and upper ocean heat content. Argo floats do a good job of measuring sub-surface temperature, but they drift around with the ocean currents. Second, this paper aims at examining both dynamic and thermodynamic anomalies of eddies, and relies on matching the temperature anomalies with eddies detected from SSH maps, thus physical consistency between SSH and the surface and sub-surface temperature is essential. Satellite observations of SSH and SST are based on different hardware platforms with different retrieval methodologies and therefore cannot guarantee physical consistency. Reanalysis datasets commonly employ certain data assimilation techniques that introduce unphysical jumps and artificial sources/sinks. Although also a reanalysis, the ECCO2 project estimates the ocean state by a least square fit of the Massachusetts Institute of Technology General Circulation Model (Marshall et al., 1997) to available satellite (e.g., altimeter for SSH, radiometer for SST) and in-situ data, including Argo floats. Using a Green function approach (Menemenlis et al., 2005), the fit is applied for a number of control parameters (e.g., the initial and boundary conditions, vertical viscosity, and bottom drag), and the model is integrated forward freely without observational data injection, hence ensuring dynamic and thermodynamic consistency along with accuracy (Wunsch et al., 2009). The ECCO2 dataset has been employed to diagnose eddy dynamics and energetics in various parts of the world oceans (Fu, 2009; Chen et al., 2014; Qiu et al., 2017; Yang et al., 2017, 2018). The following study is based on anomalies defined as the deviation from the 10-year climatological annual cycle. For simplicity hereafter we drop the word “anomaly” in this sense and reserve it to indicate the eddy anomaly, i.e., the difference of a quantity inside and outside an eddy, as elaborated in section “Eddy-Induced Anomalies.”

Eddy Identification and Tracking

Over the past decades, a number of autonomous eddy detection algorithms has been developed, based on either observational or modeling data (Chelton et al., 2007, 2011; Doglioli et al., 2007; Chaigneau et al., 2008; Nencioli et al., 2010; Dong et al., 2011; Faghmous et al., 2015; Zhang C. et al., 2015). We use the objective autonomous eddy identification algorithm developed by Faghmous et al. (2015). This algorithm first detects eddy centers by finding local SSH extrema, and then defines the outermost closed contours around the extrema as eddy boundaries. It is free of expert parameters that must be determined a priori. This algorithm provides eddy boundary information, which is of particular importance to this study in calculating eddy-induced anomalies. It should be noted that as in previous studies (e.g., Cheng et al., 2014; Faghmous et al., 2015; Liu et al., 2020), no strict distinction between mesoscale eddies and long baroclinic Rossby waves are made, since the two are both manifested as SSH extrema and they often coincide (Pingree and Sinha, 2001; Oliveira and Polito, 2013; Polito and Sato, 2015). Eddy trajectories are found by searching for similar eddies between the adjacent timesteps within a search radius determined by the maximum eddy propagation speed—the phase speed of long baroclinic Rossby waves. To account for data complexity, here we allow the eddies to disappear for 1 day in its trajectory, thus not counting its re-appearance in the next day as a different eddy. This algorithm has been successfully used in the past to detect and track mesoscale eddies in the world oceans (Putrasahan et al., 2017; Ding et al., 2018; Liu et al., 2020).

Recently, eddies maximizing below the sea surface start to attract more research interests. These subsurface eddies are being observed in various parts of the global oceans (Pelland et al., 2013; Pegliasco et al., 2015; Zhang Z. et al., 2015, 2017; Song et al., 2019; Zhu et al., 2021), and simulated in ocean models (Chiang et al., 2015; Wang, 2017; Xu et al., 2019, 2020). Only some of the subsurface eddies exert imprints on SSH, particularly the intra-thermocline eddies whose isopycnals outcrop the sea surface, and some relatively shallow sub-thermocline eddies that make contact with the mixed layer on top which transfers the eddy-induced anomalies to the surface (Assassi et al., 2016). These surface signals of subsurface eddies, however, are expected to be weak, thus the detection quality of these eddies is compromised. Deep sub-thermocline eddies with the “lens” shape are associated with anomalies of the opposite sign above and below its core, and therefore have very little surface imprint. In this work we use sea surface perturbations to detect eddies, and thereby restrict ourselves to those eddies that are detectable from the surface. Deep sub-thermocline eddies are left for future research.

Results

Detected Eddy Properties

In the 3653 days of 2008–2017, an average number of 524 cyclonic and 506 anticyclonic eddies are detected per day in the North Pacific basin, which, respectively, fall into 110,615 and 106,140 trajectories. Figure 1A shows a snapshot of the SSH field and the detected eddy boundaries. It is clear that the eddies correspond well to the SSH contours. Eddies in the Kuroshio Extension (KE, 140–170°E, 30–40°N) region are largest in size and strongest in amplitude, followed by the Subtropical Counter Current (STCC, 125°E–155°W, 15–25°N) region and along the eastern boundary and the Aleutian Islands. The distribution of eddy lifespan averaged across 1°× 1° cells is show in Figure 1B. Long-living eddies are mostly generated in the eastern basin, whereas KE and STCC eddies are short. This is presumably due to the strong baroclinic and barotropic interaction with the background flow (Ferrari and Wunsch, 2009; Waterman et al., 2011; Waterman and Jayne, 2011; Bishop, 2013; Chen et al., 2014; Yang and Liang, 2016; Yang et al., 2017, 2018; Ji et al., 2018). The trajectories of eddies longer than 180 days are shown in Figure 1C, which reveals that eddies travel longer (shorter) distances near the equator (at higher latitudes), as a result of the poleward decreasing Rossby wave speed. Eddies in the subarctic North Pacific (in the Gulf of Alaska and along the Aleutian Islands) are long-living and quasi-stationary, exhibiting the shortest travel distance (Di Lorenzo et al., 2005; Rogachev et al., 2007; Saito et al., 2014; Rogachev and Shlyk, 2018, 2019). Active eddy generation is noticeable along the eastern boundary, where coastal Kelvin waves from the tropics are transformed into Rossby waves propagating westward away from the boundary (Clarke and Shi, 1991; Fu and Qiu, 2002). These above results agree quantitatively with the work by Cheng et al. (2014) based on altimeter data, which confirms the credibility of the eddy identification and tracking results.

FIGURE 1

Figure 1. (A) Snapshot of the SSH anomaly (cm) field on January 1, 2008 (shading). Anomaly is defined as the departure from the climatological annual cycle. Areas enclosed by the red (cyan) curve represents the identified anticyclonic (cyclonic) eddies. (B) Mean eddy lifespan (days) averaged within 1°×1° grid cells organized by the generation location of eddies. (C) Trajectories of eddies with lifespan longer than 180 days detected on January 1, 2008. Red (blue) curves represent the trajectories of anticyclonic (cyclonic) eddies. (D) Mean eddy amplitude (i.e., SSHA, cm) averaged within 1°×1° grid cells organized by the instantaneous location of eddies. Shown are the absolute values averaged across cyclonic and anticyclonic eddies.

We define for each eddy an equivalent circle having the same area as the eddy, the equivalent circle’s radius as the eddy radius, and then the eddy-induced SSH anomaly (SSHA) as the difference between the SSH extrema within the eddy boundary and the average SSH outside the eddy but within three radii. This definition is similar to, but not the same as, the one employed by Frenger et al. (2013), who took the inside value as the average within two eddy radii, and the ambient value as the average within three radii. Their definition of interior value includes anomalies outside the eddy boundary, and their ambient value includes areas inside. Cheng et al. (2014) used the same inside SSH as ours, but calculated the ambient SSH as the average along the boundary. This definition suffers from contamination by residual large-scale imprints surviving the removal of the climatological means (e.g., in the STCC in Figure 1A). Our definition arguably does a better job in accurately representing both the interior and ambient values. The eddy-induced SSHA is referred to as the eddy amplitude in this work. Figure 1D shows the spatial distribution of mean eddy amplitude, which agrees well with Cheng et al. (2014) in terms of pattern and magnitude. Eddy-induced anomalies of other quantities are defined in the same way. It should be noted that in this study, in order to assess eddy-induced variability on the smallest possible scales, we do not apply any filtering criterion on e.g., eddy size or lifetime, as commonly did in literature studies (Frenger et al., 2013, 2015; Cheng et al., 2014; Chen et al., 2017; Liu et al., 2020). This results in a tremendous number of weak eddies.

To have a better overview on the detected eddies, we plot in Figure 2 the relationship between the eddies’ SSHA and their morphologic and kinematic parameters. Some previous studies (e.g., Cheng et al., 2014) investigated statistical distributions of these parameters, yet their relationship with the amplitude is seldom examined. First, Figure 2A shows that on the first order, eddy amplitude and radius exhibit a linear relationship. Large scatter, however, is present especially for moderate-amplitude eddies with radius ranging from as small as a few tens of km to as large as 350 km. The strongest eddies (| SSHA| > 60 cm) scatter less, having radius of ∼100–200 km. The most abundant weak eddies (| SSHA| < 4 cm) are at most 150 km large. Next, in Figure 2B we present the distribution of eddy deformation ratio defined as the area of the deviation between the eddy boundary and its equivalent circle relative to the eddy area following Kurian et al. (2011). Most of the strong eddies are moderately deformed, while weak eddies span over a large scatter, from near-regular to severely distorted. As for lifespan, it is found that strong eddies, mostly in the KE and STCC according to Figure 1D, are basically shorter than about a season (Figure 2C), while weak eddies can be anything between transient (a few days) and quasi-permanent (a few years) depending on their location. The translation speed of eddies is shown in Figure 2D, normalized by the theoretical value of the phase speed of long baroclinic Rossby waves at the eddies’ latitude given by Qiu et al. (1997). Note that theoretical Rossby wave speed is slightly slower than observed, thus the ratio shown in this figure is slightly higher. Large scatter is noticed in the figure, with weak eddies ranging from very slow to very fast, while strong eddies tend to travel at a speed closer to the Rossby wave speed. In terms of travel direction defined as the azimuth of the end point of a track relative to the starting point, there exists a preferred westward propagation direction (Figure 2E), same as found previously (Pingree and Sinha, 2001; Palacios and Bograd, 2005; Chelton et al., 2011). However, it is particularly interesting to note that taking all eddies into account, the dominance of westward propagation is not as obvious as for long-lived eddies, and eddies can travel to virtually every direction (Figure 2F). To further examine this point, we plot in Figure 3A the spatial distribution of mean propagation direction of all eddies, and in Figure 3B only the ones persisting for at least 150 days. The generation location of the non-westward (especially the eastward) propagating eddies largely agrees with the location of strong eddies (Figure 1C). Since the strong eddies are mostly short-lived ones (Figure 2C), raising the threshold of lifespan filters out the strong eddies that propagate to the east or other directions. The association between strong eddies and eastward propagation indicates that the simple theory of westward propagation required by conservation of potential vorticity and the beta effect (Nof, 1981; Cushman-Roisen et al., 1990) breaks up in areas with strong advection by eastward background flow (Hughes, 1996; Holland and Mitchum, 2001; Fang and Morrow, 2003; Morrow et al., 2003).

FIGURE 2

Figure 2. Scatter plots showing the statistical relationship between eddy amplitude (cm) and (A) radius (km), (B) deformation ratio, (C) lifespan (days; note the logarithmic y-axis; amplitude is taken at the eddies’ half-life), (D) the ratio between travel speed defined as the track length divided by the lifespan and the theoretical speed of long baroclinic Rossby wave (note the logarithmic y-axis), and (E) travel direction defined as the azimuth of eddy termination location relative to the generation location. Red (blue) dots represent anticyclonic (cyclonic) eddies. The radial axis indicates eddy amplitude in cm. (F) Histogram of eddy travel direction. Red stairs represent the histogram over all eddies, and blue stars are for eddies persisting for longer than 150 days.

FIGURE 3

Figure 3. Mean eddy direction defined as the azimuth of eddy termination location relative to the generation location, averaged within 1°×1° grid cells, organized by the generation location. (A) Shows mean eddy direction across all eddies in a grid cell, and (B) show mean direction across eddies persisting for longer than 150 days.

Eddy-Induced Anomalies

We next investigate thermodynamic anomalies induced by mesoscale eddies in terms of SST, isothermal layer depth (ITD), and upper ocean heat content (HCT), along with SSH. Here the ITD is calculated as the depth where seawater temperature decreases by 0.2°C from the temperature at 10 m below the surface. This definition has been widely used in literature studies (de Boyer Montégut et al., 2004; Breugem et al., 2008; Pang et al., 2019). The HCT is defined over the upper 400 m as $\int_{- 400}^{0} c ρ T d z$ , where T is temperature, c is specific heat capacity, and ρ is density. c and ρ are calculated from temperature and salinity using the polynomial algorithm of Fofonoff and Millard (1983). The choose of the upper 400 m is in accordance with many literature works (Vivier et al., 2002; Doney et al., 2007; Song et al., 2014), based on the consideration that the upper 400 m is thermodynamically the most energetic layer of the ocean. It has been shown that about 52% of surface-intensified eddies extend to less than 400 m in the vertical (Liu et al., 2020). An alternative limit of 700 m does not change our findings qualitatively.

Eddy-induced anomalies to the quantities under study [i.e., SST anomaly (SSTA), ITD anomaly (ITDA) and HCT anomaly (HCTA)] are calculated in the same way as the SSHA (as differences between the extrema inside the eddy and mean value outside the eddy but within 3 radii). The spatial distribution of mean eddy-induced thermodynamic anomalies is shown in Figure 4. On average, eddies in the Kuroshio and Oyashio Extension region are most effective in changing the SST, with mean SSTA of more than 2°C (Figure 4A). The largest SSTA appears in the confluence region between the Kuroshio and the Oyashio, reflecting the vigorous poleward heat transport by mesoscale eddies (Sugimoto and Hanawa, 2011; Seo et al., 2014; Sugimoto et al., 2014; Masunaga et al., 2016), and along the KE and Oyashio Extension currents. This is also where the eddies are the strongest (Figure 1D). Eddy-induced ITDA is the largest in the northwest, especially along the Kuril-Kamchatka Trench and around the Kuril Straits where the average | ITDA| is about 300 m (Figure 4B) and extreme values concentrate (Figure 4C). Climatologically, the background ITD in subpolar western North Pacific is known to be quite large in winter (∼1,000–3,000 m, Figure 4B contours). This is because of the strong vertical temperature inversion driven by the combination of intense surface cooling and the maintenance of the shallow thermocline by heavy precipitation and Ekman suction (de Boyer Montégut et al., 2007). The subsurface water, which is well-mixed due to the deep intense tidal mixing there (Kawasaki and Hasumi, 2010; Tanaka et al., 2010; Yagi and Yasuda, 2013; Kawasaki et al., 2021), is thus kept relatively warm, protected from the surface cooling. In this work the ITD is defined as the depth where temperature is 0.2°C lower than the SST, which means that it must be searched for beneath the large body of well-mixed warm water and is hence very deep. Mesoscale eddies traveling through these areas advect water mass from somewhere else and can thus bring along quite large ITDAs, even though their amplitudes are small. Such deep climatological ITD is not present in summer because of the surface heating and the seasonal thermocline, and so are the extreme eddy-induced ITDAs. Besides the western subpolars, the strong eddies in the KE region are also characterized with large ITDAs. We also calculated eddy-induced anomalies of mixed layer depth, which is almost identical to ITDA except in the temperature inversion regions where the mixed layer is shallow (not shown). Eddy-induced HCTA, as indicated by Figure 4D, exhibits a similar pattern with the eddy amplitude, which is not surprising since theoretically the eddy SSH anomaly is linearly related to heat content anomaly due to the thermal expansion effect (e.g., Stephenson et al., 2013).

FIGURE 4

Figure 4. Mean eddy-induced (A) SSTA (°C) (B) ITDA (m) averaged within 1°×1° grid cells, organized by the instantaneous location of eddies. Shown are the absolute values averaged across cyclonic and anticyclonic eddies. In (B), contours show the climatological mean ITD averaged across 1993–2019. Magenta contours are for winter (DJF) at 200, 500, 1,000, and 2,000 m levels, and blue contours are for summer (JJA) at 50, 100, and 150 m levels. (C) Location of eddies with | ITDA| >200 m (cyan), 500 m (blue), 1,000 m (red), and 1,500 m (yellow). (D) Same as (A,B) but for HCTA (GJ/m²).

It is a good start to examine the statistical distribution of eddy anomalies in the sense of composite mean. To do this, we follow Frenger et al. (2013) and Chen et al. (2017) to normalize the eddies according to their respective radius, but do not rotate them, with no restriction on eddy size or lifetime whatsoever. The results are shown in Figure 5. The composite-mean SSHA and SSTA confirm that generally, anticyclonic (cyclonic) eddies are associated with positive (negative) SSHA and higher (lower) temperature, which is in qualitative agreement with previous studies (Frenger et al., 2013; Faghmous et al., 2015; Chen et al., 2017), but is much weaker in amplitude as a result of the inclusion of the large number of weak eddies. The composite-mean eddy ITDA and HCTA are consistent with SSTA, indicating warmer and deeper isothermal layer containing more heat associated with anticyclonic eddies and vice versa. On average, a typical eddy with 8 cm SSHA would induce SSTA of 0.5°C, ITDA of 7 m, and HCTA of 1 GJ/m². The eddy impact maximizes near the eddy center and is basically isotropic. Cyclonic and anticyclonic eddies show general symmetry.

FIGURE 5

Figure 5. Eddy-induced (A) SSHA (cm), (B) SSTA (°C), (C) ITDA (m), (D) HCTA (GJ/m²) averaged over all eddies. The upper (lower) row is for anticyclonic (cyclonic) eddies. Eddies are scaled relative to their respective radius. Shown are the distribution of the indicated quantities within 3 times the eddy radius in each direction.

On the other hand, histograms of the parameters (Figure 6) indicate that the vast majority of eddies are very weak (P[| SSHA| < 4 cm] = 81%), and are associated with weak SSTA (P[| SSTA| < 0.4°C] = 80%), ITDA (P[| ITDA| < 20 m] = 69%), and HCTA (P[| HCTA| < 1 GJ/m²] = 64%). This is a result of the absence of eddy screening in this work. Here, for ITDA, we excluded the eddies in the western subpolar regions because of their very weak amplitude and very large ITDA. It is shown that eddy-induced ITDA elsewhere exhibits distinct seasonality (Figures 6D,E). In summer about 99% of the ITDA is below ± 20 m, while in winter this percentage is 53%. Further, the statistical relationship between eddy amplitude (SSHA) and eddy-induced SSTA is illustrated in Figure 7A. The maximum SSTA can reach ∼ ± 6°C, and average SSTA for strong eddies with |SSHA| > 40 cm is on the order of ± 2°C. Despite some scatters, a statistically significant linear relationship is still visible (linear coefficient = 0.033°C/cm). Seasonal differences are also shown in Figures 7B,C. It could be inferred that in summer, the SSHA-SSTA relation is flatter, with a larger number of strong eddies associated with weak SSTA, indicating that the mesoscale eddies are less capable at changing the SST in summer. This is presumably due to the formation of the shallow seasonal thermocline, which leads to easier SST disturbance by other forcing terms like solar radiation and cloud cover. Moreover, it is noticeable that about 20% of the eddies, mostly weak ones, have inconsistent SSHA and SSTA, i.e., anticyclonic eddies being cold eddies and vice versa. These eddies are also known as abnormal eddies (Liu Y. et al., 2021). The presence of such abnormality is most probably due to subsurface eddies. Recently, Bashmachnikov et al. (2013) and Assassi et al. (2016) showed that surface eddies have SSTAs of the same sign as the SSHA, while for subsurface eddies that do have surface imprints, the two anomalies have the opposite sign. Since subsurface eddies maximizes beneath the surface, their surface signatures are expected to be weak. Other processes, e.g., baroclinic instability, surface flooding or mixing etc., may also generate abnormal eddies (Liu Y. et al., 2021).

FIGURE 6

Figure 6. Histograms of eddy-induced (A) SSHA (cm), (B) SSTA (°C), (C) overall ITDA (m), (D) ITDA in winter (DJF, m), (E) ITDA in summer (JJA, m), (F) HCTA (GJ/m²).

FIGURE 7

Figure 7. Scatter plots showing statistical relationship between eddy amplitude (SSHA, cm) and (A–C) SSTA (°C), (D–F) ITDA (m), (G–I) HCTA (GJ/m²) over all time (1st column), winter (DJF, second column), and summer (JJA, last column). The magenta (green) line represents the average over all anticyclonic (cyclonic) eddies. The black line is the linear relationship calculated by least square fitting.

The eddy-induced ITDA, as shown in Figure 7D, can be as large as over 1,800 m, even though the western subpolar regions are excluded. The extreme values are, of course, scarce. For clarity, Figure 7D only shows | ITDA| ≤ 500 m. A closer look on winter and summer results (Figures 7E,F) reveals that most of the large-ITDA eddies are found in winter, while summer ITDAs are very weak (<50 m), regardless of eddy amplitude. The wintertime SSHA-ITDA relationship is on the first order linear, with a coefficient of 3.3 m/cm. A similar control of mesoscale features on mixed layer depth has been recently found in a modeling study for the Southern Ocean (Bachman et al., 2017). The very weak eddy-induced ITDA in summer is again suggestive of eddies’ reduced capability in changing the upper ocean temperature, and can be attributed to the same reasons as for SSTA. The excluded eddies in the western subpolar region, as mentioned about, also exhibit significant seasonality in ITDA (not shown), as a result of the disappearance of the very deep background ITD in summer. As for SSTA, abnormal SSHA-ITDA relationship is also noticed.

Finally, eddy-induced HCTA is shown as a function of eddy amplitude in Figure 7G. A better (less-scattered) linear relationship than for SSTA and ITDA is evident, with a coefficient of 0.22 GJ/cm. This is a result of HCT being a vertical integral representing the entire upper layer’s heat content, which is partially related to surface elevation due to thermal expansion. Mesoscale eddies are 3D structures (Dong et al., 2014; Liu et al., 2020), and are therefore more capable in modulating HCT than SST and ITD which depend on additional driving forces such as air-sea heat and freshwater fluxes. Seasonal contrast for eddy-induced HCTA is not obvious (Figures 7H,I), indicating less control on HCTA by the seasonally varying surface fluxes. Despite the better linear relationship, scatters do exist, and there are also some eddies with abnormal sign of HCTA.

Eddy-Induced Variabilities

In the above we have investigated the eddy-induced anomalies in terms of their mean values and statistical relationship with the eddy amplitude. However, it remains unknown how much variability do the eddies explain on different timescales. In this section we further examine this problem by reconstructing fields of eddy-induced time series. To do so, at each time step, areas within the individual eddy boundaries are filled with the anomaly values induced by the respective eddies, while areas not belonging to any eddies are set to zero. For each eddy, the eddy-induced anomaly value calculated above is assigned to the eddy’s body center, and from there the anomalies gradually reduce to zero along the radial direction toward the boundary. See Figure 8 for a schematic diagram of this “matching and filling” procedure. In this way, eddy boundaries always have zero anomaly, and large-scale signatures are prevented because the eddy-induced anomaly is already the difference between the interior maximum and the ambient mean. Here no eddy screening is made whatsoever, and there is no presumption on the shape of the eddy boundary and/or its symmetry. A simpler way, i.e., directly set eddy interior values to be the same as the original field, is subject to residual large-scale anomalies that the eddy anomalies superimpose on.

FIGURE 8

Figure 8. Schematic diagram showing the eddy reconstruction process. The eddy under study is indicated by the area enclosed by the thick black curve, with the black dot representing its body center. The two green circles show the equivalent circle and the ambient circle, i.e., a circle with 3 times the radius of the equivalent circle. When reconstructing, the eddy anomaly is assigned to the body center, and decreases linearly in the radial direction toward the boundary. Dashed curves indicate other eddies, either with positive (red) or negative (blue) anomalies. The shading in the background is to indicate the large-scale imprint left over after removing the climatology, and is not present in the reconstructed eddy field.

Conventionally, spatial filtering is more often employed to extract the mesoscale from the background. However, the filtering approach is subject to several difficulties and is therefore believed less suitable for this study. (A) First, filters are usually not capable of separating the scales in a clear-cut manner, instead the resultant scales are interdependent. Such spectral leakage can be large for some filters, for instance the universally employed simple boxcar running mean filter. See Zhou and Cheng (2021b) for more discussions about filters. (B) Second, since large meanders of the background currents, oceanic fronts, and mesoscale eddies are all mesoscale features, spatial filtering is not efficient in separating them (Zhang et al., 2019). To alleviate this problem, elongating the filtering window in the zonal direction, say, to 15°× 5°, could yield some improvement, but cleaner separation could only be achieved by using eddy detection approaches (Zhou and Cheng, 2021a). (C) Third, spatial filters are built on the assumption of homogeneity, which means the characteristics of perturbations on different scales must not change over space, and are therefore better used in a small area. Assassi et al. (2016), e.g., showed that the size of the filtering window and the length scales of the eddy and the background current are all important. In reality, length scales of eddies and the background current vary significantly over the large basin of the North Pacific. Using a moving window with adaptable size according to local eddy properties may be helpful to lessen this problem, but it is too complicated to implement and is subject to other problems due to e.g., the a priori determination of eddy properties. (D) Last, mesoscale fields of different variables extracted by the same filter are not necessarily physically consistent. For example, mesoscale eddies with closed SSHA contours do not necessarily have SSTAs with closed contours too. It is thus difficult to interpret the results involving different variables. In contrast, the matching and filling method used in this work is free from all these kinds of problems because it is based on individual eddy processes instead of statistical tools.

Figure 9A shows an example snapshot of the reconstructed SSHA, in comparison with the original field shown in Figure 1A. It is obvious that the reconstructed SSHA captures the mesoscale features but leaves out the larger-scale imprint. As a confirmation to this, the 2D spatial spectra in the KE region shown in Figure 9B exhibits spectral peaks at wavelengths of a few degrees in both the zonal and meridional directions. In the meridional there also exists a peak at about 10–15° length scale reflecting the spatial inhomogeneity of eddy activity, which would be a problem for conventional filtering approach. The original and reconstructed time series at a selected grid point (150°E, 35°N) in the KE region are shown in Figure 9C. It could be inferred that at this point, mesoscale eddies are responsible for many of the negative SSHA events and some of the positive events, both in terms of magnitude and timing. The time-domain power spectra in Figure 9D indicates that at this point, eddy-induced SSH variability is more energetic than the original signal in the high-frequency band of 2–10 days. The reason for the overestimate is because of the intermittency of eddies in the Eulerian sense. For a certain grid point, SSHA in the eddy-vacant periods are exactly zero, but when an eddy approaches, it plunges or rises abruptly. The same happens when an eddy leaves. Faster, smaller, stronger eddies cause more abruptness in the resultant reconstructed time series. Moreover, eddies sometimes disappear from the detection algorithm at one time step but reappear at the next, presumably due to background signals that interfere the detection (Faghmous et al., 2015). This introduces even more intermittency and abruptness. When estimating the spectra using Fourier transform, therefore, harmonics at the highest frequencies are overly energetic. Applying smoothing to the reconstructed time series beforehand would of course reduce the abruptness and hence the overestimated spectral power, but the time series in the eddy-occupied periods, not only around the arrivals/departures, would also be smoothed and thus the high-frequencies are still unrealistic. In the future, improvements to the spatial and temporal resolutions are expected to lessen this problem. We investigated a large number of geographical points and found that lower frequencies than 1/10 cpd are mainly free from the overestimate problem (not shown), hence the band of 2–10 days, which is known as the synoptic scale, is discarded in the following analysis.

FIGURE 9

Figure 9. (A) Snapshot of the reconstructed eddy-induced SSH (cm) field on January 1, 2008. (B) Normalized 2D power spectrum of reconstructed SSH in the Kuroshio Extension region (145°–180°E, 30°–40°N), averaged over the year 2008. (C) Time series of the original and reconstructed eddy-induced SSH (cm) at 150°E, 35°N. (D) Power spectral density of the original and reconstructed time series shown in (C).

We then investigate the horizontal distribution of mesoscale eddies’ contribution to the variability of the studied quantities. Bandpass filters are applied on the original and reconstructed fields to extract the variability on submonthly (10–30 days), intra-seasonal (30–90 days), intra-annual (90 days to 1 year), and interannual (1–5 year) timescales. Eddies’ contribution is then measured as the ratio between the standard deviations of eddy-induced and original fields. The results are shown in Figures 10, 11.

FIGURE 10

Figure 10. Ratio of standard deviation of eddy-induced and original variability for (A,C,E,G) SSHA, and (B,D,F,H) SSTA on (A,B) sub-monthly (10–30 days), (C,D) intraseasonal (30–90 days), (E,F) intra-annual (90 days–1 year), and (G,H) inter-annual (1–5 years) timescales. The black line indicates the climatological position of the KE jet at 35°N. Note the different color scale for SSTA. Contours show the standard deviation of the eddy-induced fields.

FIGURE 11

Figure 11. Ratio of standard deviation of eddy-induced and original variability for (A,B,D,F) ITDA, and (C,E,G) HCTA on (A) sub-monthly (10–30 days), (B,C) intraseasonal (30–90 days), (D,E) intra-annual (90 days–1 year), and (F,G) inter-annual (1–5 year) timescales. The black line indicates the climatological position of the KE jet at 35°N. Contours show the standard deviation of the eddy-induced fields.

For SSH, it is evident that eddy-induced variability explains much or nearly all of the submonthly SSH variability in the regions with vigorous eddy activity, i.e., along the Kuroshio and its extension, the STCC, and along the Aleutian Islands (Figure 10A). Some high percentages could also be found along the North American coast, where active Kelvin wave-Rossby wave conversion takes place (Clarke, 1983; Clarke and Shi, 1991; Fu and Qiu, 2002). In the tropical eastern Pacific, tropical instability waves are responsible for a large fraction of the submonthly SSH variability too. In the subtropical and subarctic ocean interior, where eddy activity is less pronounced, the explained standard deviation is close to zero. The pattern of explained SSH variability resembles the pattern of absolute mean eddy amplitude (Figure 1D), except that in the less pronounced eddy active regions than the KE (e.g., the STCC), the weaker eddies do explain a large fraction of total variability. Eddy contribution then decays toward sub-seasonal and longer scales except for some spots along the Kuroshio and in the Gulf of Alaska (Figure 10, left column), respectively, associated with the large meanders of the former (Sugimoto and Hanawa, 2012; Usui et al., 2013) and the quasi-permanent eddies in the latter (Di Lorenzo et al., 2005; Saito et al., 2014).

In terms of SST, eddies’ contribution on submonthly scales is the largest around Japan, especially along the Oyashio Extension Front oriented SW-NE to the west of 170°E where eddies account for as much as 40–60% of the total standard deviation (Figure 10B). The Oyashio Extension Front region is characterized with the strongest SST gradient, and the most vigorous short-term SST variability in the North Pacific (Kida et al., 2015; Zhou et al., 2015; Zhou and Cheng, 2021a). High eddy contribution is also found in the Kuroshio and Oyashio Confluence region, in the northern Japan Sea, in southeast North Pacific, and along the North American coast. Comparing the eddy contribution and the absolute mean eddy-induced SSTA (Figure 4A), it could be inferred that the eddy poleward heat transport plays a more important role in the submonthly SST variability along the Oyashio Extension front than it does in the KE and Kuroshio-Oyashio Confluence region. It is thus indicative that the Oyashio Extension front, which is a quasi-permanent feature, has a significant eddy signature in the high frequency. With increasing time scale, however, eddy contribution in the Oyashio Extension decays monotonically (Figure 10, right column). The KE along 35°N, on the other hand, exhibits relatively high percentage of eddy-explained SST variability on almost every timescale up to interannual, probably due to the temperature advection by the pinched-off eddies from the highly vibrant jet axis. The absolute value of explained variability for SST is smaller than for SSH almost everywhere on corresponding timescales, indicating the greater complexity of other processes including atmospheric forcing in driving SST variability.

Eddy contribution on ITD variability (Figure 11, left column) highlights the KE, the Kuril-Kamchatka Trench, the Aleutian Island Chain, and the STCC. Particularly, in the KE region on the north flank of the climatological KE jet at 35°N, there are the largest eddy contribution (∼100%) on every time scale up to interannual, indicating the dominant role of eddy activity on ITD variability there. South of 35°N, however, eddy-induced ITDA variability is not large. We will come back to this point later after having introduced the eddy contribution to HCTA variability. Along the Kuril-Kamchatka Trench and the Aleutian Island Chain, eddies are responsible for large fractions of ITD variability on submonthly timescales, but their contribution decreases toward the longer scales. This is most likely related to the short lifespan of eddies there, as shown in Figure 1B. The fraction in the STCC is moderate on submonthly and intra-seasonal scales but becomes too patch on longer scales. Near the eastern boundary, high percentage of ITD variability is explained by eddies north of the Vancouver Island on submonthly and intra-seasonal scales but moves to the central and southern part on intra-annual and interannual scales.

Moreover, HCT variability on submonthly scales is almost entirely explained by eddies everywhere and is therefore not shown. This is reflective of the HCT’s lack of high-frequency variability because of the huge thermal inertia of the upper ocean, and the fact that eddies can induce significant heat transport as a result of their movement (Dong et al., 2014). Such high percentage is still present in the intra-seasonal band over most of the basin except for the subarctic and the subtropics where it is on the order of 40–50% (Figure 11C). On intra-annual scales, as can be inferred from Figure 11E, eddies are responsible for over 80% of the total variability in eddy-active regions on the south flank of the KE, in the Japan Sea, and near the eastern boundary. Explained percentages remain high south of the KE on inter-annual timescales while elsewhere they reduce to less than 20% expect for some patchy spots (Figure 11G). The high contribution of eddy-induced HCTA south of the KE on every timescale up to interannual indicates that variability of the subsurface thermal structure is almost entirely determined by eddies there, instead of the low frequency basic flow, i.e., the KE jet and its southern recirculation. It is known that the KE jet, an inertial jet that maximizes below the mixed layer (Mizuno and White, 1983; Kida et al., 2015), exhibits significant variability on decadal timescales (Qiu and Chen, 2005, 2010; Taguchi et al., 2007; Ceballos et al., 2009; Sugimoto and Hanawa, 2009; Kelly et al., 2010; Sugimoto et al., 2014; Qiu et al., 2017). Therefore it is expected that the eddy contribution to HCT will further decrease on decadal timescales.

At this point one might wonder why mesoscale eddies with lifespan on the order of weeks to months can explain some of the variability on longer scales. In the KE region, this is actually a reflection of the low-frequency interaction between eddies and the jet. Extensive studies have shown that the KE jet can force violent eddy generation due to baroclinic and barotropic instability (Waterman et al., 2011; Bishop, 2013; Chen et al., 2014; Fang and Yang, 2016; Yang et al., 2017, 2018; Ji et al., 2018), and the eddies can also feed back to the jet), forming a closed feedback loop on interannual-decadal timescales. Near the eastern boundary, on the other hand, the reason is that the generation of westward-propagating Rossby waves is sensitive, on interannual timescales, to ENSO-generated SSH anomalies that travel northward as coastally-trapped Kelvin waves (e.g., Kessler, 1990; Fu and Qiu, 2002).

Now we focus on the KE region and examine the relationship between eddy-induced SST, ITD, and HCT variabilities. We have noticed that the highest area of SST variability explained by eddies lies along 35°N, whereas for ITD it is on the north flank, and for HCT it is on the south flank. This means eddies can change SST on both sides, but they are more effective in changing the ITD than HCT in the north, and the opposite is true for the south. Our results are thus indicative that the mesoscale eddies are deeper (shallower) south (north) of the mean jet axis. This conclusion is consistent with the observational study of Zhang Z. et al. (2017) and the modeling study of Xu et al. (2019), which both favored concentrated subsurface eddies south of the KE. We also verified the distribution of subsurface eddies in the ECCO2 reanalysis using the criterion proposed by Assassi et al. (2016), which relies on the sign relationship between eddy-induced SSHA and density anomalies, and found similar results (not shown). Liu et al. (2020), however, gave the opposite distribution of eddy depth around the KE axis, but they only focused on surface eddies. It is thus inferable that the southern (northern) region has shallower (deeper) surface eddies and more (less) subsurface eddies.

Conclusion and Discussion

This work examined eddy-induced surface and subsurface thermodynamic anomalies, and quantitatively studied their statistical relations with eddy amplitude. This is a step forward from literature studies on eddy dynamic parameters only, or eddy-induced SST perturbation in the sense of composite mean. Further, the multi-scale assessment on eddies’ contribution to the total variability is novel and instructive. It is now clear that mesoscale eddies can induce large anomalies in upper ocean temperature on various scales. These thermodynamic anomalies can explain a large portion of the total variability. Of the entire basin, the Oyashio Extension exhibits the largest eddy-explained SST on submonthly timescales. The Kuroshio Extension region is highlighted for the largest percentage of ITD and HCT variabilities induced by eddies on every timescale up to interannual. For submonthly timescales, particularly, the percentage of eddy-induced SST variability (∼40%) is close to the previously proposed ratio of ocean internal variability relative to the total SST variability (∼50%, Smirnov et al., 2014), indicating eddies’ dominance among all kinds of ocean internal processes at this scale. Future research on upper ocean thermodynamic variability and air-sea interaction should better acknowledge the role played by ocean mesoscale eddies.

We also used altimeter SSH observations from the AVISO project (Ducet et al., 2000), and blended satellite and buoy SST observations from the NOAA OI-SST dataset (Reynolds et al., 2007), both having the same resolution as the ECCO2, to perform the same studies on the same time range. The results are similar to the ECCO2 results presented above and are therefore omitted. Furthermore, it must be noted that the 1/4°grid size of the ECCO2 and satellite datasets can only partially resolve the heat perturbations due to mesoscale eddies. Recent studies found that the eddies redistribute heat via two major physical mechanisms: water mass trapping within the eddy body, and advective effects by the large velocities around the eddy rim (Xu et al., 2014, 2016). The trapping effect is on a larger scale that is resolvable here, whereas the advection effect needs a much higher resolution. This is thus indicative of the need of submesoscale-resolving observations, modeling results, and reanalysis for better studying eddy thermal effects.

Our results suggest that in the Kuroshio and Oyashio Extension region, mesoscale eddies contribute to about 40–50% of SST variability on monthly to intra-annual timescales. Related to this, Smirnov et al. (2014) used the statistical tool of linear inverse model and proposed that on the corresponding timescales and in the same area, SST variability induced by intrinsic oceanic processes is about 50%. Although they used different data and methodology as this work, it is still useful to compare the findings. Since mesoscale eddies are just one component of intrinsic oceanic process, that the contribution of eddies is smaller than the total oceanic intrinsic process is understandable, and also indicative of the dominant role of eddies among other oceanic processes on the intra-annual and smaller timescales. Moreover, Bishop et al. (2017) concluded, also using a statistical model, that the explained SST variability by oceanic processes in western boundary current regions increases with timescale, in contrast to our finding of monotonically decreasing eddy contribution. This suggests that as the timescale grows, other oceanic intrinsic processes than mesoscale eddies, e.g., interannual variability of western boundary currents, take over and play more important roles in air-sea interaction.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JZ performed the data analysis and wrote the initial draft. GZ proposed the main ideas, took part in the data analysis, and dramatically modified the manuscript. HL, ZL, and XC participated to the discussions and contributed to the improvement of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

GZ was supported by the National Natural Science Foundation of China (grant no. 41906001), the Natural Science Foundation of Jiangsu Province (Grant No. BK20190501), and the Fundamental Research Funds for the Central Universities (Grant No. B210202137). HL was supported by the National Natural Science Foundation of China (Grant No. 41931183). ZL was supported by the Fundamental Research Funds for the Central Universities (Grant No. B19020044) and 2020 Jiangsu Shuangchuang (Mass Innovation and Entrepreneurship) Talent Program (Grant No. 2020–30448).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Assassi, C., Morel, Y., Vandermeirsch, F., Chaigneau, A., Pegliasco, C., Morrow, R., et al. (2016). An index to distinguish surface- and subsurface-intensified vortices from surface observations. J. Phys. Oceanogr. 46, 2529–2552. doi: 10.1175/JPO-D-15-0122.1