Wind stress curl as a driving force of annual waves in the upper ocean for interpreting energetics at all latitudes

The dynamics of waves and eddies in the upper ocean plays an important role in the climate variation of tropical and subtropical regions. Previous diagnoses for annual Rossby waves in oceanic model outputs manifested zon-ally alternating signals (ZASs) in the time-averaged distributions of wind input as well as pressure-flux divergence terms in the budget equation of wave energy. This is the case when the annual mean of the wind input is estimated as the inner product of simulated velocity vector and wind stress vector in previous studies. The present study proposes a new mathematical expression for estimating the wind input that is analogous to one derived from the quasi-geostrophic potential vorticity equation. Namely, the wind input is estimated as the negative of the product of pseudo-streamfunction and wind stress curl, the latter of which is associated with the horizontal divergence of Ekman velocity. This can be interpreted as replacement of kinetic energy input with gravitational potential energy input. Pseudo-streamfunction in the present study is inverted from Ertel’s potential vorticity anomaly and is seamlessly available at all latitudes. This contrasts with the quasi-geostrophic streamfunction which is singular at the equator. The new expression enables reducing ZASs in the horizontal distributions of both wind input and pressure-flux divergence terms, without harming the qualitative advantage of energy flux vectors to indicate the group velocity of waves at all latitudes.


Introduction
Variations of surface currents in the ocean are influenced by both intrinsic dynamics and wind forcing.The input of kinetic energy through wind forcing plays a significant role in driving ocean currents and turbulent mixing within the surface mixed layer (Ferrari and Wunsch 2009).The spatial variability of winds results in Ekman pumping and suction, which affect the vertical motion of water below the surface mixed layer of the ocean.On the other hand, temporal variations in wind forcing generate both inertia-gravity waves and planetary waves within the ocean.The budget of mechanical energy in the upper layer of the global ocean was estimated in Lueck and Reid (1984).They considered the sources and sinks of energy, including wind input, heat forcing, and viscous dissipation.Their analysis revealed that the wind input is 100-500 times greater than the input of available gravitational potential energy resulting from heat forcing and precipitation/evaporation.
Previous studies put forward two mathematical expressions for estimating the energy input from wind forcing.One is the classical expression which is based on the budget of kinetic energy to calculate the wind input by taking the dot product of surface velocity vector and wind stress vector.Using this expression, Wunsch (1998;Scott and Xu (2009) estimated that the global energy input from wind to geostrophic flows at the sea surface was 0.85-1.0TW (1 TW = 10 12 watts).This may be compared with the result of Alford (2001) who estimated that the energy input from wind forcing to near-inertial waves in the upper ocean amounted to about 0.29 TW.Wang and Huang (2004) estimated the wind energy input to the Ekman layer as 3 TW.Xu et al. (2016) used the classical expression to estimate the wind input to mesoscale eddies in the ocean, revealing a significant suppression of oceanic eddies by atmospheric winds.Another approach for estimating the wind input involves quasi-geostrophic dynamics.In this expression, the wind input is calculated as the negative of the product between the quasi-geostrophic streamfunction (whose vertical variation is associated with buoyancy) and the wind stress curl (which captures the horizontal divergence of Ekman velocity).This second expression was used in some studies such as the paper by Berloff and McWilliams (1999a).They investigated the spatial and temporal variability of mid-latitude oceanic gyres generated by steady winds in 1.5-layer and 2-layer quasi-geostrophic models to study a marginal stability problem, highlighting the instability of steady western boundary currents at moderate Reynolds numbers.
As a first step, the present study focuses on climatological wind forcing on ocean currents.Many studies noted the effect of wind stress curl on annual baroclinic Rossby waves (RWs) originating from eastern boundary regions at mid-latitudes.Krauss and Wuebber (1982) suggested that wind stress curl in the eastern boundary regions of the North Atlantic Ocean effectively drives RWs.Cummins et al. (1986) investigated how RWs are influenced by the annual cycle of wind stress curl in the North Pacific Ocean.They noted clear wave crests and troughs in the anomaly field of sea surface height (SSH) in regions between 20° and 44° N, indicating that RWs originated from the eastern boundary regions.While the crests and troughs of SSH near the eastern boundary were aligned parallel to the meridional direction, a southeast-northwest inclination was exhibited further away from the boundary.As noted by Schopf et al. (1981), the alternation in direction could be attributed to variations in wave speed with latitude depending on the deformation radius change.The presence of high meridional wavenumbers contributes to the generation of meridional group velocity.In the Southern Hemisphere, Reason et al. (1987) demonstrated using a numerical experiment that annual RWs are generated by the maxima of wind stress curl in three regions at 25° S or 38° S around the eastern boundary.
We note that the estimation of wind input rates for oceanic annual waves, as demonstrated in the Indian Ocean by Li and Aiki (2020), exhibited zonally alternating signals (referred to as ZASs in the present manuscript) that can complicate the understanding of the annual mean energy budget.These ZASs occur near the eastern boundary when the period of RWs is close to that of wind forcing, which is often unavoidable in model experiments using annual wind forcing.The objective of this study is to investigate the presence of ZASs in wind input across other oceans, explore the factors contributing to ZASs, and propose a method to mitigate their effects for a more comprehensive understanding of upper ocean dynamics.This investigation will be valuable for future studies aiming to gain insight into air-sea interactions from an energy perspective.The structure of this manuscript is as follows.In Sect.2, we present the results of a model experiment utilizing a shallow-water equation system associated with the first baroclinic mode forced by climatological wind forcing in the Pacific Ocean.In Sect.3, we examine mathematical expressions for estimating wind input, modifying the seamless diagnostic scheme to facilitate the analysis of annual wave energy transfer and budget in the upper layer of the Pacific Ocean.Finally, Sect. 4 provides a summary of our findings and highlights directions for future research.

Materials and methods
We have performed a 20-year experiment in the Pacific Ocean using a shallow water equation system associated with the first baroclinic mode.Investigation for the second and greater baroclinic modes should be done in a future study following the result of the present study.A total of 2400 snapshots have been stored in 20 years with an interval of 3.04 (= 365/120) days.The model integration has reached a climatological equilibrium within 20 years, and all 120 snapshots in the last year have been used in the following investigation.The present study has used the monthly climatology of wind forcing with subtracting a time mean component, as in Li and Aiki (2020), to investigate the dynamics of annual waves.We shall show for the Pacific Ocean that ZASs appear in the energy input by wind forcing.

Model setup
The model domain of the present study is from 50° S to 60° N and from 120° E to 70° W with realistic coastlines.The model equations are discretized in spherical coordinates using Arakawa's C-grid with a halfdegree resolution in both the zonal and meridional directions.Basic equations in spherical coordinates are shown in Appendix A. All equations in the main body of the present manuscript are written in Cartesian coordinates for simplicity, where the variables x and y denote the zonal and meridional directions, respectively, with values increasing eastward and northward.We have simulated the time evolution of zonal velocity u , meridional velocity v , and layer thickness h using a shallow-water equation system to read, where the symbol f is the Coriolis parameter.The sym- bols F u and F v in Eqs.(1a) and (1b), respectively, are eddy viscosity terms as represented using the Smagorinsky et al. (1965)  (1a) Gravity wave speed c in Eqs.(1a) and (1b) has been esti- mated to be 2.71 m/s for the first baroclinic mode using the vertical profiles of annual mean salinity and temperature from the World Ocean Atlas (WOA).In Eqs.(1a)-(1c), the symbol H 1 = 1107.5m is the reference thickness associated with the first baroclinic mode, estimated as based on the vertical mode decomposition.The non-dimensional parameter, α = 0.39, has been used to represent the wind-induced momentum input to the mixed-layer for the first baroclinic mode.See Appendix A in Li and Aiki (2020) for the details of the vertical mode decomposition.The climatological annual mean of the mixed-layer depth, H mix = 45.5 m, has been estimated based on a temperature difference of 0.2 °C using Argo data in the tropical Pacific Ocean, and the symbol H bottom = 4100 m is the reference depth from the sea surface to the bottom of the tropical Pacific Ocean.
The zonal phase velocity of equatorial waves on an equatorial β-plane may be written as where ω is wave frequency, k is zonal wavenumber, n is meridional mode number, and β = ∂f ∂y at the equator (Matsuno 1966).Substitution of n = −1 to Eq. (1d) yields ω/k = c representing equatorial Kelvin waves (KWs).Substitution of a set of n = 1 and ω = k = 0 to Eq. (1d) yields ω/k = −c/3 representing long Rossby waves in the first meridional mode.This corresponds to the observational result of oceanic RWs (Chelton and Schlax 1996).Let L = 15, 000 km be the zonal distance of the Pacific Ocean at the equator from 145° E to 80° W. It takes L c and 3L c for equatorial KWs and long equatorial RWs, respectively, to travel eastward or westward in the Pacific Ocean basin.The sum of these periods, L c + 3L c , is referred to as the basin mode period and is estimated to be 0.70 years for the first baroclinic mode.

Basic results
We explain the seasonal sequence of simulated quantities in the last year of the climatological experiment.Figure 1 shows snapshots in the middle of January, April, July, and October.The signals of zonal velocity (color shading in the left panels of Fig. 1) deflect toward lower latitudes as RWs propagate westward.The zonal wind stress (contours in the left panels of Fig. 1) is much stronger in January and July than in April and October, indicating the significant influence of the monsoon in Southeast Asia and Northeast Australia.In January, the wind stress is westward at 0°-20° N and 40°-20° S, and eastward at 20°-50° N and 20° S-0°, while the wind direction is reversed in July. (1d) Geopotential is written as, whose signals are more widely spread in the meridional direction of the Pacific Ocean (as shown by color shading in the right panels of Fig. 1) than that of zonal velocity.This indicates that gravitational potential energy is greater than kinetic energy in mid-and high-latitude regions, that is consistent with the quasi-geostrophic dynamics.The ZASs of geopotential in the eastern boundary regions at mid-latitude are generated by wind stress curl corresponding to the simulation result of Cummins et al. (1986).The variations in the group velocity depending on the deformation radius lead to the meridional propagation of geopotential (Schopf et al. 1981).In January, a negative and strong geopotential anomaly area is generated in low-latitude regions, as indicated by the red box in Fig. 1b.The negative signals propagate to the western boundary regions in July.
Potential vorticity anomaly may be written as Contours in the right panels of Fig. 1 represent pseudostreamfunction which is derived by inverting q using the following equation, where ∇ 2 is the horizontal Laplacian operator.This equa- tion has been solved using the Successive Over-Relaxation (SOR) method (Perrin 1997) with a coastal boundary condition of ϕ = 0 .Note that pseudo-streamfunction ϕ is available at all latitudes, in contrast to quasi-geostrophic streamfunction which may be written as p f .The signals of pseudo-streamfunction in Fig. 1 are positively correlated with that of geopotential in the Northern Hemisphere and negatively correlated in the Southern Hemisphere.The significant anomalies of pseudo-streamfunction are concentrated in the western equatorial regions.

(1e)
Figure 2 shows a set of Hovmöller diagrams in both mid-latitude and equatorial regions.Zonal velocity (color shading in the left panels of Fig. 2) exhibit clear westward propagating signals at both 35° N and the equator, indicating the influence of RWs.While zonal velocity at 35° N manifests a seasonal cycle, it takes about 2 years to travel from 130° to 140° W. On the other hand, zonal wind stress (contours in the left panels of Fig. 2) at 35° N as well shows a seasonal cycle.Zonal wind stress is positive from November to April, and negative from May to October.These are revisited later in the manuscript.At the equator, the regions of greater magnitude in zonal velocity appear on the western side of the basin where zonal velocity is eastward from June to November and westward from December to May.Geopotential (color shading in the right panels of Fig. 2) as well shows westward propagating signals.At 35° N, geopotential and pseudo-streamfunction (contours in the right panels of Fig. 2) both show a seasonal cycle with a similar period as zonal velocity.At the equator, the regions of greater values in pseudo-streamfunction anomaly are found in the western Pacific Ocean.These are consistent with the model results of Qiu (2003).We note that correlation between geopotential and pseudo-streamfunction is not found in equatorial regions where geostrophy is not robust.Pseudo-streamfunction at the equatorial may capture the signals of mixed Rossby-gravity (Yanai) waves associated with the semi-annual component of wind forcing, to be explored in a future study following the spectral analysis approach of Song and Aiki (2021).

Results and discussion
We investigate the budget of wave energy in this section.All results in this section have been derived by taking a time average of 120 snapshots (except for those otherwise noted) through the last year of model integration.

Classical expression
Manipulation of shallow water Eqs.(1a)-(1c) yields an equation for the budget of wave energy to read, where the second line is pressure-flux divergence, the third line is wind input rate, and the last line is dissipation (3a) ∂ ∂t rate (as well as buoyancy production rate).In this classical expression, the wind input is estimated as the inner product of simulated velocity vector and wind stress vector, as explained in Sect. 1 of this manuscript.The third line of Eq. (3a) can be interpreted as the input of kinetic energy through wind forcing.This expression of the wind input has often been used in previous studies (Wunsch 1998;Scott and Xu 2009;Ferrari and Wunsch 2009;Zhai et al. 2012;Renault et al. 2016).See Aiki et al. (2016) for the separate budget equations of kinetic energy and gravitational potential energy to validate expressions for the forcing for each.
We have plotted the time-averaged distribution of some terms in Eq. (3a) to investigate the energy budget of annual waves in the Pacific Ocean.Despite representing the annual mean, Fig. 3a, b manifests sign-indefinite signals in the wind input and the pressure-flux divergence, respectively, that we refer to as ZASs.The patterns of ZASs in both hemispheres converge closer to the equator as moving from east to west.The dissipation rate in Fig. 3c exhibits greater values in the western equatorial region.We note that the classical expression of the wind input yields distinct ZASs in regions near the eastern boundary around 30°-40°N and 40°-30° S in the Pacific Ocean.The results correspond to McCreary et al. (1987) that the annual RWs vanishes around 40°N in the North Pacific Ocean and the observations (Qiu et al. 1997).In the Northern (Southern) Hemisphere, the isolines of these signals are in the northeast-southwest (southeast-northwest) direction that is symmetric about the equator.In the Indian Ocean as well, ZASs appear in the wind input when estimated using the classical expression (Li and Aiki 2020).Li and Aiki (2020) did not conduct an in-depth investigation into the origin of ZASs, apart from their general speculation that ZASs in the Indian Ocean were believed to result from a phase relationship between ocean flows and wind forcing.The present study focuses on these signals that may disturb the understanding of the energy budget of annual waves.
To investigate the origin of these signals, we plot a Hovmöller diagram in Fig. 4a at 35° N between 140° and 130° W, representing a typical region for ZASs.The time series of the classical wind input rate shows westward shifting with a seasonal cycle.The pattern of the wind input originates from a phase relationship between ocean flows and wind forcing, as illustrated in Fig. 2a

The modified (phase-resistant) expression of the present study
As explained in Sect. 1 of this manuscript concerning the quasi-geostrophic dynamics, wind input may also be estimated as the negative of the product of streamfunction and wind stress curl (Scott and Straub 1998;Berloff andMcWilliams 1999a, 1999b).Streamfunction in the quasigeostrophic dynamics is associated with (the vertical integral of) temperature anomaly (through the hydrostatic balance).Wind stress curl is associated with the horizontal divergence of Ekman velocity (and thus the vertical convergence of velocity through the incompressibility condition).The product of streamfunction and wind stress curl is closely linked to the budget of gravitational potential energy, understanding vertical integration by parts.This should not be confused with whether an injected energy is efficiently stored in a kinetic energy reservoir or a gravitational potential energy reservoir according to relationship between wavelength and the deformation radius.
With an intent to introduce the idea of the quasi-geostrophic expression mentioned above, we propose to modify the wind input term as in the third line of the following equation to read, wind is wind stress curl.Pseudo-streamfunction ϕ has been defined in Eq. ( 2b) of the present manuscript, and ϕ x = ∂ x ϕ and ϕ y = ∂ y ϕ .Pressure-flux divergence term at the second line of Eq. (3b) has also been modified from that in Eq. (3a) to maintain mathematical consistency.The third line of Eq. (3b) represents the modified version of the wind input term, following the expression in the quasi-geostrophic literature.In contrast to quasi-geostrophic streamfunction which may be written as p f , pseudostreamfunction ϕ is available at all latitudes that is a cor- nerstone in the present study.
Figure 5a, b shows the modified wind input and the modified pressure-flux divergence of annual wave energy estimated by the modified expression in Eq. (3b).The modified wind input (Fig. 5a) becomes more concentrated with small anomaly areas at the equatorial region.The modified pressure-flux divergence (Fig. 5b) exhibits a similar pattern to the modified wind input.It can be said that the new (modified) expression provides a phase-resistant feature for estimating both the wind input and the pressure-flux divergence.
Terms on the third line of Eq. (3b) may be better explained if we decompose the horizontal component of velocity as, where superscripts "ekm", "geo", and "igw" stand for velocity components associated with Ekman flow, geostrophic flow, and inertia-gravity waves, respectively.
One way to define the Ekman velocity is to use the reference thickness to read, which is directed perpendicular to wind forcing vector.The result is that there is no wind input to the (3b) ∂ ∂t Ekman component of oceanic velocity in a single-layer framework, which is consistent with a classical framework with vertical profiles wherein wind input to Ekman spiral flow at the sea surface (Wang and Huang 2004) is canceled in sharp by the vertical integral of dissipation associated with turbulent viscosity, as detailed in Appendix B of the present manuscript.
Geostrophic velocity is defined in the present study using pseudo-streamfunction to read, (5b) which is available at all latitudes.Geostrophic velocity satisfies the following, where the second line has been derived using Eq. ( 6a), the third line has been derived using ϕ ≈ p f , the last line has been derived using Eq. ( 5a).This makes it possible for the sum of additional terms in the modified expression to offset the wind input rate in the classical expression in such a way as to reduce ZASs in the annual mean estimation. (6b) The last line of Eq. ( 6b) is analogous to the last term of Eq. (3b), allowing for to be interpreted as a variant of thickness forcing F h which leads to the source/sink of gravitational potential energy.
On the other hand, velocity associated with inertiagravity waves, coastal KWs, and equatorial KWs is characterized by no perturbation of potential vorticity as well as pseudo-streamfunction, where q has been defined in Eq. ( 2a).Substitution of Eq. ( 7a) to the third line in Eq. (3b) yields, which retrieves the classical expression of wind input.
To summarize, it is only for geostrophic velocity that the modified expression of wind input can be considered as representing the input of gravitational potential energy while canceling out the input of kinetic energy.To investigate the elements of reducing ZASs, we plot Fig. 4b, the Hovmöller diagram of which is the sum of additional terms of the modified wind input expression with respect to the classical expression.The sum of additional terms in the modified expression manifests evident consistency with the westward shifting pattern of the classical expression, except that they have opposite signs.
To investigate the details of the ZAS reduction in the Pacific Ocean, we have estimated the magnitude of the wind input associated with the classical expression using as shown in Fig. 6a where the areas of greater orders in magnitude near 30°-40°N and near the coastlines of the North American continent, corresponding to regions where ZASs are significant.We have also estimated the magnitude of the modified expression using as shown in Fig. 6b.Both terms (9a) and (9b) are in a unit of W/m 2 .Signals in Fig. 6b are smaller by one to two orders of magnitude than that in Fig. 6a in the mid-and low-latitude regions of the Pacific Ocean, especially in the regions of ZASs.To summarize the phase-resistant feature is given by cancelation within each of the second and last terms on the left-hand side of Eq. ( 6b) as shown by its first equality.Such remedy is not required for the product of streamfunction and wind stress curl.

Validation of the modified energy flux
We have plotted energy fluxes to help better understand the annual mean budget of wave energy in the Pacific Ocean.There are three schemes for illustrating the transfer routes of wave energy based on the group velocity, as explained in Appendix C of the present manuscript.Below we show the distributions of the energy fluxes based on the three schemes, as their performance in the Pacific Ocean has not been compared in the previous studies.
Figure 7 shows energy flux vectors associated with the pressure-flux scheme (Cummins and Oey 1997), and the Orlanski-Sheldon scheme (Orlanski and Sheldon 1993, hereafter OS93 scheme) following expressions (16a) and (16b), respectively, in Appendix C of the present manuscript.The energy fluxes of both schemes are compatible with the wave energy budget Eq.(3a). Figure 7a for the pressure-flux scheme manifests westward flux at the equator and eastward flux in regions between 8° S and 8° N.This scheme assumes only inertia-gravity waves but not mid-latitude RWs.On the other hand, Fig. 7b for the OS93 scheme is able to display the westward transfer of wave energy by off-equatorial RWs.The distinguishing feature of the OS93 scheme is the lack of values in the vicinity of the equator and coastlines, associated with equatorial and coastal KWs, respectively.Figure 7 illustrates the lack of validity in the pressure-flux and OS93 schemes for determining the group velocity of annual waves in the Pacific Ocean.This result is in line with similar diagnoses for annual waves in the Indian Ocean (Li and Aiki 2020) and the Atlantic Ocean (Song and Aiki 2020).
Figure 8a shows the annual mean distribution of the energy flux as estimated by the Aiki-Greatbatch-Claus scheme (Aiki et al. 2017, hereafter AGC17 scheme), following expression (16c) in Appendix C of the present manuscript.The AGC17 energy flux is compatible with the classical wind input in the wave energy budget Eq.(3a).The AGC17 energy flux (Fig. 8a) can reveal both the westward flux of energy associated with mid-latitude RWs and the eastward flux of energy associated with equatorial KWs.In the Pacific Ocean, equatorial KWs traverse toward the eastern boundary, undergo reflection, and subsequently diverge into two packets of westward RWs.These wave packets then coalesce into a single packet of westward propagating waves at 160° W as indicated by Toyoda et al. (2021Toyoda et al. ( , 2023)).We need to modify the AGC17 energy flux in such a way as to be compatible with the modified wind input in the wave energy budget Eq.(3b).This is referred to as the modified AGC17 scheme (hereafter, AGC17M scheme), following expression (16d) in Appendix C of the present manuscript.Figure 8b shows the annual mean distribution of the AGC17M energy flux.The striking resemblance between Fig. 8a, b implies that modifying both the wind input and the energy flux vectors has minimum impact on the qualitative superiority of the AGC17 scheme in terms of seamless connection between equatorial and off-equatorial regions.Figure 9 shows the time evolution of zonal flux and wind input rate during the last year of the model experiment at both 8° N and the equator.This figure is designed to facilitate a more comprehensive evaluation of the disparities between the AGC17 energy flux (left panels of Fig. 9) and the AGC17M energy flux (right panels of Fig. 9).Zonal flux (color shading) in all panels of Fig. 9 exhibits the evident signals of westward transfers associated with RWs. Figure 9c, d exhibits the weak signals of eastward transfers in the eastern domain associated with equatorial KWs.There is almost no difference between the zonal flux of wave energy (color shading in Fig. 9) between the two schemes.This confirms that the new expression in the present study continues the superiority of the AGC17 scheme in the Pacific Ocean.At both 8° N and the equator, the wind input rate of the classical expression (contours in Fig. 9a, c) shows a much stronger value generally than the modified wind input rate (contours in Fig. 9b, d).At 8° N, significant values in the modified wind input are distributed in the central and eastern basins.At the equator, the modified wind input exhibits significant negative values in the western boundary regions, indicating that layer thickness anomalies are dumped by the horizontal divergence of Ekman velocity.

Conclusions
The present study has focused on the annual mean energy budget in the upper layer as a basis for understanding the dynamics of large-scale waves.In the classical expression of the energy budget, wind input is estimated as the inner product of horizontal velocity vector and wind stress vector.Application of this expression to the energy diagnosis of annual waves yields ZASs in the annual mean wind input as shown by Li and Aiki (2020) for the Indian Ocean and the present study for the Pacific Ocean.We have investigated the reason why ZASs appear in the classical expression.Velocity associated with oceanic annual RWs contains both temporal and zonal variations  7 except for a the AGC17 scheme and b the AGC17M scheme.The AGC17 scheme assumes the group velocity of all waves at all latitudes while wind stress varies mainly in time to represent a seasonal cycle.The phase relationship between velocity and wind stress causes ZASs which disturb the understanding of the annual mean wind input.
We note that, in quasi-geostrophic dynamics, wind input is determined as the negative of the product of streamfunction and wind stress curl.Geostrophic streamfunction is associated with (the vertical integral of ) temperature anomaly (through the hydrostatic balance), while wind stress curl is associated with the horizontal divergence of wind-induced Ekman velocity (which may be rewritten as the vertical convergence of velocity using the incompressibility condition).The product of streamfunction and wind stress curl is closely linked to the budget of gravitational potential energy (after vertical integration by parts).This quasi-geostrophic expression of the wind input is singular at the equator, so the present study has modified it using pseudo-streamfunction which is obtained from inverting the anomaly of Ertel's potential vorticity.The modified expression of wind input takes the form of representing the input of gravitational potential energy with canceling out the input of kinetic energy.Because it uses pseudo-streamfunction, the modified (phase-resistant) expression of the wind input is applicable to waves at all latitudes.
To maintain the mathematical consistency between equations for the budget of wave energy, modification of the wind input affects the expression of the wave energy flux.The AGC17 energy flux is compatible with the classical wind input in the wave energy budget Eq.(3a), while the AGC17M energy flux is compatible with the modified energy input in the wave energy budget Eq.(3b).The striking resemblance between Fig. 8a, b implies that modifying both the wind input and the energy flux vectors has minimum impact on the qualitative superiority of the AGC17 scheme in representing the group velocity of waves at all latitudes.
The result of the present study may be applied to the investigation of interannual variations in such a way as to reduce the complexity associated with a seasonal cycle in oceanic flows.We plan to utilize a model with more realistic setup and higher vertical resolution in the future.The realistic model would enable us to separate geostrophic flow and Ekman flow.Whether it is still effective to reduce ZASs will be checked.Goddard and Philander (2000) pointed out that sea surface temperature (SST) and available potential energy are closely correlated, implying a connection between the wind input of wave energy and El Niño events.Toyoda et al. (2021Toyoda et al. ( , 2023) ) diagnosed an ocean reanalysis dataset in terms of the temporal evolution of wave energy associated with ENSO.They provided an energetic view for investigating ENSO which had often been understood in terms of raw quantities including SST, SSH, and thermocline depth in previous studies.In future studies, the energy perspective is expected to be used more extensively to enhance our understanding of the role of waves in tropical climate variations.

Appendix A: Equations in spherical coordinates
All equations in the main body of the present manuscript are written in a Cartesian coordinate system, for simplicity.Both the model code of shallow water equations and the diagnosis code of energy quantities have been written in a spherical coordinate system, which is explained below using and θ as longitude and latitude coordinates, respectively.
The shallow-water Eqs. ( 1a)-(1c) in Sect.2.1 may be expressed in spherical coordinates to read, where R is the radius of the Earth.Manipulation of Eqs.(10a)-(10c) yields an equation for the budget of wave energy to read, which corresponds to Eq. (3a).The pressure-flux divergence and wind-forcing terms in Eq. (11a) may be rewritten in such a way as to yield a term associated with wind stress curl to read, which corresponds to Eq. (3b) and curl tau = 1 ∂ ∂t ∂ ∂t H 1 , (16a) H 1 ρ 0 up, H 1 ρ 0 vp, of the present study yields Fig. 7a.This scheme cannot indicate the group velocity of mid-latitude RWs, to be consistent with the results of Li and Aiki (2020) for the Indian Ocean and Song and Aiki (2020) for the Atlantic Ocean.
The OS93 scheme may be used to investigate the energy flux of mid-latitude RWs, that is written as, where subscripts x and y indicate partial differentia- tion in the zonal and meridional directions, respectively.Application of this scheme to the result of the model experiment yields Fig. 7b.This scheme is singular at the equator, and can indicate neither the energy flux of equatorial waves nor that of coastal KWs.These problems were resolved in the AGC17 scheme which is seamlessly applicable to waves at all latitudes with the following expression, where ϕ is defined by Eq. (2b) in the present manuscript.Application of this scheme to the result of the model experiment yields Fig. 8a.The AGC17 scheme was used for investigating waves in the Indian Ocean (Ogata and (16b) Aiki 2019;Li and Aiki 2020, 2022, 2024;Li et al. 2021), the Atlantic Ocean (Song and Aiki 2020, 2021, 2023;Song et al. 2023a, b), the Pacific Ocean (Toyoda et al. 2021(Toyoda et al. , 2023)), and the atmosphere (Aiki et al. 2021).In Toyoda et al. (2021Toyoda et al. ( , 2023) ) and Song and Aiki (2023), the horizontal variation of gravity wave speed c was consid- ered in their calculation of the pseudo-streamfunction.
To be consistent with the modified (phase-resistant) expression of the wind input in the present study, we propose to modify the expression of the AGC17 energy flux as follows, which is referred to as the AGC17M scheme.The zonal component of the AGC17M energy flux is shown in Fig. 8b.Comparing Figs.8a, b and 9, we confirm that the modification of the energy flux maintains the qualitative superiority of the AGC17 scheme both in the horizontal distribution and time evolution through the year.

(
See figure on next page.)Fig. 1 Snapshots through the last year of the basic climatological model experiment.The snapshots are taken in the middle of a, b January, c, d April, e, f July, and g, h October.Left panels a, c, e, and g color shading represents simulated zonal velocity u (with a unit of m/s) and contours represent zonal wind stress τ x wind (with a contour interval of 0.02 N/m 2 , all positive values are shown by solid lines and all negative values are shown by dashed lines), and arrows represent wind stress vectors.Right panels b, d, f, and h color shading represents simulated geopotential anomaly p c which has been rescaled to be compared with the magnitude of velocity (with a unit of m/s) and contours represent pseudo-streamfunction ϕ inverted using Eq.(2b) (with a contour interval of 4000 m 2 /s, all positive values are shown in solid lines and all negative values are shown in dashed lines)

Fig. 2
Fig. 2 Hovmöller diagrams for the results in the last year of the climatological model experiment.The research areas are at a, b 35° N, 140°-130° W and c, d the equator, 120° E-70° W. Left panels a, c: color shading represents simulated zonal velocity u (with a unit of m/s) and contours represent zonal wind stress τ x wind (with a contour interval of a 0.004 N/m 2 , and c 0.02 N/m 2 , all positive values are shown in solid lines and all negative values are shown in dashed lines).Right panels b, d: color shading represents simulated geopotential anomaly p c which has been rescaled to be compared with the magnitude of velocity (with a unit of m/s) and contours represent pseudo-streamfunction ϕ inverted from Eq. (2b) (with a contour interval of b 50 m 2 /s, and d 4000 m 2 /s, all positive values are shown in solid lines and all negative values are shown in dashed lines) for the zonal component.Taking the time average of signals in the red box areas of Fig. 4a yields positive values.Likewise, taking the time average of signals in the green box areas of Fig. 4a yields negative values.These lead to ZASs in the annual mean field of the wind input.

Fig. 3
Fig. 3 Annual mean energy quantities as estimated by the classical expression.Color shading in each panel shows a wind input rate, b pressure-flux divergence, and c energy dissipation rate.All three terms are with a unit of W/m 2 and are based on Eq. (3a)

Fig. 4
Fig. 4 Hovmöller diagrams associated with wind input during the last year of the climatological model experiment.This analysis is at 35° N from 140° to 130° W. Color shading represents a wind input rate in the classical expression, which is written as uτ x wind + vτ y wind in Eq. (3a), and b the sum of additional terms −ϕcurl tau + ϕ y τ x wind − ϕ x τ y wind , associated with the modified wind input in Eq. (3b).Both terms are with a unit of W/m 2

Fig. 5
Fig. 5 Modified annual mean energy quantities as estimated by the modified (phase-resistant) expression in the present study.Color shading in each panel shows a modified wind input rate, and b modified pressure-flux divergence.Both terms are with a unit of W/ m 2 and are based on Eq. (3b)

Fig. 6
Fig. 6 Comparison between the two expressions of wind input rate in terms of magnitudes.a The classical wind input rate as measured by term (9a), and b the modified wind input rate as measured by term (9b).The comparison is based on the result of the climatological model experiment during the last year.Both magnitudes are with a unit of W/m 2

Fig. 7 Fig. 8
Fig. 7 The annual mean flux of wave energy as estimated by different diagnosis schemes.a The pressure-flux scheme assuming the group velocity of inertia-gravity waves and b the OS93 scheme assuming the group velocity of quasi-geostrophic waves.Color shading shows zonal energy flux (with a unit of W/m).Arrows represent energy flux vectors, with thin arrows for the energy flux in the range of 200-1000 W/m, and thick arrows for greater ranges.The result is based on the climatological model experiment during the last year

Fig. 9
Fig. 9 Hovmöller diagrams for examining the annual mean flux of wave energy.The analyses are at a, b 8° N and c, d the equator.Color shading represents zonal energy flux (with a unit of W/m) estimated by a, c the AGC17 scheme and b, d the AGC17M scheme without a time average.Contours represent the wind input rate to oceanic waves (with a contour interval of 10 -4 W/m 2 ) as estimated by a, c the classical expression and b, d the modified (phase-resistant) expression where all positive values are shown in solid contours and all negative values are shown in dashed contours.The diagrams represent the last year of the climatological model experiment