A Mechanistic Analysis of Tropical Pacific Dynamic Sea Level in GFDL-OM4 under OMIP-I and OMIP-II Forcings

The sea level over the tropical Pacific is a key indicator reflecting vertically integrated heat distribution over the ocean. Here we use the Geophysical Fluid Dynamics Laboratory OM4 (GFDL-OM4) global ocean-sea ice model forced by both the CORE and JRA55-do atmospheric states (OMIP-I and OMIP-II) to evaluate the model performance and biases compared against available observations. We find persisting mean state dynamic sea level (DSL) bias along 9◦N even with updated wind forcing in JRA55-do relative to CORE. The mean state bias is related to biases in wind stress forcing and geostrophic currents 5 in the 4◦N to 9◦N latitudinal band. The simulation forced by JRA55-do significantly reduces the bias in DSL trend over the northern tropical Pacific relative to CORE. In the CORE forcing, the anomalous westerly wind trend in the eastern tropical Pacific causes an underestimated DSL trend across the entire Pacific basin along 10◦N. The simulation forced by JRA55-do significantly reduces the bias in DSL trend over the northern tropical Pacific relative to CORE. We also identify a bias in the easterly wind trend along 20◦N in both JRA55-do and CORE, thus motivating future improvement. In JRA55-do, an accurate 10 Rossby wave initiated in the eastern tropical Pacific at seasonal time scale corrects a biased seasonal variability of the northern equatorial counter-current in the CORE simulation. Both CORE and JRA55-do generate realistic DSL variation during El Niño. We find an asymmetry in the DSL pattern on two sides of the equator is strongly related to wind stress curl that follows the sea level pressure evolution during El Niño.

extends over years . The OMIP protocol is detailed in Griffies et al. (2016) and Tsujino et al. (2020), with the use of two forcing datasets allowing us to assess robustness of simulated features and to better attribute biases. Following Tsujino et al. (2020), all simulations are spun-up by running for five cycles of the respective forcing datasets, over which time the upper ocean reaches a quasi-equilibrium state. We then present our analysis for the sixth forcing cycle.

Observationally based datasets 70
Total sea level changes can be derived from satellite altimetry. The Copernicus Marine Environment Monitoring Service (CMEMS) provides a gridded sea level dataset by combining multi-mission altimeter satellite data since 1993 (https://resources. marine.copernicus.eu/). The monthly sea surface height (SSH) anomaly from the geoid is used in this study. Following the definition of DSL according to OMIP Gregory et al., 2019), the area-weighted global mean is removed from the observational data and model data to make them comparable. 75 For steric sea level anomaly and the depth-integrated density changes, we use the EN4 dataset from the Met Office Hadley Centre, which provides quality-controlled subsurface ocean temperature and salinity profiles based on objective analyses (Good et al., 2013) (https://www.metoffice.gov.uk/hadobs/en4/). The gridded temperature and salinity profiles are used to derive density changes at all available grid cells. We also calculate the density changes solely caused by the thermal expansion (Roquet et al., 2015). 80 To investigate the possible bias of the wind forcing, we use the Wave and Anemometer-based Sea Surface Wind (WASwind) data (Tokinaga and Xie, 2011). The WASwind provides a global coverage of zonal and meridional wind stress based on wind observation from the International Comprehensive Ocean-Atmosphere Dataset at the monthly frequency with 4 • by 4 • grid resolution from 1950-2009. Through height-correction for the anemometer-measured winds, the WASwind is not subjected to the spurious upward trend due to increases of anemometer height in the ship-based measurement. The dataset also incorporates 85 the estimated winds from wind wave height. We find this dataset to be suitable for our analysis of the tropical Pacific, and it complements the assessment from Taboada et al. (2018), who focused on upwelling patterns in the global ocean.

Analysis framework
Following Gill and Niller (1973), we separate time tendencies in sea level, η, into mass and steric contributions (see equation (2) in (Yin et al., 2010)) ∂η(x, y, t) ∂t = 1 ∂ρ(x, y, z, t) ∂t dz. (1) In this equation, g is the gravitational acceleration, p a is the pressure at the ocean surface due to atmospheric mass loading (which is zero in this study when discussing DSL variation in the model), and p b is the pressure at the ocean bottom. The first term on the right-hand side measures local sea level change due to changes in mass within a seawater column. In the second right-hand side term, ρ is the in-situ seawater density, with density changes integrated over the depth of a water column (from ocean bottom at z = −H to surface at z = η and as normalized by the reference density ρ 0 = 1025 kg m −3 ), yielding sea 95 level changes from density (steric) effects. The minus sign on the steric term arises since decreases in density, as from ocean warming, lead to increases in sea level. As noted earlier, the global mean sea level time series is subtracted from all regional sea levels from observations and models to allow for direct comparisons of the resulting DSL.
Based on the 1993-2017 observational data of regional mean sea level from CMEMS and steric sea level from EN4, we find that the sea level change over this period is dominated by steric effects in the tropical Pacific (defined as the region within the 100 20 • S-20 • N zonal band inside the Pacific basin). The regional averaged steric signal accounts for more than 75% of the variance in the total sea level change over the eastern Pacific (180 • -eastern boundary) and accounts for more than 85% of the variance in the western Pacific (120 • E -180 • ) [ figure 1]. Particularly, the steric signal within the upper 400 meters accounts for more than 95% of the variance in the total steric signal in both the eastern and western tropical Pacific. This result shows the central role of density changes in the upper 400 meters, which relates to the thermocline depth changes, in accounting for patterns of 105 sea level variability. The residual of the sea level variance is related to the mass component (equation (1)). These results are consistent with earlier analysis from CORE simulations documented by Griffies et al. (2014).
Density changes in the tropical Pacific are dominated by temperature changes in the upper 400 meters [ figure 1]. For the upper 400 meters, such changes can arise from surface heat fluxes as well as lateral and vertical ocean heat transport. The surface heat flux provides a thermodynamical forcing that directly causes a local thermosteric sea level change (sea level 110 changes due to thermal expansion) through diabatic heating. Surface wind forcing, on the other hand, induces a dynamical effect related to Ekman transport that causes light surface waters to diverge or converge, which in turn modifies sea level.
Surface wind stress curl causes a Sverdrup transport associated with both Ekman transport near the surface and geostrophic transport below. The internal wave propogation, like Rossby wave and Kelvin wave, at the seasonal and interannual time scale also have large effect on the sea level changes in the tropical Pacific. All these effects can have very different contributions to 115 sea level at different time scales and different spatial scales. In this study, we aim to characterize mechanisms that cause sea level variations and trends, and determine reasons for simulation biases.

Time Mean Dynamic Sea Level
In this section we focus on the time mean DSL patterns, with the time mean computed over the years 1993-2007 that are common to both JRA55-do (OMIP-II), CORE (OMIP-I), and observations [ figure 2]. This bias pattern is consistent with 120 Griffies et al. (2014) and Tsujino et al. (2020). The only processing step difference with the previous studies is the Gaussian smoothing of the observational data. There are two reasons why smoothing is not applicable for this study. One is the larger Rossby radius of deformation in the tropical region compared to the mid-latitudes. Therefore, the eddies over the tropical in the observational data. The following analyses will concentrate on the tropical Pacific region (defined as the region within the 20 • S-20 • N zonal band inside the Pacific basin).  Because the bias shows largely zonal structures, we investigate the zonal mean DSL and wind stress in the tropical Pacific. Figure 4a shows that the DSL in the simulation is generally lower than the observation by about 0.02 meter except near 10 • N.
The meridional gradient of the zonal mean DSL, which is highly correlated with the narrow zonal currents in the tropical Pacific, is comparable with the observation in the Southern Hemisphere. The Northern Hemisphere, on the other hand, shows underestimated DSL gradient in two zonal bands, 9 • N-20 • N and 4 • N-9 • N due to the lack of DSL trough around 9 • N and DSL 140 ridge at 4 • N. Between 9 • N-20 • N, the down-gradient DSL toward the equator is associated with the strength of north equatorial current (NEC). Between 4 • N-9 • N, the up-gradient DSL toward the equator is associated with the strength of north equatorial counter-current (NECC). Therefore, the smaller or missing DSL gradient in both simulations could lead to underestimated NEC and NECC.
where EB is the eastern boundaries of the tropical Pacific at each latitude, τ is the wind stress vector, β = df /dy is the meridional derivative of the Coriolis parameter, f , and ρ 0 = 1025 kg m −3 is the reference density.

The role of zonal currents in the ocean heat budget
To analyze the impact of zonal currents on heat, we consider two boxes within the tropical Pacific and assess the mean heat as shown in figure 5, and express the heat budget as In this equation, V ∂Q ∂t dV is the heat content changes in a box, which is related to the thermosteric sea level changes, 170 and F srf is the area integrated net surface heat flux. T m represents the volume mean temperature inside a box. For Wbox, T m = T mw and for Ebox, T m = T me . We calculate the heat advective transport term following the analysis of Lee et al. (2004) and Ray et al. (2018). The residual term accounts for transport from small-scale processes including vertical mixing, sub-grid scale processes, with these processes not of concern here and are small compared to the role of currents. is consistent with the finding of diabatic heating at the surface controlling the heat movement over the ocean (Holmes et al., 2019). Despite the large net surface heat flux in the Ebox, the heat is not stored in the eastern tropical Pacific but flushed to 180 the western tropical Pacific. We find that, with larger net heat flux over the Ebox, the heat advection across the 180 • transect is larger in CORE than JRA55-do. For the Wbox, the main heat loss occurs at the 120 • E transect which advects the heat out with the magnitude a little less than half of the heat coming from the eastern tropical Pacific. Heat loss also occurs across the other three boundaries of the Wbox but considerably smaller than across the 120 • E transect. Besides some small differences in surface heat flux and the heat advection at 180 • transect, we do not see significant differences in the mean state heat budget

Dynamic Sea Level Trend
Over the 1993-2012 period, satellite altimetry has shown a significant sea level rise which is associated to the warming trend 225 in the western tropical Pacific [ figure 8]. The trend has been extensively studied due to its significant asymmetric zonal pattern over the tropical Pacific which is related to global warming rates (Peyser et al., 2016;Merrifield, 2011).

Ekman layer response
To help reveal mechanisms for the underestimation of the DSL trend, we investigate the wind stress forcing in CORE and JRA55-do and compare with the WASwind observational product. We calculate the wind stress trend bias by subtracting the  same zonal band but is roughly three times smaller than the CORE simulation [ figure 12]. The smaller wind stress curl trend bias with JRA55-do is mainly due to the missing westerly trend bias in the 10 • S to 10 • N region.
The westerly trend bias in CORE forcing data is a result of the multiplicative factor (R s ) applied to the vector winds in 260 the reanalysis data by Large and Yeager (2009). The multiplicative factor is determined by the five-year mean ratio between reanalysis data and observational data. Since the factor is designed to make the mean state of wind amplitude in reanalysis data  We conclude from this analysis that the positive wind stress curl trend bias west of 150 • W, found in both CORE and JRA55do, is mainly due to the easterly wind stress trend bias. Additionally, the positive wind stress curl trend bias east of 150 • W in 270 the CORE forcing is due to the combined effect of easterly and westerly trend biases. The positive wind stress curl trend bias in the zonal band creates an artificial Ekman suction that causes the DSL trend in the simulations to be biased low.

Barotropic response
We now examine how the barotropic geostrophic response affects the DSL trend.

Rossby wave propagation
The seasonal variation of DSL in the 2 • N to 10 • N zonal band is a result of the Rossby wave propagation that is generated in the eastern and central tropical Pacific by the local wind stress curl (McPhaden et al., 1988). Figure 15 shows the spatial pattern of wind stress curl anomaly and the DSL anomaly in the eastern tropical Pacific between 2 • N to 10 • N.
The Hovmöller diagram of the meridional mean from 2 • N to 10 • N shows the propagation of this seasonal DSL signal and   The missed timing in CORE is related to the weaker Ekman pumping/suction east of 90 • W throughout the year, thus causing significant DSL biases when comparing to observations. The timing of the JRA55-do simulation is relatively close to observations but with underestimated amplitudes in both Ekman pumping and DSL. This analysis shows the importance of 320 resolving dynamics near the ocean boundary as they can strongly affect the basin-scale DSL variation on seasonal time scales.
Outside the initiation region, both forcing data show weaker Ekman suction during JJA and weaker Ekman pumping during DJF between 150 • W and 180 • when compared to observations [ figure 16b,d]. The weaker forcing does not affect the propagation of the DSL signal, but it causes the DSL signal in both simulations to be biased low west of 150 • W due to the lack of continuous external forcing.

325
The vertical Ekman-induced velocity in the Ekman layer is given by In the tropics, a vertical velocity in the Ekman layer can be generated with zero wind stress curl but strong zonal wind due to the large β effect. Therefore, the weaker Ekman suction during JJA and weaker pumping during DJF could be related to the

335
The other interesting fact is the compensating effect between zonal wind stress and wind stress curl on the Ekman mechanism in this zonal band. In other words, the first two terms on the right hand side of equation 4 which represents the wind stress curl induced Ekman effect compensates with the last term which represents the zonal wind stress induced Ekman effect. Figure   17 quantitatively visualizes this compensating effect perfectly, which shows the zonal wind related Ekman contribution and wind stress curl related Ekman contribution. During JJA, the westerly wind dominates the zonal band due to cross-equatorial    (Vinogradov et al., 2008).

Zonal currents
Due to the significant seasonal fluctuations in DSL within the tropical Pacific, we use the same box budget as in Figure 5 to study the corresponding zonal current variations between Wbox and Ebox. The zonal current anomaly is strongly affected by the meridional DSL gradient. In the 2 • N to 10 • N zonal band, we find one of the largest DSL variance. The corresponding 365 meridional DSL gradient changes, in response to the Rossby wave propagation, can have a large impact on the seasonal changes of the NECC. Figure 19a,b shows the seasonal amplitude and phase defined by the largest values of the monthly climatology and the corresponding month, respectively. The SEC and the EUC show comparable and dominating roles on seasonal heat advection between Ebox and Wbox. The seasonal amplitude of heat advection in the CORE simulation is larger than with JRA55-do for all currents, which is consistent with larger wind stress forcing and DSL amplitude in CORE.

370
All seasonal phases agreed between the two simulations except for the NECC and SEC. The SEC difference between the two simulations is mainly due to the small difference in the sub-seasonal signal that causes the difference in phase [ figure 19c].
For the NECC, the difference in phase can also be seen in the volume transport, which peaks in June for CORE and November for JRA55-do [ figure 19b,c]. The simulation forced by JRA55-do shows a better agreement with the observed NECC which peaks in December (Johnson et al., 2002).

375
The zonal current analysis confirms the importance of an accurate simulation of the DSL at the seasonal time scale. This analysis also demonstrates the crucial role of the surface wind stress timing near the eastern tropical Pacific, which influences  ]. This current reversal, however, does not exist in the observations. In the simulations, the current reversals 380 results from the underestimated mean state of the currents that leads to the incorrect seasonal heat and mass transport direction.
However, due to the small contribution in volume and heat budgets from these two currents, the model simulation is not greatly affected by this reversal at the seasonal time scale. ]. To determine the composite, we remove the mean, trend, and seasonal signals in all of the time series.

Oscillator theories
We first focus on the zonal mean of the DSL variation in the tropical Pacific region during El Niño. A clear recharge-discharge oscillator affecting DSL can be seen evolving during El Niño [figure 21a,c,e] (Jin, 1997). According to the oscillator theory, the warm water volume continues to be charged into the equatorial region before it reaches the peak of an El Niño. It then 400 discharges to higher latitude after the peak as a result of the Sverdrup transport response to wind forcing over the tropical Pacific. The SST lags the DSL positive anomaly before the peak and the negative anomaly after the peak of El Niño. This lag means that the subsurface processes build-up the warm water volume before the SST starts to warm up during the charging stage. During the discharging stage, the subsurface processes release the warm water volume before the SST starts to cool down.

405
Consistent with the theory, DSL reaches the peak at the end of Year 1 for both simulations during the charging stage which increase the warm water volume near the equator and quickly changes to a negative anomaly while the higher latitude changes from a negative to positive anomaly, which indicates the start of the discharging stage [ figure 21a,c]. The lag of the SST is also clear in both simulations. The CORE simulation shows a larger DSL bias which is consistent with the bias shown at other time scales [ figure 21b,d]. The recharge-discharge pattern can also be seen during the longer period composite  [ figure   410 21e].
We separate the El Niño period into four different stages [ figure 22]. The first is the initiation stage calculated as the mean of seven months to twelve months prior to the maximum ONI (February, Year1 to July, Year1). In this stage we see the positive DSL anomaly initiated in the central tropical Pacific along the equator. The second stage is the mean of six months prior to the maximum ONI value (August, Year1 to January, Year2). The positive DSL anomaly propagates eastward due to the 415 propagation of an equatorial downwelling Kelvin wave while the negative DSL anomaly strengthens in the western tropical Pacific in both hemispheres (Zebiak and Cane, 1987). Based on the ENSO oscillator theories, the negative DSL anomaly in 32 https://doi.org/10.5194/gmd-2020-374 Preprint. Discussion started: 30 November 2020 c Author(s) 2020. CC BY 4.0 License. the western tropical Pacific is a result of two phenomena. One is the upwelling Rossby wave propagating westward on two sides of the equator reaching the western boundary based on the delayed oscillator theory (Zebiak and Cane, 1987). The other is related to the shoaling of the off-equator thermocline based on the western tropical Pacific oscillator theory (Weisberg and 420 Wang, 1997). The latter also leads to the increase of sea level pressure due to decreasing SST and the associated easterly wind at the equator in the third stage.
The third stage is the mean of six months following the maximum ONI value (January, Year2 to June, Year2). The positive anomaly in the eastern Pacific starts to dissipate through a coastal Kelvin wave moving heat toward the poles (Johnson and O'Brien, 1990). The reflected downwelling Rossby wave on two sides of the equator, on the other hand, tries to move the heat 425 back to the warm pool in the western Pacific (Picaut et al., 1997). Like in stage 2, the DSL drop in stage 3 over the western tropical Pacific is also a result of two phenomena. Based on the delay oscillator theory, the reflected upwelling Rossby wave related DSL drop in the western tropical Pacific starts showing at the equator. On the other hand, the western tropical Pacific oscillator theory emphasizes the importance of equatorial easterly winds, due to the positive sea level pressure off-equator, on generating the negative DSL anomaly at the equator. This stage is also referring to the discharging stage in the recharge-430 discharge oscillator theory where the warm water volume is discharged poleward from the equator. The final stage is calculated as the mean of seven months to twelve months following the maximum ONI value (July, Year2 to December, Year2). In the final stage, the downwelling Rossby wave on two-sides of the equator shows in the western tropical Pacific for both simulations.
In general, the model simulations show great resemblance in terms of spatial patterns and temporal changes. Both simulations demonstrate the oscillator theories well. However, asymmetry in the DSL changes on two sides of the equator seem missing 435 from the theories mentioned above during El Niño. The asymmetric pattern, like the negative DSL signal in the southwestern tropical Pacific, is related to other mechanisms dominating the DSL signal over certain regions, which reduces the symmetry along the equator.

Asymmetry in the dynamic sea level along equator
The negative DSL anomaly and the asymmetry of DSL in the southwestern tropical Pacific are not due to the lack of discharging 440 of warm water from the equator, but instead due to the drop in DSL signal that is intensified after the peak ONI ( (a,c,e,g) 1958-2007 and (b,d,f,h) 1993-2007. The black dashed box shows the region of asymmetric sea level change in the southwestern tropical Pacific at stage 3.
turning creates a negative wind stress curl in the southwestern tropical Pacific that causes Ekman suction and the associated negative DSL anomaly.
After the ONI peak, a low-pressure center is established in the south-central tropical Pacific [ figure 24e,f]. This low-pressure center, closer to the southwest tropical Pacific than in stage 2, causes a stronger negative wind stress curl than in stage 2 due to the increased pressure gradient that results in stronger Ekman suction and DSL drop. Eventually, the dropping DSL signal 455 subsides due to decreased and reversed pressure contrast related to the Southern Oscillation [ figure 24g,h]. Since the composite based on observations only includes four El Niño events, the different sea level pressure amplitudes could be readily dominated by the strongest events. Therefore, the amplitude changes between longer  and shorter periods (1993-2007) should not be seen as a long term trend but instead a sign of possible decadal variability of the ENSO/SOI strength.   figure 24], which is a process described by the western tropical Pacific oscillator theory (Weisberg and Wang, 1997 biases in the JRA55-do simulation. An asymmetry in the DSL change during an El Niño event on two sides of the equator has not been explained in the existing oscillator theories. We find the asymmetry in the DSL pattern is strongly related to wind 515 stress curl that follows the sea level pressure evolution during El Niño. An atmospheric sea level pressure gradient created by a pressure drop in the southeastern tropical Pacific and a pressure rise in Australia as well as the sea level pressure drop in the cold tongue region are the main contributors to the negative wind stress curl that drives the DSL drop. The high correlation between wind stress curl and DSL shows the dominant influence of local wind forcing on the DSL signal. The analysis demonstrates the importance of the external forcing on the off-equator DSL changes during an El Niño event. 520

Recommendations for key bias reductions
To reduce the time mean bias of DSL in future OMIP simulations over the tropical Pacific, accurate zonal wind stress shear near the ITCZ is crucial. Especially over the 4 • N to 9 • N zonal band, the biases in the time mean of the DSL can further affect the mean ocean current strength and cause an artificial seasonal ocean current reversal.
For the bias in the long-term trend of DSL, the JRA55-do has significantly improved the westerly wind bias at the eastern 525 Pacific near the equatorial region (10 • S to 10 • N) that results in the reduction of DSL bias along 10 • N. To further reduce this DSL bias along 10 • N, we find easterly trend biases exist along 20 • N in both the JRA55-do and CORE forcings that create the artificial positive wind stress curl along 10 • N, resulting in the negative sea level bias. The multiplicative factor applied to correct the wind stress is highly related to the wind stress bias. A careful evaluation of the wind stress field across different time scales after applying the factor is necessary to reduce the bias of the simulated DSL.

530
For the seasonal variation of the DSL, the dominant external forcing is region-specific. Over the deep tropics, from the equator to around 10 • N/S, ocean dynamics are the dominating factor which was mainly initiated at the eastern boundary near 90 • W due to Ekman suction/pumping induced by local wind stress curl. Therefore, capturing the timing and strength of the local wind stress curl is crucial to better simulate the seasonal variation of the DSL. The accurate DSL variation can further improve the simulated ocean current in seasonal time scale, like the improved NECC in the JRA55-do comparing to the CORE 535 forced simulation. However, the continuously underestimated DSL amplitude in the JRA55-do simulation, which is related to the underestimated wind stress curl in the same zonal band needs further improvement. Besides the local boundary wind forcing near 90 • W, the underestimated seasonal zonal wind stress variability between 150 • W and 180 • in both JRA55-do and CORE causes the underestimated DSL signal throughout the year. Due to the compensation effect of zonal wind stress and wind stress curl on the strength of Ekman suction/pumping in the tropics, a separated evaluation of both zonal wind stress and 540 wind stress curl is needed in the future to avoid the compensated bias shown in the CORE forcing data.
For interannual variability, the El Niño related DSL variation is well represented in JRA55-do and CORE for both the timing and the spatial pattern. However, we still see a amplitude bias for DSL in both CORE and JRA55-do forced simulations. The JRA55-do does improve the amplitude bias when comparing to the CORE.
The multiple time scale analyses in this study provide a comprehensive view of important variability, changes, and the 545 associated biases in the tropical Pacific. Any improvement of the biases mentioned in this study will be very helpful to further a mechanistic understanding of tropical Pacific DSL patterns using OMIP simulations, and by extension their related coupled climate models. Author contributions. CWH, JY, SMG proposed and led this evaluation study. CWH performed the analyses including processing the model outputs and observational data. RD produced the model outputs. All authors contributed to the writing and editing processes.
Grid Federation. We thank the EU Copernicus Marine Service Information providing the Sea surface height data from multiple satellite altimetry measurements. We thank the Met Office Hadley Centre archiving the temperature and salinity data. We thank Dr. Hiroki Tokinaga on processing and archiving the WASwind data.