A Feasibility Study of Three-Dimensional Empirical Orthogonal Functions From the NASA JPL Ocean General Circulation Model: Computing, Visualization and Interpretation

2.1 Temperature data from the NASA JPL OGCM

We use the monthly temperature data from the NASA JPL non-Boussinesq ocean general circulation model (OGCM) from January 1950 to December 2003 (Song and Hou (2006), Song et al. (2010)). The horizontal resolution of the JPL OGCM is –degree latitude-longitude and was later aggregated to 1-degree as in Shen et al. (2017). After aggregating the data, the files were separated into MatLab files according to their month and year. Figure 1 shows one file of 3D ocean temperature data for January 1950.

Figure 1

The data encompasses the entire global ocean with land represented by a nonnumber (NaN) grid cell and temperatures as a °C grid cell. There are a total of 32 depth layers whose depth boundaries are given in Table 1. As an example take the first boundary from Table 1 as 5. This means the top layer extends from 0 to 5 meters. Similarly, the second layer is from 5 to 10 meters, etc.

This paper uses the JPL non-Boussinesq (mass conserving) OGCM as it well portrays ocean dynamics. The model was based on both the S-Coordinate Rutgers University Model (Song and Haidvogel (1994)), and the Regional Ocean Model System (Shchepetkin and McWilliams (2003), Haidvogel et al. (2000)). The JPL non-Boussineq OGCM served as a resolution to the uncertainties in altimetry sea-surface-height (SSH) and ocean bottom pressure (OBP) in Boussineq approximations (Song and Hou (2006), Song et al. (2010)) identified by Huang and Jin (2002). Both SSH and OBP’s consistency within the model have been validated by Moon et al. (2013) and Song and Colberg (2011). Because the OGCM output is dependable in its representation of ocean dynamics, we find it appropriate for our EOFs.

2.2 Formulation of the space-time data matrix

The OGCM output is an array written in 3D space grid boxes separated monthly for each year. Our time spans from January 1950 to December 2003. If we take January as an example it will have 54 years from 1950 to 2003; the same can be said for February, March, ……, December. Because of seasonal variation, we compute EOFs for each month from January, February, ……, December therefore we will compute 12 sets of EOFs. Let us first consider the data on a single layer with a given depth and time. We use January in the following description.

Denote the data matrix of one year and layer by T_year,layer and the final space-time data matrix by X. The gridded temperature data (T_{longitude,latitude}) for the first year (Y₁₉₅₀) and first layer (l₁) form a 360 × 180 matrix shown below:

The row data of the matrix are ordered from south to north starting at latitude –89.5° and ending at latitude 89.5° totaling 180 columns. The column data are ordered according to longitude from east to west starting at longitude 0° ending at longitude 359° totaling 360 rows. Namely, the first 180 points in the first row of the matrix represent 89.5°S to 89.5°N for a longitude 0° as shown by the first row of the matrix. The next row will show the same latitudes, 89.5°S to 89.5°N, for the next longitude 1°. This continues until longitude reaches 359°, the last row of the data matrix. Flattening this matrix by rows gives us a row vector of length 360 × 180 = 64,800.

We transpose this row vector ${\overrightarrow{T}}_{Y_{1950},l_1}$ to a column vector ${\overrightarrow{T}}_{l_1,Y_{1950}}$ and aggregate all the column vectors from 1950 to 2003 to form a space-time data matrix for the first layer of January. The 2D space-time data matrix for a single layer, denoted by $X_{Ja{n}_{2D}}$, looks like the following:

Anomalies can be computed from this matrix and multiplied by the area weight of the grid box i depending on its latitude in radians φ_i. These area weights reflect the geometric sizes of their corresponding grid boxes.

The result of the multiplication is the area-weighted space-time anomaly matrix. Conveniently the structure of this matrix allows for the removal of land/NaN data by removing rows as every year contains the same grid point of land data. For example, if January of 1950 had a land point at T_0,–89.5, then January 1951–2003 will have the same land point at T_0,–89.5. Only after the removal of land data can the area-weighted space-time anomaly matrix be fed into SVD, or any other method of computing EOFs, and produce the 2D EOFs layer by layer.

We can build on the idea of computing 2D EOFs to compute 3D EOFs by stacking 32 layers on top of each other as shown below.

Here,

is the total number of 3D grid boxes for the OGCM output. Following the same way, we formulate our February space-time data matrix, and so on, since our EOFs are computed for each month from January to December. This space-time data matrix for January is one of the 12 matrices used to compute 3D EOFs as outlined in the next section (see equation 7).

4.1 Eigenvalues and variances explained

EOFs, PCs, SVD eigenvalues, and eigenvalues of the temporal covariance matrix were computed for each month from January to December. Altogether 12 sets of EOFs were generated. Figure 3(a) shows the first eigenvalue of the temporal covariance matrix for the weighted temperature anomalies for each month from January to December. These eigenvalues represent variances of the space-time data projection onto their corresponding EOF1 for each month. The peaks in January and August imply large variances in both months, relative to their neighboring months. Our examination of all the EOF patterns has shown that those two peaks are likely the manifestation of the relatively stronger equatorial upwelling in January and August.

Figure 3

Another important feature of Figure 3(a) is the obvious difference between January and December. The difference is around 3% and is likely an indication of the presence of a temporal trend in the space-time temperature data. If there were no temporal trend the annual cycle would prevail showing a 12-month cycle instead, and hence December’s variance would be similar in size to January’s. Visually this would be observed as a smoother transition between December to January. Later on, we further quantitatively verify our assumption here by solving the eigenvalue problem for the de-trended data. Indeed, the annual cycle of λ₁ appears.

Figures 3(b) and (c) show the scree plot for both peak months January and August. The scree plot plots every eigenvalue for their respective month as a percentage variance or cumulative variance. The purpose is to quantify the variance explained by each mode and to identify how many modes are needed to explain a significant amount of variance. The blue dotted line explains the percentage variance of each EOF mode relative to the total variance, i.e.,

The red dotted line is the cumulative percentage variance relative to the total variance, defined as follows

Comparing Figure 3(b) and (c), January shows about p₁ = 30% variance explained by the first mode while August shows around 33%. Their second modes drop to 11% and 6% for January and August respectively. Their third modes are around 5% for both January and August. Figure 3(b) shows that for January, eigenvalues λ₄ and λ₅ are close to each other and have a very small difference p₄ – p₅. EOF errors are large because they are proportional to the inverse of the difference of two neighborhood eigenvalues, i.e., 1/(λ₄ – λ₅), according to North’s rule of thumb (North et al. (1982)). This implies that the mode mixture between EOF4 and EOF5 is stronger, and the geometric pattern of EOF4 may not have a meaningful physical interpretation due to the uncertainty caused by the mode mixing. For August, this mode mixing occurs from EOF3. Therefore, in this paper, we focus on the physical interpretations of EOF1 and EOF2. Up to mode 2, the cumulative percentage variance q₂ reaches slightly over 40% for January and slightly below 40% for August. Meaning a large number of variances are explained by EOF1 and EOF2. Considering the case of historical data reconstruction, like Shen et al. (2017), a number of modes may be used such that 80% or 90% of variances is represented. For our study here on the NASA JPL OCGM, it takes about 42 modes to represent 90% of cumulative variance for both January and August.

4.2 Physical interpretation of EOF1: Equatorial upwelling

Figure 4 shows our 3D EOF1 for January and August at different depth levels. We intend to use these EOF patterns to explain the ocean dynamic features associated with EOF1, PC1, and λ₁. Recall that the first mode of January and August explains 30% and 33% variance respectively. This means that a third of the cumulative variance in the data for each month is explained by the first mode alone. 2D EOF results from area-weighted anomaly data often show that the first mode is either a manifestation of a temporal trend or ENSO (Shen et al. (2017)). Interestingly enough, the first mode of our 3D EOF1 is not only related to the temporal trend, but also the vertical structure of ocean dynamics known as the equatorial upwelling.

Figure 4

Equatorial upwelling can be observed in Figure 4 as the horizontal divide of cool temperatures a little north of the equator with warm temperatures straddling it. Weisberg and Qiao (2000) and American Meteorological Society (2012) define equatorial upwelling as the rising of water along the equator occurring typically from 200 meters to the surface. Near the surface, the equatorial surface waters are being deflected by wind and Ekman flow. Thus, the surface waters are pushed away from the equator, and colder deep ocean waters are pumped up or more formally the deep ocean water upwells (Wyrtki and Eldyn (1982)). Our EOF1 provides a more precise quantitative description of the upwelling. Figures 4 and 5 show a strong upwelling starting around the depth of 100 meters that stops around 400 meters (specifically, see Figure 4(c) – (h)), and Figure 5(c)–(d)).

Figure 5

It is important to mention that a variety of signals can be explored in Figure 4. Despite the focus on equatorial upwelling in this section, there are other notable signals such as the significant warming pattern from surface to 400 meters in the Southern Ocean between 60◦S and 30°S. This is a widely known result of ocean warming due to global warming (see the satellite altimetry data from Song and Colberg (2011)). This is one of many signals to dissect from the first 3D EOF mode, but we will mainly focus on equatorial upwelling for the remainder of this subsection.

The signal strength of the equatorial upwelling process varies with depths. Looking at the surface layer for both January and August in Figure 4(a)–(b), the signal seems relatively weak, if we treat our 3D EOF1 as the upwelling mode. January shows a small line of cooler waters a little west of 120°E to approximately 170°E. August is less obvious with a smaller diffused line of cool water between 130°E to 170°E. Although the line of cool waters is narrow and small on the map, the region still covers a large equatorial area because Figure 4 is in the 2D Mercator map projection.

The upwelling region gets larger in deeper layers, particularly at 150 meters we see a very clear divide between 120°E and 180°E in Figure 4. The strong upwelling continues at 200 meters for both months, but January shows a stronger signal as expected from its higher eigenvalue (see Figure 3). At 300 meters, the two warmer strips off the equator begin to fade away. The warm-cool-warm sandwich pattern of the upwelling signal starts to disappear and eventually has little reminiscence at 400 meters suggesting the vertical sandwich has a depth of about 400 meters. The length in the tropical Pacific in the east-west direction stretches about 13,000 kilometers from 120°E to 120°W. The “bread pieces” of warm waters in both south and north are thicker than the thin “cheese filling” in the middle, as indicated by Figures 4 and 5.

To show how the upwelling signal changes by depth from the perspective of cross-section, we have made a north-south cross-sectional cut along a meridional line of 160°E, shown in Figure 5. To first identify upwelling in Figure 5, we look for a cold section of water just north of the equator with warmer waters around it. Both January and August show cold waters at the surface between the equator and 10° north. Figure 5(a) and (b) suggest that upwelling stops before 500 meters as the 3D EOF1 is quite uniform and shows little spatial variation for layers deeper than 500 meters. These observations motivated us to focus instead on depths from the surface to 400 meters where we can see more clearly a sandwich pattern of an equatorial cold (the narrow blue region at the equator) wrapped by two pieces of off-equatorial warm waters (the red regions). At about 100 meters, the warmer sections are distinct, continuing until they reach their strongest reading at around 150 meters. The warmer patches start to become thinner at 200 meters with the southern half between 10°S and 0° almost entirely gone for August at 250 meters and January at 300 meters. There is a small region in the southern patch around 20°S that persists to 400 meters along with the northern warm patch at 20°N. The entire upwelling region is limited to a zonal bend between 20°S and 20°N. However, this is only the EOF1 pattern, which reflects a temporally mean state. For a particular month, the equatorial upwelling region may be much larger, reaching 30°S and 30°N (Figure 6).

Figure 6

Recall EOFs are just a statistical tool used to show spatial patterns. We can identify the dynamics within the spatial patterns, and quantify their variance using EOFs. Equatorial upwelling has not been widely explored in previous EOF calculations therefore it is difficult to compare these results to other calculations such as Meyers et al. (1999), Mu et al. (2018), and Shen et al. (2017), whose results are favorable to ENSO. However, equatorial upwelling is an observable dynamic as seen in Figure 6. By visualizing the OGCM data we are able to identify what dynamics to expect of equatorial upwelling in EOF1.

Figure 6 plots the OGCM output for the strongest upwelling months and years: January 1993, January 2003, August 1990, and August 2003. We can observe that upwelling for all dates occurs in both the Pacific and Atlantic oceans. The upwelling signal is the strongest in the western and the central tropical Pacific, in contrast to the ENSO warm pattern over the eastern and central tropical Pacific. The upwelling over the western tropical Pacific in January 1993 (see Figure 6(a)) has a quantified warm-cold-warm sandwich pattern as 20°C – 12°C – 20°C, covering from 30°S to 35°N.

Focusing on January 1993 in Figure 6 the equatorial upwelling patterns in the Atlantic are also limited to the western part of the Atlantic ocean. The southern part of the sandwich bread piece encompasses a small area. In the Indian Ocean, the northern part of the sandwich bread piece almost does not exist. The upwelling signal is much weaker than that in the Pacific and Atlantic. The same can be said for the upwelling in all other months shown in Figure 6 when examining the water temperature distribution at the depth of 200 meters.

The high-latitude green regions of Figure 6 show another vertical structure of ocean dynamics: Ocean ventilation. This is the process in which surface water, recently in contact with the atmosphere, sinks, transporting carbon and oxygen to deeper ocean layers (Khatiwala et al. (2012) Naveira Garabato et al. (2017). Regions such as the Labrador Sea are more sensitive to ocean ventilation as its cold winter air will cause surface waters to lose heat to the atmosphere becoming cooler, denser, less buoyant, and consequently sinking deeper (Marshall and Schott (1999) Koelling et al. (2022)). These Labrador sea waters (LSW) that have undergone such convection can be seen in Figure 7.

Figure 7

This cross-sectional figure is a zonal mean, i.e., an average from longitude 0 to 360°. We should expect to see homogeneous values for depths up to 2500 meters at latitudes relative to the Labrador sea according to Yashayaev and Loder (2017). Observing Figure 7 from 55°N to 60°N, we can identify LSW occurring from 50 meters to 3000 meters. This range is much deeper than what was suggested from recent papers like Koelling et al. (2022) and Yashayaev and Loder (2017).

The cross-sectional zonal mean map of Figure 7 also shows equatorial upwelling. The August upwelling signal at about 100 meters depth appears to be stronger than that of January according to Figure 7(c) and (d).

4.3 EOF2 and the spatial pattern of El Niño

Figure 8 shows the eigenvalue λ₂ of the temporal covariance matrix for every month. In contrast to λ₁ in Figure 3, there is no longer a large jump from December to January in λ₂. As previously mentioned the smooth transition is due to an annual cycle.

Figure 8

Figure 9 shows the EOF2 pattern at five selected layers: 5, 50, 100, 150, and 200 meters. The warm pattern in the eastern tropical Pacific is the well-known ENSO mode, a phenomenon that occurs mostly in the winter months of the northern Hemisphere as seen in NOAA CPC (2005). The winter ENSO pattern corresponds to the large eigenvalue λ₂ in winter months. Particularly December, November, and January have the three largest λ₂ values according to Figure 8. This is expected as higher eigenvalues imply larger variance and vice versa. The figure also shows the three months with the smallest λ₂ values are August, July, and September indicating there is little chance for an ENSO signal to occur in these months.

Figure 9

As expected from SST studies and textbooks, such as Shen and Somerville (2019) and Shen et al. (2017), Figure 9 shows January has an apparent ENSO signal at its surface. The January ENSO pattern vertically exists at the surface and continues to deeper layers. The pattern starts to wither at 150 meters and is entirely dispersed at and deeper than 200 meters (see Figure 9(g) and (i)).

At 100 meters (see Figure 9(e)), we see a large continuous blue region over the western tropical Pacific that corresponds to the cold anomalies of ENSO underneath the surface as shown in Shen et al. (2017) and similar to Roemmich and Gilson (2011). This cold anomaly stretches east-west over 7,000 km and north-south over 2,000 km, almost twice the area of the United States. This body of cold anomaly water continues to exist at 200 meters (see Figure 9(i)). Our EOF2 data shows that the body begins around 50 meters and ends around 250 meters. Also at the 100 meters (see Figure 9(e)), we see the warm tongue in the eastern tropical Pacific split by a small strip of cool waters just north of the equator. This might be part of the signal of equatorial upwelling.

Although the SST ENSO events rarely occur in August, our 3D EOF2 for August still shows a warm tongue pattern over the eastern tropical Pacific. Despite that, if we were to compare either column in Figure 9 it is evident the ENSO signal in August is less coherent compared to January. This may mean that the 3D ENSO dynamics continuously exist throughout the year and do not completely die out in the summer.

For the August EOF2, in the region south of the equator, the cold waters are more persistent than in January. To better understand what depth these cooler waters reach, a zonal cross-section from the surface to the bottom of the ocean (see Figure 10) is taken at 10°S. This cross-sectional map shows consistency in the behavior displayed by Figure 9(j).

Figure 10

The vertical map (Figure 10) cuts through the conventional SST ENSO pattern over the tropical Pacific. Figure 10(c) shows that the eastern tropical Pacific’s warm water body is limited to the top 150 meters, but the western tropical Pacific’s cold water body approximately spans 50 to 250 meters. In January there is also a cold-warm-cold pattern appearing in the Indian Ocean located within the top 250 meters as seen in Figure 10(c). However, this Indian ocean signal is not obvious for the August EOF2.

To show a more robust signal of ENSO, we have plotted the meridional mean from latitude 10°S to 10°N in Figure 11. Figure 11(c) shows a clear blue-red separation marked by the slanted yellow strip. The slope of the slanted yellow line is approximately –2 meters per longitude. This helps quantify the 3D structure of the El Niño dynamics, such as the thickness of the eastern tropical Pacific’s warm water to the surface, and the thickness of the western tropical Pacific’s cold water underneath the surface. January ENSO’s warm pattern ends at approximately 150 meters while August ENSO’s warm pattern ends at around 120 meters (see Figure 11(c) and (d)).

Figure 11

4.4 PCs, trend and de-trended

This sub-section describes the temporal patterns of the 3D EOF-PC analysis. Principal components (PCs) are the eigenvectors of the temporal covariance matrix and are also the column vectors of the orthogonal V matrix. PCs are on a time scale and show temporal patterns.

PC1 shown in Figures 12(a) and (b) show an upward trend of temperature anomalies. Then what would be the cause of this upward trend? Does the association of EOF1 with equatorial upwelling imply that the upwelling contains some heat from global warming? If it does contains some heat from global warming, how much? The detailed mechanisms of heat transfer and distribution are still to be investigated, even if the association between upwelling and global warming is true. Some literature has already appeared to address these questions. For example, Song and Colberg (2011) explored the existence of warming in the deep ocean; Wang et al. (2022) studied some specifics of global warming on equatorial upwelling. Studies like these may serve as the basis for further quantitative investigations on whether or how upwelling and global warming are correlated.

Figure 12

PC2 of January is supposed to be a time indicator of ENSO occurrence. This interpretation works well for the SST analysis on ENSO. When regarding 3D ocean, we may need to consider the time delay of heat transfer in the top 1500 meters. The El Niño pattern in the August EOF2 has already implied this time delay (see Figure 9). It is unknown how long it takes for the heat accumulation in the eastern tropical Pacific from the surface to 150 meters and for the heat loss in the western tropical Pacific in the deep layers from about 50 to 250 meters. Figure 12(c) shows that the minimum of January PC2 is about –0.32 at January 1990. According to NOAA CPC (2001) Oceanic Niño Index (ONI) based on SST, the northern hemispheric winter of 1989–1990 was a neutral year, while 1988–1989 winter was a La Niña year, 1990–1991 winter was neutral too, and 1991–1992 winter was an El Niño year.

The second minimum value of January PC2 (shown in Figure 12(c)) is about 0.30 occurred on January 1975, which was La Niña month. The maximum January PC2 is about 0.31 occurred in January 2003. The northern hemispheric winter of 2002–2003 had an El Niño, while the previous winter was neutral. The above analysis of the maxima and minima implies that PC2 from the 3D covariance matrix, although associated with El Niño, is not a direct indicator of the surface El Niño time. Further analysis may be needed to interpret PC2 and its association with El Niño.

A similar analysis may be done for the maxima and minima of August PC2. Figure 12(d) shows that the global maximum of the August PC2 is about 0.25 occurred in August 1979, which was a neutral month. The global minimum of the August PC2 is about –0.32 occurred in August 2003, the same year, i.e., 2003, as the global maximum of the January PC2. As discussed earlier, the northern hemispheric winter of 2002–2003 was in El Niño condition according to SST, but the summer of 2003 was in neutral condition.

The January PC3 shown in Figure 12(e) has a V-shape, going down from the 1950s to the 1970s, and then going up from the 1970s to the 1990s. This shape agrees with that of the Atlantic Multi-decadal Oscillation (AMO) index time series, and also that of the spatial average time series of the northern hemispheric average temperature, cumulative Southern Oscillation Index (SOI), and the cumulative weighted SOI1 (CSOI1) (see Figs. 4.2 and 4.3 of Shen and Somerville (2019)). We limit ourselves in claiming similar temporal patterns among these parameters. Whether the similarity is supported by climate dynamics may be worth further investigation.

August’s PC3 shown in Figure 12(f) has its global maximum in 1989 and several local minima in the 1960s and 1970s. This temporal pattern appears to agree with that of the North Atlantic Oscillation (NAO) (Hurrell (1995)). Again, we present here only the similarity of the temporal patterns. The climate dynamic support for this similarity needs further investigation.

The trend in PC1 motivated us to explore the EOF-PC analysis for the de-trended anomalies. This was done by subtracting the predicted values of the simple linear regression from the anomalies. The same volume weight is applied to the de-trended anomalies as described in the methodology section. Finally, the temporal covariance matrix can be calculated from the weighted de-trended anomalies to obtain EOFs, PCs, and eigenvalues.

By removing the trend we hope to gain insight as to why there is such a large jump of λ₁ between December and January shown in Figure 3. Figure 13 displays the resulting first eigenvalues ${\lambda }_1^{(d)}$ of the temporal covariance matrix from the de-trended anomalies for each month from January to December.

Figure 13

The shape of ${\lambda }_1^{(d)}$ for the de-trended anomalies is almost the same as that of λ₂ before the de-trending shown in Figure 8, although the values of ${\lambda }_1^{(d)}$ and λ₂ are not exactly the same. The large gap between December to January in λ₁ is gone in ${\lambda }_1^{(d)}$ (see Figures 3 and 13 respectively), implying the trend may have caused the gap. On that account the temporal trend in the 12-month time period from January to December is reflected in λ₁, corresponding to the equatorial upwelling mode.

In removing the linear trend the ENSO mode becomes the first mode. Comparing the PCs before de-trending the anomalies (Figure 12) and after de-trending the anomalies (Figure 14), it is clear that there is a mode shift: De-trended PC1 corresponds to PC2 before de-trend, and de-trended PC2 corresponds to PC3 before de-trend. This mode shift due to de-trending often occurs in EOF analysis (see Chapters 9 and 10 of Shen and Somerville (2019)).

Figure 14

4.5 Compare the 3D EOF1 to the 2D EOF1

As a way to compare results with other publications and reinforce the value of computing 3D EOFs, we computed 2D EOFs. It is important to note that 2D EOFs are orthonormal on each layer, in contrast, 3D EOFs are orthonormal only for the entire ocean domain of all layers, and not orthonormal on each layer. Figure 15 shows the 2D EOF1 for the top layer based on the anomalies, not the de-trended anomalies, from the NASA JPL OGCM.

Figure 15

To reiterate the 2D EOF calculation was only applied to the surface layer so ocean dynamics interconnection with other layers is not considered. Consequently, ENSO, not the equatorial upwelling, appears in the first mode of 2D EOFs, similar to Shen et al. (2017), Meyers et al. (1999), and Mu et al. (2018). The spatial pattern of Figure 15’s 2D EOF1 is very similar to that of the 3D EOF2 at the top layer shown in Figure 9(a). However, the data for Figure 9(a) is not normalized to one and is not orthogonal to the 3D EOF2 at the same layer.