Carlos Carbonel
Universidad Nacional Mayor de San Marcos (UNMSM), Lima, Peru.
DOI: 10.4236/ojms.2026.164012   PDF    HTML   XML   6 Downloads   37 Views   Citations

Abstract

The dynamics of the upper ocean layer off Peru during El Niño Modoki 2002-2003 is investigated numerically using a vertically integrated $1\frac{1}{2}$ layer stabilized finite element model for the open and limited coastal ocean. The model solves the motion and continuity hydrodynamic equations coupled with the advection–diffusion transport equation for temperature. Solutions are forced by time dependent regional heat fluxes (short wave radiation and long wave back radiation) and wind data derived from the “European Centre for Medium-Range Weather Forecasts (ECMWF)”. The role of the regional wind and heat fluxes on the time dependent dynamics is evaluated. The problem investigated here is whether the regional forcings and the ocean hydro thermodynamical response in space and time are connected. Numerical solutions are found in an oceanic region between the equator and 30˚S and between 70˚W and 100˚W. The model is successful in reproducing time dependent features of El Niño 2002-2003, such as the Sea Surface Temperature (SST) field, which compares favorably with observed patterns in terms of spatial distribution during the time evolution of El Niño. A Mean Absolute Error of 0.56˚C is evaluated. Dynamically, the model results demonstrate that the SST field evolution during 2002-2003 is a causal effect of the regional forcings, mainly influenced by the heat fluxes. The results were possible due to a particular observed thermal feature of the study region: the rate of change of the Sea Surface Temperature and the net incident heat flux are linked during 2002-2003.

Keywords

Numerical Modelling, Regional Forcings, Heat Fluxes, Eastern South Pacific Ocean, Finite Elements

Share and Cite:

1. Introduction

In the eastern South Pacific Ocean (Figure 1), the simultaneous occurrence of colder waters near the coast (upwelling) and a southward warmer intrusion offshore (during El Niño) configures a scenario with strong SST gradients. The classical definition of the El Niño phenomenon is related to SST anomalies and the intrusion of warmer temperatures along the coastal eastern Pacific Ocean.

Figure 1. The study region (0˚S - 30˚S, 70˚W to 100˚W) is located in the eastern South Pacific Ocean.

But during 2002-2003, a non-canonical El Niño called El Niño Modoki took place, following reported contributions [1] [2]. Also, some authors describe it as an Eastward Central Pacific El Niño [3] [4].

Along the equator, during 2002-2003, El Niño Modoki is characterized by a warming of the NINO3 (150˚W - 90˚W; 5˚N - 5˚S) and NINO4 (150˚E - 150˚W; 5˚N - 5˚S) regions of the order of 1˚C. Otherwise, along the Perú coast, important SST anomalies were observed. During autumn 2002, SST anomalies of up to 2.5 ˚C appeared off Chicama (7.8˚S), diminishing towards the south and north. Also, it was possible to observe a reactivation of the upwelling in very coastal zones between 6˚S and 11˚S (Punta Falsa and Huacho). The 15˚C isotherm was located between 30 and 50 m, below its usual depth to the north of 9˚S [5]. In summer 2003 along the coast (3.5˚S - 12˚S), in January, the deepening of isotherms continued, but in February ascended because of the reactivation of coastal upwelling [6].

A challenge of El Niño Modoki is the fact that its generating mechanisms are still largely unknown [7]. Otherwise, it is well known that the Sea Surface Temperature is very dependent on radiation from the sun [8]. In the oceans, the SST variability is the dynamic result of the energy obtained from the sun and the current circulation in the deep and coastal ocean. A possible link between heat fluxes and warming during El Niño Modoki needs to be reviewed carefully.

In particular, the exchange of heat fluxes during warming El Niño events in the Pacific Ocean was not well discussed in the literature. But an important possibility associated with the warming effects during El Niño was indicated by Weare in 1983 [9]. Additionally, a study analysis [10] supports the view that surface heat flux variations have a role in the meridional warming of SST during ENSO development. This new conclusion differs from previous studies (e.g., [11]). Also, the comparison of East-West Short Wave Radiation (SWR) and SST [8] indicates that in the Eastern Pacific the SWR is higher during the initiation of El Niño, reaching higher SST.

Also, the influence of atmospheric heat fluxes has not been considered in numerical simulations in its real magnitude. But some signs of its possible influence were observed in some works. In the pioneering simulation of Philander and Siegel (1985) for El Niño 1982-1983 [12], constant surface heating was considered (short wave and long wave back radiation), obtaining high SST distributions with intrusion of warm waters off the South American coast. Moreover, Jia, Wells and Rowe (1990) [13] also simulated El Niño 1982-1983 but without heat input at the surface from the atmosphere. Their results show the persistence of cold waters off the South American coast (which was not understood), indicating implicitly that an important forcing was not included.

The aforementioned signals, associated with atmospheric heat fluxes, are a strong reason to explore the influence of atmospheric forcings, including solar radiation, on the evolution of sea surface temperature and currents in the eastern Pacific Ocean. The warming processes of the eastern Pacific Ocean during El Niño Modoki were not studied dynamically and need to be reviewed carefully. The present paper is a little step in this direction.

In the present study, the role of only time dependent regional atmospheric forcings will be investigated numerically, evaluating the time evolution of the resulting SST variability in space and time and the current fields in the domain. In this way, the importance of the regional atmospheric forcings could be really understood. The numerical study is performed covering the region between the equator and 30˚S and between 100˚W and 70˚W (Figure 1) with an unstructured finite element mesh (Figure 2).

The use of complicated models can obscure the basic physical processes at work in them, whereas the adequate formulation of simple models can show clearly the importance of particular processes. Here, a simple model ($1\frac{1}{2}$ layer) of limited area domain is presented to describe the upper layer dynamics off Perú during El Niño Modoki 2002-2003. Numerically, the presented model has some technical advantages. The finite element Petrov-Galerkin formulation model uses a stabilizing operator in space and time to improve the classical Galerkin approach [14] [15]. This methodology allows accurate solutions, precluding the use of artificial filters, mass diffusion and dampers, commonly used in several popular models. The numerical solutions are obtained strictly from the finite element formulation of the governing equations in the domain. An additional advantage is the treatment of the open boundaries which is based on weakly reflective open boundary conditions derived from the characteristic equations [16]. The effectiveness of those kinds of open boundary conditions for coastal circulation problems has been previously reported [17] [18].

In the present paper, some of the results are the following: the model forced by regional forcings is able to describe the time evolution of the SST field which compares favorably with observations. The results verified the connection between the forcings and the thermal response of the regional domain.

The article is organized as follows: Section 2 describes the model and its parametrizations. Section 3 presents the results of the modelling of ocean dynamics during 2002-2003 and Section 4 presents the findings of the study.

2. The Model

The used ocean model (a gravity reduced $1\frac{1}{2}$ layer model) has one active layer of density $\rho$ (entire layer above the thermocline), overlaying a deep layer of density ${\rho }^l$ , where horizontal gradients are assumed to vanish [15]. A Cartesian coordinates system (x, y) oriented positively to the East and North directions respectively is considered. The non-permanent vertical integrated upper layer governing equations of motion, continuity and heat

$\partial u/\partial t+u\nabla u+g\left(\sigma \nabla h+\left(h\theta /2\overline{\gamma }\right)\nabla T\right)+fk\times u-A_H{\nabla }^{2}u-\left({\tau }^s-{\tau }^i\right)/\rho h=0$(1)

$\partial h/\partial t+u\nabla \left(uh\right)-w_e=0$ (2)

$\partial T/\partial t+u\nabla T-K_H{\nabla }^{2}T-\left(Q_S-Q_I\right)/\rho c_{p}h=0$ (3)

The governing Equations (1), (2) and (3) describe the time evolution of the velocity $u=\left(u,v\right)$ , the upper layer thickness h and temperature T of the upper active layer in each point of the spatial domain. The Nabla operator is represented by $\nabla =\left(\partial /\partial x,\partial /\partial y\right)$ . The lower layer has a constant temperature T^l. The upper layer describes the entire layer above the thermocline and due dynamical reasons can become small during upwelling motions and also is able to deepen during downwelling events.

In the above equations, Q_S is the heat input at the ocean surface, Q_I is the heat loss at the interface, ${\tau }^s=\left({\tau }_x^s,{\tau }_y^s\right)$ is the wind stress at the surface, ${\tau }^i=\left({\tau }_x^i,{\tau }_y^i\right)$ is the stress at the interface, $k$ is the unit vector in the vertical direction and f is the coriolis parameter (evaluated in function of the spatial coordinates x, y). The entrainment velocity is w_e, c_p is the specific heat and the parameter $\sigma =1-\gamma$ where $\gamma =\rho /{\rho }^l$ and $\overline{\gamma }=\gamma /\left(\gamma -\sigma \right)$ . The upper layer density varies according to $\rho ={\rho }^{\text{l}}\left[1-\alpha \left(T-T^l\right)\right]$ where α is the coefficient of thermal expansion. The entrainment velocity w_e cools the upper layer and prevent the interface of surfacing in regions of intense upwelling [19],

$w_e=\left[{\left(H_e-h\right)}^2+\left(H_e-h\right)\left|H_e-h\right|\right]/2t_e{H}_e$ (4)

where H_e is the entrainment referential depth and t_e is a specified relaxation time.

To completely define the proposed ocean model we need to characterize the applied external sources of energy, the time dependent boundary conditions and the initial state from which the relevant hydrothermodynamical response will be obtained. The hydrodynamics of the model is forced by the wind stresses evaluated from the wind velocity fields using the quadratic relations ${\tau }_x^s={\rho }_{air}{c}_D{W}_x\left|W\right|$ , ${\tau }_y^s={\rho }_{air}{c}_D{W}_y\left|W\right|$ and $c_D=\left(0.509+0.065\left|W\right|\right)/1000$ [20]. The terms ${\tau }_x^i,{\tau }_y^i$ , are parameterized equal to ρhru and ρhrv respectively, where r is a constant.

The thermodynamic of the model is forced by the short wave radiation Q_sw and the resultant outgoing longwave radiation Q_lw. The heat balance of the model is associated to the surface heat input Q_S defined as

$Q_S=Q_{sw}-Q_{lw}-Q_{Lat}-Q_{Sen}$ (5)

where Q_Sen is the sensible heat and Q_Lat is the latent heat. The sensible heat is calculated from the bulk formula $Q_{sen}={\rho }_{air}{c}_{pair}{C}_S\left(T-T_{air}\right)\left|W\right|$ , where c_pair is the specific heat of the air, C_S an empirical coefficient, and T and T_air are the sea surface temperature and air temperature respectively. The latent heat is estimated from $Q_{Lat}={\rho }_{air}{C}_{L}L\left(q-q_{air}\right)\left|W\right|$ , where L is the latent heat of vaporization, C_L an empirical coefficient and q, q_air are the saturation specific humidity of the sea surface and air respectively. For a constant lower layer temperature T^l, the heat flux across the interface is given by:

$Q_I=k_I\left(T-T^l\right)/h$ (6)

in this way the gain or loss of heat across this interface, depends on the dynamical convergence or divergence of the flows $\omega$ in the upper layer, as expressed by the parameter $k_I=\omega H$ .

Concerning the boundary conditions it is assumed that land type boundaries give rise to the classical non-slip conditions

$u=v=0$ (7)

and for upper thickness and for the temperature, homogeneous conditions are prescribed. On the sea side boundaries we prescribe a weakly reflective open boundary condition, based on the characteristic method. This is done assuming that on a normal direction (x_n) to the open boundary

$\text{d}{R}^{\pm }/\text{d}t=0$ (8)

where the invariant $R^{\pm }=u_n\pm 2c$ are valid, along the characteristics ${\left(\text{d}x/\text{d}t\right)}^{\pm }=u_n\pm c$ normal to the open boundary. u_n represents the velocity component normal to the boundary and $c=\sqrt{g\sigma h}$ is the celerity wave. The weakly reflective conditions on the open boundary are defined by the in-going characteristic of the presented equations. For the temperature at the open boundaries, homogeneous condition is assumed.

Symmetric forms are always desirable because they possess better stability properties, an important feature to be used in the design of approximate solutions of advection dominated problems and shallow water wave problems [21].

For the upper layer dynamics problem treated in this paper a complete symmetrization is not achieved. Nevertheless, the hydrodynamic part of the equations governing the ocean circulation could be symmetrized. Furthermore, the whole set of unknown variables retains the same physical dimensions, which means that the transformed equations are adequately scaled and a new set of variables are defined as $d=2\sqrt{g\sigma h}$ and $b=c^{0}T/T^u$ which have the dimensions of velocity. Here T^u is a referential upper layer temperature and $c^o=\sqrt{g\sigma H}$ is a referential celerity. After the change of variables, the set of equations for the unknown vector V could be written as:

$\underset{\_}{I}\partial V/\partial t+\underset{\_}{A}\partial V/\partial x+\underset{\_}{B}\partial V/\partial y+\underset{\_}{C}V+\underset{\_}{D}\left({\partial }^{2}V/\partial x^2+{\partial }^{2}V/\partial y^2\right)+F=0$ (9)

where

$\underset{\_}{A}=\left(\begin{array}{cccc}u& 0& c& s\\ 0& u& 0& 0\\ c& 0& u& 0\\ 0& 0& 0& u\end{array}\right)$ , $\underset{\_}{B}=\left(\begin{array}{cccc}v& 0& 0& 0\\ 0& v& c& s\\ 0& c& v& 0\\ 0& 0& 0& v\end{array}\right)$ , $\underset{\_}{C}=\left(\begin{array}{cccc}0& -f& 0& 0\\ f& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$ (10)

$F=\left(\begin{array}{c}-{\tau }_x/\rho h\\ -{\tau }_y/\rho h\\ -w_{e}g\sigma /c\\ \left(Q_I-Q_S\right)c^o/\rho c_{p}h{T}^u\end{array}\right)$ , $V=\left(\begin{array}{c}u\\ v\\ d\\ b\end{array}\right)$ , $\underset{\_}{I}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$ (11)

$\underset{\_}{D}=\left(\begin{array}{cccc}-A_H& 0& 0& 0\\ 0& -A_H& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -K_H\end{array}\right)$ (12)

and $s=gh\theta T/\left(2\overline{\gamma }{c}^o\right)$ .

A space-time finite element partition is considered. The time interval [0, T_f] is divided into time sub-intervals $I_n=\left(t_{n+1}-t_n\right)=\Delta t$ , for $t\in \left[0,T_f\right]$ , where the time levels t_n, t_n+1 belong to an ordered partition of time levels $t_0=0<t_1<\cdots <t_n<t_{n+1}<\cdots <t_F=T_f$ .

The spatial domain Ω is partitioned into N sub-domains Ω_e. For each level of time n, the domain space-time integration is the Slab $S_n=\Omega \times I_n$ with border ${\Gamma }_n=\Gamma \times I_n$ , so that the Slab is composed of N space-time elements $S_n^e={\Omega }_e\times I_n$ .

A finite-dimensional trial solution space $U_n^h$ and a test function space ${\widehat{U}}_j^n$ are defined. For the presented model, continuous interpolation in space and time is adopted. The Petrov-Galerkin approximate solution for the ocean model problem is the vector $V^h\in U_n^h$ (of piecewise polynomials) which satisfies the following variational formulation for all ${\widehat{V}}^h\in {\widehat{U}}_n^h$

${\displaystyle {\int }_{S_n}{\widehat{V}}^h{R}^h\text{d}\Omega \text{d}t}+{\displaystyle {\sum }_{e=1}^N{\displaystyle {\int }_{S_n^e}{\overline{\Psi }}_e{G}^h{R}^h\text{d}\Omega \text{d}t}}=0$ for $n=1,2,3,\cdots$ (13)

where the residual vector R^h is

$R^h=\underset{\_}{I}\partial V^h/\partial t+\underset{\_}{A}\partial V^h/\partial x+\underset{\_}{B}\partial V^h/\partial y+\underset{\_}{C}{V}^h+\underset{\_}{D}\left({\partial }^2{V}^h/\partial x^2+{\partial }^2{V}^h/\partial y^2\right)+F$ (14)

and the space-time operator G^h is

$G^h=\underset{\_}{I}\partial {\widehat{V}}^h/\partial t+\underset{\_}{A}\partial {\widehat{V}}^h/\partial x+\underset{\_}{B}\partial {\widehat{V}}^h/\partial y$ (15)

For an equation system, the ${\overline{\Psi }}_e$ definition is even more complicated due to the presence of multiple wave components. For the ocean problem focused in this paper

${\overline{\Psi }}_e=1/\sqrt{\left(1/\Delta t^2\right)\underset{\_}{I}+\left(1/l_e^2\right)\left\{{\underset{\_}{A}}^2+{\underset{\_}{B}}^2\right\}}$ (16)

where $l_e=\sqrt{{\left(\partial {\xi }_i/\partial x\right)}_{\text{max}}^2+{\left(\partial {\xi }_i/\partial y\right)}_{\text{max}}^2}$ is the representative length of the triangular spatial element Ω_e, and ξ_i are the local co-ordinates. The first integral in the variational formulation (13) is the Galerkin term. The second integral represent the Petrov-Galerkin contribution, where ${\overline{\Psi }}_e$ contain the intrinsic time scale free parameters [22]. The resulting system of algebraic equations is solved using the direct method of Gauss (L-U descomposition). The program code is written in Fortran.

3. Numerical Experiments

Numerical experiments are performed to clarify the relative importance of regional forcings in the resulting time dependent evolution of the upper layer dynamics.

The study area extends in the north-south direction from equator to the 30˚S and in the east-west direction from the 70˚W to 100˚W. This region was represented by an unstructured mesh of 5684 surface triangular elements and 2994 nodes, as illustrated in Figure 2.

Near the coastal boundaries, a finer mesh was employed which gradually becomes coarser and coarser into the offshore region. The smaller elements are located along the coastline.

3.1. Model Forcings

Ten days averaged meteorological data are extracted from ECMWF (ERA-Interim), on a 0.5˚ × 0.5˚ grid to determine the forcing functions of wind velocities W_x, W_y and Q_S (Q_S = Q_sw − Q_lw − Q_lat − Q_sen) for the $1\frac{1}{2}$ layer model. To calculate the forcing function of the model, it is used ten days averaged fields for the surface incomming short wave radiation Q_sw, the outgoing longwave radiation Q_lw, the wind velocity components W_x, W_y the air temperature T_air and derived air specific humidity q_air.

Figure 2. Finite element mesh of the study area between the equator and 30˚S and between the 70˚W to 100˚W.

For the finite element mesh, the discrete fields W_x, W_y, Q_sw, Q_lw, T_air, q_air used as input for the model, were obtained interpolating linearly the data at each node of the FEM mesh. It is necessary to remark that the surface heat fluxes Q_Sen and Q_Lat, were obtained using the upper layer water temperature T calculated by the model. The observed SST fields are used only in the subsections 3.4 and 3.6 for comparison with the calculated results, evaluating in this way the performance of the model. There are not remote forcings in the model.

3.2. Parameters of the Model

The parameters of the model are mostly physically realistic values. The referential densities of the upper, lower layer and air ρ^u = 1022.5 kg·m⁻³, ρ^l = 1025.5 kg·m⁻³ and ρ_air = 1 kg·m⁻³, are representative of observed conditions in the region. The coefficient A_H is 600 m²·s⁻¹, K_H is 150 m²·s⁻¹ and r = 0.7 × 10⁻⁵. The adopted entrainment time scale t_e = 1/4 day, ensures that the upper layer thickness does not vanish in coastal upwelling areas. The parameters C_L, C_S for the evaluation of the heat fluxes Q_Lat and Q_Sen, range in the literature from 0.0008 to 0.002 [23]-[25], here were taken equal to 0.0008. The specific heat of the air is c_pa = 1004 J·˚C⁻¹·kg⁻¹. The saturation specific humidity q and q_air are evaluated in function of the saturation vapor pressure e and e_air [25]. The latent heat of evaporation is fixed equal to L = 2.44 × 10⁶ J/kg. The value of the coefficient of thermal expansion is α = 2.4 × 10⁻⁴ ˚C⁻¹. The time step is defined as Δt = 3 hrs.

3.3. Spin up of the Model. Initial Conditions

The initial conditions of ocean real flows is a modelling problem. In the present paper, to obtain realistic dynamic initial conditions at time t = 0, for the field variables, the ocean model was spin up from a state of rest and constant h, T fields and integrated forward in time for a period of 1 year using the monthly forcing fields from January to December of a canonical year (Similarly to the previous reference article [15]). The fields at the end of the year were taken to initialize the simulation of El Niño 2002-2003 which was started at the beginning of January 2002 and continued up to December 2003.

3.4. Simulation of SST Changes during 2002-2003

The model simulation is started in the month of January 2002. Figures 3-8 show the SST solutions at various stages from 2002 until July 2003. In those figures, the left column shows the instantaneous fields calculated in the middle of the month, while the right column shows the observed fields corresponding to the second set of data (10 days) of observed data for each month.

The resulting instantaneous patterns of the temperature are similar to the observed mean patterns. The model simulation results reproduce favourably the time evolution of the observed SST fields during El Niño 2002-2003. It is important to indicate that the model is able to grasp the main aspects of the upwelling dynamics after the warming phases. The cool water along the coast arises by the divergence of surface waters (see e.g. mid-June, mid-July 2002, Figure 4 and Figure 5) forced by the increasing wind velocities.

Figure 3. Model SST (left) and observed SST (right) from January to March 2002 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figure 4. Model SST (left) and observed SST (right) from April to June 2002 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figure 5. Model SST (left) and observed SST (right) from July to September 2002 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figure 6. Model SST (left) and observed SST (right) from October to December 2002 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figure 7. Model SST (left) and observed SST (right) from January to March 2003 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figure 8. Model SST (left) and observed SST (right) from April to June 2003 (arranged from top to bottom). The SST has a contour interval of 0.5˚C.

Figures 9-11 show the velocity fields and surface layer thicknesses on certain specific days. It can be observed that the greatest layer thicknesses are always located in the southwest region. Exceptionally, in the summer of 2002, near the coast, a deepening of the surface layer is also observed, which is associated with surface countercurrents. Along the coast, the velocity field shows the coastal boundary layer effect (Figure 9).

Many models are better in reproducing warm conditions or tend to be warmer than the normal conditions. But, the presented $1\frac{1}{2}$ layer model shows a gradual change from warm to cool phase without numerical instabilities, due particular advantages of the numerical formulation, responding directly to forcings model.

Figure 12 illustrates the warming response of the marine region. The ten-day means of the observed SST in the study region and the integrated means of the instantaneous calculated model temperature every 5 days, during El Niño Modoki 2002-2003, are represented. The calculated model results follow the time variability of observation means, validating the time evolution of the thermal behavior in the region. The results show no numerical problems that can cause divergence of the time dependent solution. There is a temperature difference between model and observed values, particularly during the warming phases, which could be a consequence of not including forcings and interactions. The Mean Absolute Error (MAE) of the calculated time series is equal to 0.56˚C.

The comparison (Figure 12) verifies dynamically that the SST changes are linked to the regional wind and heat fluxes forcings.

Figure 9. Some instantaneous calculated maps of velocity field V and upper layer thickness field h at around mid-March 2002.

Figure 10. Some instantaneous calculated maps of velocity field V and upper layer thickness field h at around mid-June 2002.

Figure 11. Some instantaneous calculated maps of velocity field V and upper layer thickness field h at around mid-August 2002.

Figure 12. The integrated ten-day means of the observed Temperature (T_obs from ECMWF data source) and model instantaneous solution of Temperature (T_cal) every 5 days for the study area during 2002-2003.

3.5. A Regional Thermal Feature during El Niño 2002-2003

A particular thermal feature was observed in the study region (0˚S to 30˚S and 100˚W to 70˚W) during 2002-2003. Based on the extracted observations (ECMWF ERA-Interim), Figure 13 shows the integrated monthly means of the rate of change of observed sea surface temperature (∂T/dt) and the net incident heat flux (Q_n = Q_sw − Q_lw) for the study area.

Figure 13 suggests a link between Q_n and ∂T/∂t, insinuating the importance of atmospheric forcing through the incident heatflux Q_n during the warming changes in the eastern coastal Pacific Ocean region. The mentioned feature was a motivation of the present study.

Figure 13. The rate of change of observed sea surface temperature (∂T/dt) and the net incident heat flux Q_n, for the study area during January 2002 to December 2003 (ECMWF-ERA Interim data source).

3.6. The Role of the Initial Conditions

Figure 14. Case with a colder initial field state. The integrated ten-day means of the observed Temperature (T_obs from ECMWF data source) and model instantaneous solution of Temperature (T_cal) every 5 days for the study area during 2002-2003.

To determine the role of the initial conditions on the simulation, a test is carried out. For this purpose, the initial conditions (section 3.3) used in the main run were modified, changing only the temperature field by a factor of 0.8. That is, a colder initial SST field was imposed. Figure 14 shows the calculated integrated SST of the region compared to the observations during 2002-2003. The model takes about 3 months to get a better approach to the observed SST means of the region. The Mean Absolute Error (MAE) of the calculated time serie is equal to 0.61˚C. The results indicate that is necessary a very carefully definition of the dynamical pattern in the initial conditions.

4. Summary and Conclusions

A vertically integrated $1\frac{1}{2}$ layer stabilized finite element model has been considered to understand the surface variability produced by only regional forcings during El Niño 2002-2003 in the eastern South Pacific Ocean. The model includes the hydrodynamics and thermodynamics of the upper ocean layer in a limited-area domain. The formulation uses a stabilizing operator in space and time to improve the classical Galerkin approach. This methodology allows accurate solutions.

The model solutions are forced by ten-day-averaged regional meteorological data, extracted from ECMWF. The time evolution results for the SST field compare well with the observed SST evolution during 2002-2003. The southward migration of the warm water tongue to the south, during the two phases of El Niño, was reproduced. During the warming phases, southward velocity components are induced and a coastal countercurrent is developed.

The source of the warming in the model region is the net incident heat flux and the promising results were possible due to a particular observed thermal feature of the region: the rate of change of the Sea Surface Temperature (∂T/dt) and the net incident heat flux (Q_n) are linked during 2002-2003.

The evolution of the integrated SST of the model in the study region follows the time variability of integrated observed SST, validating the time evolution of the thermal behavior in the region. The results show no significant accumulation or loss of heat in the system that can cause divergence of the time dependent solution.

It could be remark that it has been developed a $1\frac{1}{2}$ layer ocean model that is numerically adequated to simulate many aspects of the time evolution of the coastal ocean response off Peru. Additionally, the model has a simple structure (one dynamic layer) to understand the surface variability. In this way, it was possible to determine the real role of the regional forcings off Peru. Particularly, the importance of the influence of the heat fluxes in the development and decline of the warming phases during the El Niño event 2002-2003.

The presented model of the coastal ocean off Peru generates the SST variabilities during El Niño Modoki 2002-2003 by means of a surface atmosphere-ocean interaction.

Acknowledgements

This research was supported by “Fondo para la Innovación, Ciencia y Tecnología (FINCYT)”, Grant No. ECIP-1-P-033-14 and sponsored by the Instituto del Mar del Peru. Thanks to C. Arellano and C. Quispe for the data support and to J. Tam and D. Gutierrez for the valuable comments and coordination.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Ashok, K., Behera, S.K., Rao, S.A., Weng, H. and Yamagata, T. (2007) El Niño Modoki and Its Possible Teleconnection. Journal of Geophysical Research: Oceans, 112, C11007.[CrossRef]
[2]	Wang, C. and Wang, X. (2013) Classifying El Niño Modoki I and II by Different Impacts on Rainfall in Southern China and Typhoon Tracks. Journal of Climate, 26, 1322-1338.[CrossRef]
[3]	Yu, J. and Kim, S.T. (2010) Three Evolution Patterns of Central-Pacific El Niño. Geophysical Research Letters, 37, L08706.[CrossRef]
[4]	Zhang, W., Wang, Y., Jin, F., Stuecker, M.F. and Turner, A.G. (2015) Impact of Different El Niño Types on the El Niño/Iod Relationship. Geophysical Research Letters, 42, 8570-8576.[CrossRef]
[5]	Moron, O. (2011) Oceanographic Conditions during Autumn 2002: Demersal Cruise RV olaya 0205. Informe Instituto del Mar del Perú, 38, 41-45.
[6]	Tenorio, J. andMoron, O. (2011) Oceanographic Aspects of the Northern Peruvian Sea during Summer 2003: RV olaya 0301-02 Cruise. Informe Instituto del Mar del Perú, 38, 53-58.
[7]	Liang, X.S., Xu, F., Rong, Y., Zhang, R., Tang, X. and Zhang, F. (2021) El Niño Modoki Can Be Mostly Predicted More than 10 Years Ahead of Time. Scientific Reports, 11, Article No. 17860.[CrossRef] [PubMed]
[8]	Guo, B., Subrahmanyam, M.V. and Huang, Y. (2021) Shortwave Radiation and Sea Surface Temperature Variations over East and West Tropical Pacific Ocean. Open Access Library Journal, 5, 1-9.
[9]	Weare, B.C. (1983) The Possible Link between Net Surface Heating and El Niño. Science, 221, 947-949.[CrossRef] [PubMed]
[10]	Colasacco-Thum, N. (2015) An Investigation of Surface Heat Fluxes during El Niño-Souther Oscillation(Enso) Evolution in Reanalyses. Master’s Thesis, University Wisconsin-Madison.
[11]	Wang, W. and McPhaden, M.J. (2000) The Surface-Layer Heat Balance in the Equatorial Pacific Ocean. Part II: Interannual Variability. Journal of Physical Oceanography, 30, 2989-3008.[CrossRef]
[12]	Philander, S.G.H. and Seigel, A.D. (1985) Chapter 33 Simulation of El Niño of 1982-1983. Elsevier Oceanography Series, 40, 517-541.[CrossRef]
[13]	Jia, Y.L., Wells, N.C. and Rowe, M.A. (1990) A Simulation of the 1982-1983 El Niño. Journal of Geophysical Research: Oceans, 95, 5395-5403.[CrossRef]
[14]	Carbonel H., C.A.A., Galeão, A.C. and Loula, A.D. (2000) Numerical Study of Petrov-Galerkin Formulations for the Shallow Water Wave Equations. Journal of the Brazilian Society of Mechanical Sciences, 22, 231-247.[CrossRef]
[15]	Carbonel H, C.A.A. and Galeão, A.C.N.R. (2010) A Finite Element Model for the Ocean Circulation Driven by Wind and Atmospheric Heat Flux. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 32, 44-53.[CrossRef]
[16]	Verboom, G.K., Stelling, G.S. and Officier, M.J. (1982) Boundary Conditions for the Shallow Water Equations. Engineering Applications of Computational Hydraulics, 1, 230-262.
[17]	Carbonel, C.A.A.H. (2003) Modelling of Upwelling-Downwelling Cycles Caused by Variable Wind in a Very Sensitive Coastal System. Continental Shelf Research, 23, 1559-1578.[CrossRef]
[18]	Carbonel, C.A.A.H. and Galeão, A.C.N. (2007) A Stabilized Finite Element Model for the Hydrothermodynamical Simulation of the Rio De Janeiro Coastal Ocean. Communications in Numerical Methods in Engineering, 23, 521-534.[CrossRef]
[19]	McCreary,, J.P. and Kundu, P.K. (1988) A Numerical Investigation of the Somali Current during the Southwest Monsoon. Journal of Marine Research, 46, 25-58.[CrossRef]
[20]	Enriquez, A.G. and Friehe, C.A. (1997) Bulk Parameterization of Momentum, Heat, and Moisture Fluxes over a Coastal Upwelling Area. Journal of Geophysical Research: Oceans, 102, 5781-5798.[CrossRef]
[21]	Shakib, F. (1989) Finite Element Analysis of the Compressible Euler and Navier-Stokes Equations. Ph.D. Thesis, Stanford University.
[22]	Hughes, T.J.R. and Mallet, M. (1986) A New Finite Element Formulation for Computational Fluid Dynamics: III. The Generalized Streamline Operator for Multidimensional Advective-Diffusive Systems. Computer Methods in Applied Mechanics and Engineering, 58, 305-328.[CrossRef]
[23]	Bunker, A.F. (1976) Computations of Surface Energy Flux and Annual Air-Sea Interaction Cycles of the North Atlantic Ocean. Monthly Weather Review, 104, 1122-1140.[CrossRef]
[24]	Kraus, E.B. and Businger, J.A. (1994) Atmosphere-Ocean Interaction. Oxford Monographs on Geology A. Oxford University Press.
[25]	Rosenfeld, L.K., Schwing, F.B., Garfield, N. and Tracy, D.E. (1994) Bifurcated Flow from an Upwelling Center: A Cold Water Source for Monterey Bay. Continental Shelf Research, 14, 931-964.[CrossRef]