{'fraction_non_alphanumeric': 0.024030209406110538, 'fraction_numerical': 0.013525575008582218, 'mean_word_length': 4.9574642126789366, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulation of the Pacific equatorial thermocline with an ocean general circulation model

The reproduction of equatorial thermocline characteristics by numerical models is of great importance for seasonal forecasts, and it is also a key factor in operational ocean forecasting. Due to difficulties in the numerical description of adiabatic eddy processes, the model thermocline is usually diffused and shifted in depth. In this work, we perform a numerical experimental study of the sensitivity of the Pacific equatorial thermocline to horizontal and vertical mixing parameterizations in the INMIO ocean general circulation model. It has been shown that the sharpness and location of the thermocline can be improved by a combined parameterization of eddy mixing and isopycnal diffusion, as well as by proper tuning of the background coefficient of vertical diffusion.


Introduction
The ocean thermocline is one of the key objects whose reproduction largely determines the adequacy of an ocean general circulation model. In the middle latitudes, its physics has been studied in sufficient detail and generally does not cause problems for modelers. On the other hand, the equatorial thermocline is more difficult to reproduce. At the same time, its parameters significantly affect the structure of horizontal and vertical heat fluxes and, consequently, are key to the model climatic system characteristics. This primarily applies to the intensity and duration of the El Niño -Southern Oscillation phases [1]. In [2], an analysis of the equatorial heat budget was carried out, which showed that, depending on the sharpness of the model thermocline, the westward or vertical heat transfer becomes the dominant factor in the development of El Niño phases. Knowledge of this development is of great importance both as a source of long-term predictability for seasonal forecasts, and as a factor of short-term operational ocean forecasting. It directly influences the system of equatorial currents and upwellings that affect the state of the ocean, including the mid-latitudes [3,4].
In nature, the equatorial thermocline of the Pacific Ocean represents a very sharp temperature jump. Its depth, usually defined by the position of the 20°C isotherm, gradually increases from the eastern to the western side of the ocean. A peculiar feature of the Pacific equatorial region is the dominant role of adiabatic processes in the ocean circulation [1]. Thus, even small values of model diffusion, both induced by numerical schemes and explicitly specified, can introduce significant distortions into the structure of the thermocline. Models often experience a complex of problems in this area, which is called the "diffuse thermocline". First of all, this is expressed in a decrease of the vertical temperature gradient and, accordingly, in an increase of the interval of depths at which the In addition, the model solution may manifest an excessively sharp immersion of the thermocline when moving along the equator to the west. In the meridional section, the isotherms can strongly bend upward, which is manifested in a surface cold bias and in the corresponding distortion of the heat fluxes balance. In [5], this was interpreted as the model underestimation of the intensity of tropical instability waves, which are responsible for the flattening of isotherms through heat transfer into the region of the East Pacific cold equatorial "tongue". With the INMIO model used in that work, this bias appeared despite the high spatial resolution of the model, generally sufficient for the description of the equatorial eddy dynamics.
Many works studied the formation mechanisms of the equatorial thermocline by means of idealized examples, e.g. [6,7]. With the development of numerical modeling, a number of studies have appeared that improve the characteristics of full three-dimensional numerical models in the equatorial Pacific, mainly concerning the parameterizations of vertical mixing [2,8]. The goal of the present work is to experimentally study the sensitivity of the equatorial thermocline to the settings of the horizontal and vertical mixing schemes in a numerical model with fixed advection schemes. For this, we use the nonglobal model configuration of the isolated Pacific Ocean that has realistic boundary conditions and allows us to save computer resources during a series of experiments.

Model and experiment configuration
The INMIO model [5,9] approximates the 3DPEM system of three-dimensional ocean dynamics equations with the Boussinesq and hydrostatics approximations by the finite volume method on a Btype grid in vertical z-coordinates with an explicit description of surface fluxes and a nonlinear kinematic free surface condition. When approximating horizontal differential operators, only explicit time schemes are applied, which has made it possible to efficiently parallelize the model under the control of the CMF software package [10].
The Pacific Ocean model calculation domain is confined by solid boundaries from the north in the Bering Strait and from the west on the 116ºE meridian. In the east, the periodic boundary condition is set in the Drake Strait, allowing the waters flowing into the strait to appear from a "window" in the western boundary at the same latitudes, thus mimicking the Antarctic Circumpolar Current. A regular latitude-longitude grid is used with a horizontal resolution of 0.25 degrees, which allows partial resolving of mesoscale eddies. The vertical discretization includes 49 horizons with a step from 6 m in the upper layer of the ocean to 250 m in deep waters.
One of the key ideas of the modern general circulation model development is the most possible explicit description of the eddy dynamics. In particular, this implies the deprecation of the secondorder viscosity operator ("Laplace viscosity"), which was previously used mainly to ensure numerical stability and to simulate non-resolved lateral boundary layers. The transition to such low-viscosity mode of calculations requires providing self-consistency of the model numerical schemes at the difference level. This was performed for the INMIO model in previous works [11,5]. In the present study, this allowed us to take settings of the turbulent exchange coefficients similar to the eddypermitting configuration ORCA025 [12], which has the same horizontal resolution.
In the control experiment, the horizontal diffusion coefficient for heat and salt is taken equal to 300 m 2 /s at the equator and scaled to the poles in proportion to the square root of the grid cell area. In the momentum equations, the Laplace viscosity is not used. To ensure numerical stability, the biharmonic filter is used with a coefficient of -1.5·10 11 m 4 /s at the equator, scaled in proportion to the cell area to the power of 3/2 and with the adaptive addition according to the Smagorinsky method (in the formulation [13] with coefficient C = 1).
The background coefficients of vertical viscosity and diffusion are 10 -4 and 10 -5 m 2 /s. The shear vertical mixing is parameterized according to the Munk -Anderson scheme with maximum values equal to 10 -2 and 10 -3 m 2 /s, respectively. The time step is 10 minutes, while the barotropic velocities are calculated in a nested loop with a step of 20 seconds. To approximate the advection of momentum and tracers, we use the centered difference scheme in space and the leapfrog scheme in time. On solid lateral boundaries, the conditions of no normal fluxes of mass, heat, salt, and momentum are used. All experiments were performed for 20 years of model time, during which the fluxes at the atmosphere-ocean interface were defined by the atmospheric conditions and bulk formulas of the CORE-I protocol [14]. The model starts from the state of rest and the annual mean temperature and salinity fields of the World Ocean Atlas 2009 (WOA09) climatology [15,16]. The focus of this work is on the cold tongue of the East Pacific equatorial waters, which spreads westward from the coast of South America. Therefore, in the polar regions to simulate the sea ice we use a simple thermodynamic model [17], which does not require a large amount of computation.

Experiment results
In this section, we present the results of the series of numerical experiments. Continuing to each subsequent experiment, we will change one model parameter each time. All presented model fields will be annual means for the last (20 th ) year of experiments. Figure 1 shows the equatorial distribution of the temperature averaged in the band from 1ºS to 1ºN. Plate 1(a) presents the annual mean data of WOA2009. The thermocline is pronounced over most of the considered interval of 140ºE -80ºW. The vertical temperature gradient (hereinafter, the average between the isotherms of 16 -24ºC at 160ºW) is 0.129 ºC/m. The depth of the 20°C isotherm gradually increases from 40 m in the east to 180 m in the west, while at 160°W it is 150 m.
The control experiment described above is presented in plate 1(b). It can be seen that the thermocline is substantially diffused, especially in the eastern part of the region. The vertical gradient is 0.072 ºC/m. In the course of the experiment the upper layer of the ocean warms excessively and, hence, the depth of the thermocline is much higher than that observed in nature. From 110ºW to 85ºW, the surface temperature is overestimated by 1-1.5ºС. The isotherms of the upper part of the thermocline sink more sharply than in the observational data. The depth of the thermocline (isotherm of 20ºС at 160ºW) is 207 m.
The configuration of the control experiment is generally similar to the configuration of the calculations carried out in [5]. In that work, similar biases of the diffuse thermocline were encountered. Now, to overcome them, let us pay attention primarily to the mechanism of model heat diffusion. Since the adiabatic heat transport processes dominate in the region under consideration, the parameterization of their explicitly not resolved component is of great importance. It is necessary to describe the process of the available potential energy (APE) transfer into the kinetic energy of subgrid eddies, in particular, into the subgrid component of tropical instability waves. As mentioned in [18], the baroclinic instability plays a significant role in their generation. Typically, these processes are described by using one or another modification of the Gent-McWilliams closure [19]. In the present work, we use this scheme in the skew formulation of [20] thanks to its numerical stability, as well as to the ease of its implementation together with the isopycnal diffusion scheme [21], which is acknowledged as more physically justified than the simple horizontal Laplace diffusion. In our case, the procedure [22] is applied to taper the components of the mixing tensor in regions where the slope of isoneutral surfaces exceeds the value of S max = 0.01.
Thus, experiment (c) differs from experiment (b) in that the horizontal Laplace diffusion scheme is replaced by the parameterization [20] with the same diffusion coefficient. As a result, the vertical temperature gradient of the thermocline increased to 0.075 ºC/m. Its depth decreased to 201 m, while the improvement in the position of the isotherms of 14-16ºC in the lower part of the considered depth interval is more noticeable. In (c) and subsequent plates, the solid lines show the results of the current experiment, while the dashed ones repeat the previous experiment for convenience of comparison.
In experiment (d), compared to the one (c), the Smagorinsky coefficient C of [13] was increased from 1 to 1.5. This led to a slight increase in the thermocline depth (by 1-2 meters), while the vertical gradient change was negligible. The plate does not allow one to clearly see these changes, but to preserve the sequence of experiments and for consistency with Figure 3, we will keep it.
Further, in experiment (e), the nominal (equatorial) value of the diffusion coefficient was increased from 300 to 500 m 2 /s. This led to a decrease in the slope of the isotherms in the upper part of the thermocline, which better corresponds to the observations. Besides, it is consistent with a qualitative In Figure 2, the same experiments are illustrated by the zonal temperature distribution near the equator on the 140ºW meridian. One can see an overestimation of the thermocline depth by the model, as well as the convexity of the upper isotherms, which apparently indicates an underestimated intensity of the heat exchange of the cold tongue with the surrounding waters. The most noticeable effect in correcting these biases under the conditions of our experiments was exerted by the transition to the parameterization [20] and by the decrease in the background vertical diffusion coefficient. The performed tuning of the model made it possible to bring the surface temperature at this point of the equator closer to the observational data, thereby improving the temperature contrast between the waters above and below the thermocline, which is important for the reproduction of the ENSO phases.

Conclusions
We performed a series of numerical experiments with the eddy-permitting INMIO model to study the sensitivity of the Pacific equatorial thermocline to the settings and characteristics of turbulent mixing parameterizations. Although the problem of a diffused thermocline is still open, the results of this work may help us better understand the mechanisms of this bias and possible ways to overcome them, as well as directions for improving the ocean research and forecasting model systems. The results of our calculations have shown that the utilization of the Gent -McWilliams eddy mixing parameterization and the Redi isopycnal diffusion in a combined formulation [20] reduces the thermocline reproduction errors by increasing its sharpness and decreasing the depth. It also allows one to adjust the slope of the isotherms. Calculations with a small background vertical diffusion have made it possible to improve the adiabaticity of the model solution and, as a result, to reduce the near-surface stratification biases. At the same time, a decrease in the maximum values of the diffusion coefficient has led to an improvement in the characteristics of the equatorial thermocline in the western part of the Pacific Ocean, but caused an unacceptable increase in the slope of the isotherms. Also, it was found that the solution is somewhat sensitive to the background viscosity parameters and the isopycnals slope tapering.