Oceanic general circulation models are robust and reliable tools for studying the circulation and dynamics of the oceans (Grötzner et al., 1999; Zavala-Hidalgo et al., 2003; Meza-Padilla et al., 2019; Morey et al., 2020; Olvera-Prado et al., 2023). Oceanic general circulation models are robust and reliable tools for studying the circulation and dynamics of the oceans. The characteristics of the numerical models allow for obtaining data continuously. Depending on the physics the models can resolve and the available computing resources to execute them, it is possible to obtain long-term data (several decades) with higher spatial and temporal resolutions than those typical in observations. Specifically, free-running simulations provide a physically consistent description of ocean dynamics that satisfies the primitive equations, allowing a better understanding of ocean dynamics (Morey et al., 2020). This work used the outputs of a 21-year free-running simulation of the GoM hydrodynamics performed with the HYbrid Coordinate Ocean Model (HYCOM), whose characteristics can be reviewed in Olvera-Prado et al. (2023). The simulation covers the period 1992-2012; the horizontal domain covers the region [98°W, 77°W] × [18°N, 32°N], with a 1/25° of horizontal resolution and 27 hybrid vertical layers. The validation of such numerical simulation showed that the simulated GoM hydrodynamics and hydrography are consistent with those obtained using ocean altimetry and in-situ observations (Olvera-Prado et al., 2023).

The most important variability modes were estimated from daily isopycnic layer-thickness anomalies of a surface and a deep layer via the Empirical Orthogonal Functions (EOFs) (Storch and Zwiers, 1999). The surface layer comprises the layers with density values from 17.25 to 26.52, ranging from 0-250 m approximately; the deep layer comprises the isopycnal 27.74 and ranges from 1000-1800 m approximately. The isopycnic layers’ selection is based on the representation of the deep GoM basins as a two-layer system described by Hamilton et al. (2018), in which the 6°C isotherm divides the upper (0− ∼1000 m) and lower (below ∼1000 − 1200 m) layers. In this work, the combined thickness of the surface and deep layers does not cover the total depth of the ocean, so they are not in physical contact; the reason for doing this was to prevent the surface and deep layer fluctuations from mirroring each other, allowing exploring their vertical connection. In addition, the thickness choice of the surface and deep layers made it possible to capture the most energetic component of the variability in each layer, accentuate thickness fluctuations, and reduce the ocean bottom effects on the deep fluctuations.

The variability analysis focuses on the mesoscale, i.e., ocean signals with space scales of 50-500 km and timescales of 10-100 days (CTOH, 2013). Thus, the data were filtered to enhance the mesoscale signal using a low-pass Lanczos filter considering different cutoff periods and a two-dimensional convolution filter considering different space windows. The following study cases were analyzed:

However, only the results for Case 2 are shown since it is the most representative of the mesoscale. The trend and annual cycle were removed from the time series of layer-thickness at each grid point, and the correlation matrix was constructed for the surface and deep layers separately. The correlation matrix is the covariance matrix of the standardized anomalies of the data (Wilks, 2019). The variability analyses intend to identify the intercorrelation between the thickness anomalies in the same layer and between the surface and deep layers (characterized by the correlation matrix); they do not intend to identify the strongest variations of thickness anomalies (characterized by the covariance matrix). The following analyses were performed:

In the separated and coupled EOF analyses, the corresponding vector time series X → ( x → , t ) can be expanded in terms of its spatial patterns, e →   k ( x → ) , and their time-varying coefficients, α k ( t ) , also called the principal components (PCs),

where x → denotes the spatial coordinates, t denotes time, and k runs over the dimension of X → ( x → , t ) . Since e →   k is orthogonal to each other, the regions well explained by each spatial pattern are different. A common approach to identify the regions well explained by a particular set of EOFs { e →   m } is to compute its proportion of explained variance (PEV),

where var denotes variance. Similar to e →   k , PEV can also be displayed as a map. Since the PEV represents the proportion of the total variance of the data explained by a particular set of EOFs, the spatial patterns are significant only in the regions with large PEV; these regions are statistically meaningful and were used to describe the spatial variability of the southern GoM. Generally, the regions with large EOF amplitude coincide with those with large PEV.

The predictability timescales of southern GoM variability were estimated by analyzing the characteristic properties of the PCs of each EOF analysis using time series analysis techniques (Storch and Zwiers, 1999; Box et al., 2015; Brooks, 2019). It was supposed that the predictability skill of southern GoM variability depends on the persistence of the PCs. The PCs can be grouped into processes with long or short memory depending on the persistence of their current values. In long-memory processes, the current value tends to persist during long periods, resulting in processes with long-term persistence; in short-memory processes, the values persist during shorter periods (Storch and Zwiers, 1999; Brooks, 2019). Thus, the skillfulness of the persistence forecasts is greater for long-memory processes than for short-memory processes and is highly determined by their corresponding auto-correlation functions (Storch and Zwiers, 1999).

The predictability timescale for each mode of southern GoM variability was estimated using a measure of the decorrelation time τ of the corresponding PC for each EOF analysis,

where ρ (k) is the auto-correlation function of the corresponding PC at discrete lag k. This work supposes that autoregressive processes of order p, or AR(p) processes, can describe the PCs. AR(p) processes are an important class of time series models (Brooks, 2019), widely used in climate research because they represent discretized versions of ordinary differential equations (Storch and Zwiers, 1999) and because they are helpful in modeling and predicting a variable of interest using only its own and past values, without the necessity of specifying a theoretical model of its behavior (Brooks, 2019). With this approach, an AR(p) process was fitted to each PC following Box et al. (2015), the theoretical auto-correlation functions of the fitted AR processes were obtained by solving the Yule-Walker equations following Storch and Zwiers (1999), and their predictability timescales were estimated using Equation (3).

Trenberth (1985) proposed the decorrelation time given by Equation (3) as an indicator of the persistence timescale of the corresponding time series and its associated processes. DelSole (2001) showed that Equation (3) arises naturally in statistical sampling theory and provides an accurate predictability estimate for oscillatory auto-correlation functions. Buckley et al. (2019) used this definition of the decorrelation time to estimate the predictability timescales for sea surface temperature and upper-ocean heat content in the North Atlantic. Because τ is based solely on the local auto-correlation function, it can be interpreted as a lower bound on the predictability timescale since skillful predictions, at least as long as τ, can be made using linear regression, where the inclusion of additional predictors can extend the predictability timescale (Buckley et al., 2019).

The stochastic origin of variability is based on the stochastic null hypothesis of climate variability proposed by Hasselmann (1976), which says that the long-term (or low-frequency) variability of a system can be understood as the integrated response of it driven by high-frequency fluctuations. The stochastic climate model of Hasselmann can be derived via a statistical dynamical model in terms of a fast (e.g., atmosphere) and a slow (e.g., ocean) component of the climate system (Storch and Zwiers, 1999), resulting in that, to first order, the long-term variability can be explained by an AR(1) process.

The AR(1) process is commonly taken as the null hypothesis to explain long-term variability (Grötzner et al., 1999; Storch and Zwiers, 1999; Dommenget and Latif, 2002). Several hypothesis tests have been developed to investigate the adequacy of this null hypothesis in describing long-term variability. The following two hypothesis tests were used to investigate the stochastic origin of southern GoM variability, considering the PCs of each EOF analysis (subsection 3.4). This work fitted an AR(1) process to each PC following Box et al. (2015). The fitted AR(1) process represents the null AR(1) process, and its spectrum represents the corresponding null AR(1) spectrum, that is, the reference spectrum against which the PC spectrum was compared and used to perform the hypothesis tests.

This subsection analyzes the correlation between the variability of the surface and deep layers for different time lags between them. The objective is to investigate the timescales in which surface and deep variability are strongly connected, for which a coupled EOF analysis was performed considering lags from 1000 to −1000 days of the surface layer relative to the deep layer. For example, a lag of 250 days means that the surface layer is shifted forward in time 250 days relative to the deep layer. For the lags of ±1000 days, the resulting time series corresponds to a reduction of 13% in the total length of the analyzed data.

Figure 2 shows the contained variance in the first four modes considering different time lags of the surface layer relative to the deep layer (the results are very similar for the other study cases mentioned in subsection 2.1); the figure only shows the results when the set of modes is not degenerate according to the rule of thumb of North et al. (1982). The contained variance is taken as a measure of the coupling’s strength between surface and deep variability for a given lag: a high variance corresponds to strong coupling, whereas a small variance corresponds to weak coupling. The percentage of the total variance contained in the first four modes varies little (from 39-44%) for lags from 1000 to -1000 days of the surface layer relative to the deep layer; the southern GoM variability is very complex, resulting in a non-evident coupling between its surface and deep variability. The most intense coupling for a particular state of surface variability at a given time is found with the concurrent and quasi-concurrent states of deep variability. The coupling between surface variability and deep variability is slightly higher for previous states of deep variability than for subsequent states of it, which could suggest that deep variability precedes surface variability and that the deep layer is more likely to influence the subsequent evolution of the surface layer than the other way around. Nonetheless, since the suggested direction of information transmission between the layers was obtained through purely statistical analysis, additional specific analyses are required to establish causality relationships between the deep and surface layers.

In order to perform an analysis of coupled variability in the southern GoM resulting in coupled modes with the highest possible contained variance, a coupled EOF analysis with a zero-day lag between the surface and deep layers was considered. Given the chosen zero-day lag, the coupled analysis describes concurrent spatial variability between the surface and deep layers; it identifies surface and deep regions spatially connected and the characteristics of that connection.

The distribution of contained variance in each mode for the separated and coupled EOF analyses ( Table 1 ) indicates that southern GoM variability is a complex process. The total variance is distributed in many modes, and no single mode contains a considerable portion of it. For the separated analyses, the first mode contains a quarter of the total variance; the contained variance in the first mode is slightly smaller for the coupled analysis. The decrease in contained variance for the higher-order modes is smaller for the surface layer, then for the deep layer, and ultimately for the coupled layers; such behavior can be transferred to the ability to define the corresponding variability. Because the contained variance in the first four modes is significant for the separated and coupled EOF analyses, it is assumed that these modes describe the predominant southern GoM variability well and that the study of it can be adequately addressed via EOFs decomposition.

What characteristics would the spatial patterns obtained from the isopycnic layer-thickness field have? In the surface layer, the isopycnic layer-thickness contours are a good proxy of circulation since they closely follow the flow streamlines; thus, the EOFs reflect where the main circulation variability occurs. In the deep layer, the isopycnic layer-thickness contours are not a good circulation proxy (Cushman-Roisin et al., 1990; Olvera-Prado et al., 2023). Nonetheless, since the isopycnic layer-thickness structure throughout the water column responds to pressure gradients, layer-thickness variability (and thus the associated EOFs) is expected to reflect regions with high circulation variability in each layer. Moreover, due to the weak deep GoM stratification, a bathymetry-influenced deep variability, strongly governed by the conservation of potential vorticity, is also expected. Accordingly, a description of surface and deep circulation variability of the southern GoM is first presented ( Figure 3 ) and then referred to when describing the EOFs ( Figure 4 ).

Figure 3 shows a monthly climatology of vertically averaged currents in a surface layer (0-250 m) and a deep layer (1000-1500 m) for the 21 years of simulation. The brown lines in the deep layer maps represent planetary potential vorticity ( f / H ) contours, where f is the local Coriolis parameter, and H is the bottom depth. In the surface layer, the currents’ speed is approximately three times greater than in the deep layer, with significant changes in the circulation patterns throughout the year, in contrast to the more persistent deep circulation. The changes in surface circulation patterns (panels of the first column in Figure 3 ) can be described as follows:

The circulation patterns in the deep layer (panels of the second column in Figure 3 ) have small variations throughout the year, reaching their minimum intensity during summer. The circulation is aligned with the planetary potential vorticity contours with a cyclonic rotation sense: west of 95.0 ∘ W , the currents go from north to south, and east of 95.0 ∘ W , the currents go from south to north. The circulation in the southeast region, where the bathymetry is abrupt, is weak, with few changes throughout the year. The changes in deep circulation patterns can be described as follows:

Figure 4 shows the first four EOFs of the surface and deep layers for the separated and coupled EOF analyses. The EOF amplitude is shown in color scale, and yellow contours represent the PEV of each EOF. Due to the thickness data processing (subsection 2.1), the EOFs are assumed to represent the mesoscale variability in each layer adequately.

First, the separated variability modes are described (panels of the first and second columns in Figure 4 ). The first four EOFs in both layers are well-defined and recognizable, and together they occupy the entire basin. As expected, the EOFs identify regions with high circulation variability in both layers ( Figure 3 ); they rank variability modes based on their contained variance but are not expected to identify individual dynamical modes (Storch and Zwiers, 1999; Monahan et al., 2009). In the deep layer, the bathymetry influences the variability; layer-thickness variations align with the planetary potential vorticity contours and intensify along bottom slopes. The spatial variability modes are described below:

The spatial patterns of coupled variability (panels of the third and fourth columns in Figure 4 ) differ from those of separated variability (panels of the first and second columns in Figure 4 ), which could indicate a significant connection between the surface and deep layers. By comparing the spatial structure and time evolution of coupled and separated variability, it is possible to obtain insights into the dynamics of southern GoM variability throughout the water column. The first four EOFs in the surface layer (panels of the third column in Figure 4 ) are well-defined, recognizable, and significant (they have spatial regions with appreciable PEV). However, only the first two modes in the deep layer (panels ofthe fourth column in Figure 4 ) are significant; the third and fourth modes lack spatial regions with appreciable PEV, indicating that higher-order surface variability has no connection with deep variability. Results also show no vertical correspondence between surface and deep layer-thickness variations; the fluctuations of deep isopycnals do not mirror those of surface isopycnals. The spatial variability modes are described below:

The predictability timescales of southern GoM variability were estimated by considering the decorrelation time (subsection 3.2) of each PC for the separated and coupled EOF analyses. The fitting of an AR model to each PC resulted in an AR(5) for each. Figure 5 shows the predictability timescales for the first four variability modes of the surface, deep, and coupled layers; the predictability timescales of the PCs are representative of the characteristic timescales of the corresponding spatial variability. The shortest predictability timescales are for the surface layer, with values from 1.4 to 10.6 months, which reminds the dominant timescales of mesoscale variability in the southern GoM (Vázquez de la Cerda et al., 2005; Pérez-Brunius et al., 2013; Zavala-Sansón et al., 2017). The deep layer has the longest predictability timescales, with almost three times the scales of the surface layer, from 3.0 to 30.5 months. Deep variability is more persistent and changes in longer timescales than surface variability, in agreement with the temporal change of surface and deep currents (monthly climatology of currents shown in Figure 3 ).

The predictability timescales of the coupled layers depend on the coupling’s strength of the surface and deep layers. The predictability timescales for the coupled layers lie between those for the surface and deep layers but closer to those of the deep layer, with values from 3.4 to 27.2 months ( Figure 5 ). Such behavior indicates that the memory of the coupled system resides in the deep layer but only for the first variability mode. By coupling a system composed of a fast (the surface ocean layer in this study) and a slow (the deep ocean layer in this study) component, the predictability of the fast component can be extended by coupling it with the slow component (Hasselmann, 1976; Grötzner et al., 1999; Dommenget and Latif, 2002). Nevertheless, the intensity of the coupled system’s high-frequency internal noise could weaken the coupling’s strength and drastically reduce its predictability timescale. The predictability timescales for the coupled system have the following behavior ( Figure 5 ):

• First mode. The long predictability timescale for this mode suggests an intense surface-deep coupling, in which both layers evolve conjointly with a predictability timescale of 27.2 months. The predictability timescale in the surface layer was substantially extended by a factor of 2.6 compared with that of the separated surface layer, resulting in a very similar timescale to that of the separated deep layer. The predictability timescale in the deep layer is very similar to that of the separated deeplayer (30.5 months).

• Higher-order modes. The surface-deep coupling is not strong enough to extend the predictability timescales in the surface layer with respect to those obtained in the separated analysis. Such weak surface-deep coupling resulted in predictability timescales in the deep layer similar to those of the separated surface layer.

The previous analysis indicated a connection between the surface and deep layers in the southern GoM, in which the coupled variability (variability throughout the water column) is more similar to the deep variability, at least for the principal variability mode. In order to explore further the connection of coupled variability with surface and deep variability, the similarity between each coupled PC and each surface and deep PC was examined considering two standard similarity measures, the root mean square similarity and Pearson’s correlation, where 1 indicates the maximum similarity (Cassisi et al., 2012). Figure 6 shows the similarity analysis, which helps reveal the coupling’s strength between coupled and separated variability. The first three coupled PCs have a clear and dominant association with a separated PC. The coupled PC 1 is mainly associated with the deep PC 1. The coupled PC 2 is highly associated with the deep PC 2. The coupled PC 3 is mainly associated with the surface PC 2; from the previous analysis of spatial variability, it is an expected result since the third coupled mode in the surface layer ( Figure 4K ) is practically the same as the second separated mode in the surface layer ( Figure 4E ).

The connection between coupled and separated variability is also appreciable in the time domain. Figure 7 shows the time series of the first three PCs of coupled layers with their corresponding associated PCs of separated layers. Their temporal nature is described below.

• The PCs for the first variability mode of the coupled and deep layers ( Figure 7A ) are very similar, as are their corresponding spatial patterns ( Figures 4D, B ). The PCs evolve in long timescales, with negative or positive values maintained for extended periods, although their wiggles do not always coincide; they are very time persistent and have the longest predictability timescale among all the variability modes ( Figure 5 ).

• The PC for the second coupled variability mode evolves with more accentuated short-term variations and has a high similitude with the deep PC 2 ( Figure 7B ). Such short-term variations in both PCs resulted in a low time persistence and short predictability timescales ( Figure 5 ) of their associated spatial patterns ( Figures 4F, G, H ).

• The PC for the third coupled mode is practically the same as for the second surface mode ( Figure 7C ); their corresponding spatial patterns ( Figures 4K, E ) are also practically the same. Due to the evolution with high-frequency variations of these modes, their time persistence is low, and their predictability timescales are short ( Figure 5 ).

There is a similarity between the variability modes obtained in the coupled and separated layers analyses. Nevertheless, some coupled modes differ from all those obtained in the separated layers analyses, suggesting that these result from the surface-deep coupling. Examples of this coupling are the spatial pattern of the first mode in the surface layer ( Figure 4C ) and the patterns of the second mode in both layers ( Figures 4G, H ).

As a preamble to analyze the stochastic origin of southern GoM variability, the distinctive timescales of surface, deep, and coupled variability of the southern GoM were revisited using the spectra of their PCs ( Figure 8 ). Figure 8 shows the sample spectra (cyan lines) and multitaper spectra (blue lines) of the first four PCs for the separated and coupled EOF analyses, with their corresponding null AR(1) spectra (solid red lines) and rejection limits at the 0.01 level (dashed red lines). The multitaper method consistently estimates a process’s true PSD by averaging a set of modified periodograms obtained using several tapers as windowing functions, reducing spectral leakage (Thomson, 1982). In terms of the PSD of a process, the longer the memory of the process, the higher its values for frequencies approaching zero. The principal variability mode for the separated and coupled layers ( Figures 8A-C ) varies over long timescales. The PSD of the PC increases with decreasing frequency (a typical red noise process), leading to the associated spatial patterns ( Figures 4A-D ) being long-term persistent. The long-term persistence of the principal mode is stronger in the deep layer ( Figure 8B ) than in the surface layer ( Figure 8A ), showing the well-known fact that the deep layer evolves in longer timescales than the surface layer. The PSD increase of the principal mode of the coupled layers ( Figure 8C ) is very similar to that of the deep layer ( Figure 8B ), indicating that the timescales of coupled variability are similar to those of deep variability.

In agreement with the predictability timescales analysis ( Figure 5 ), the higher-order variability is not long-term persistent; the PSD of the associated PCs does not consistently increase with decreasing frequency ( Figures 8D-L ). In the surface layer ( Figures 8D, G, J ), a significant portion of the PSD is concentrated in the three-month to one-year period band, with a small decrease for periods longer than one year. The PSD increase in the deep layer ( Figures 8E, H, K ) has a small bump in the one-month to three-month period band and a pronounced decrease for periods longer than one year. For the coupled layers ( Figures 8F, I, L ), a large portion of the PSD is concentrated in the three-month to one-year period band, with a small decrease for periods longer than one year. In summary, higher-order coupled modes vary on similar timescales to higher-order surface modes but without an associated correspondence between the modes of the same order.

When considering coupled variability of the southern GoM, there is a separation between its principal and higher-order modes ( Figures 2 and 4 − 8 ). Its principal variability (timescales and long-term persistence) is mainly associated with its deep dynamics; the deep layer strongly determines the principal variability throughout the water column. Its higher-order variability has surface and deep dynamics characteristics, with no clear association with one of them.

Now it is the turn to explore the stochastic origin of the long-term variability of the southern GoM, considering the two hypothesis tests described in subsection 2.3. First, the analysis was performed considering the hypothesis test described by Wilks (2019). For this, the spectra of the PCs and the spectra of the AR(1) models fitted to each of them (the null AR(1) spectra of the PCs) shown in Figure 8 were used. All the PC spectra do not rise above their AR(1) rejection limits at frequencies lower than 0.5 cpy (periods longer than two years), but they do rise above them at frequencies higher than 0.5 cpy (periods shorter than two years). Thus, the simplest stochastic model adequately describes only the long-term variability of the surface, deep, and coupled layers in the southern GoM. The high thermal inertia of the isopycnal layers acts as a memory, integrating short-term random variations and carrying energy to the long term in the manner described by Hasselmann (1976), producing a reddened spectrum consistent with an AR(1) process at frequencies lower than 0.5 cpy. On the other hand, the PC spectra deviate from their corresponding null AR(1) spectra at frequencies higher than 0.5 cpy. Southern GoM variability in timescales from some weeks to a couple of years is strongly affected by the episodic occurrence of different mesoscale circulation processes (Vázquez de la Cerda et al., 2005; Pérez-Brunius et al., 2013; Zavala-Sansón et al., 2017; Guerrero et al., 2020), which could explain the PC spectra deviations described above.

The results of the hypothesis test proposed by Dommenget and Latif (2002) are shown in Table 2 , considering null AR(1) and AR(2) spectra and a rejection limit at the 0.05 level. The PSD increase of each first PC is consistent with an AR(1) process, indicating that their source of long-term variability can be associated only with short-term random variations. The principal variability mode of the surface, deep, and coupled layers agrees with the stochastic null hypothesis of Hasselmann (1976); they can be described using a linear dynamics approach in terms of a fast and a slow component, with the involved dynamical processes grouped in one of those components (Storch and Zwiers, 1999). The random origin of the principal variability of the southern GoM does not mean that it is not predictable at all; only a completely random process such as AR(0) is not predictable at all. The performed fitting of higher-order AR models to the PCs to estimate their predictability timescales ( Figure 5 ) does not invalidate the results concerning the stochastic origin of southern GoM variability, as these methodologies have different objectives. A precise fitting of an AR model to a time series attempts to describe most of its wiggles, but this does not necessarily correspond to improving the description of its dynamics.

The higher-order variability modes are not consistent with AR(1) processes ( Table 2 ). These modes contain a more complex variability structure than the dominant ones, linear dynamics cannot explain them, and higher-order AR processes are needed to describe them. The different sources of southern GoM variability can be involved in such behavior, like surface-layer processes (e.g., atmosphere-ocean fluxes, advection and circulation of adjacent regions, and entrainment of sub-layers), a range of circulation features with a strong interaction between them (Zavala-Sansón et al., 2017; Guerrero et al., 2020), and topographic effects ruling to a certain extent the circulation throughout the water column (Pérez-Brunius et al., 2013; Zavala-Sansón, 2019). However, the higher-order variability can be adequately described using AR(2) processes, indicating that their dynamics is somewhat more complex than a linear one.