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

Existing oceanic studies on either data reconstruction or dynamics often used 2-dimensional empirical orthogonal functions (EOF) for sea surface temperature (SST) and for deep layers. However, large-scale oceanic dynamics, such as equatorial ocean upwelling and arctic ocean ventilation, imply the existence of strong covariance among the temperatures and other parameters between different layers. These ocean dynamics are not best represented in the isolated 2-dimensional layer-by-layer calculations, while the 3-dimensional EOFs have a clear advantage. The purpose of this paper is to demonstrate 3D EOF calculations based on the NASA Jet Propulsion Laboratory (JPL) ocean general circulation model (OGCM) from surface to 5,500 meters depth, with 33 depth layers, 1-degree latitude and longitude spatial resolution, and monthly temporal resolution. We also present visualizations of the 3D EOFs and make physical interpretations of the first two EOFs. Our 3D EOF results demonstrate that (i) the 3D spatial pattern of equatorial ocean upwelling is mainly reflected in the first EOF mode and has its most variabilities within the depth layer between 100 and 400 meters, (ii) the 3D El Niño Southern Oscillation (ENSO) dynamic pattern is mainly reflected in the second EOF mode and is mostly confined from surface to the depth of 150 meters, and (iii) the lead eigenvalue from the 3D EOF calculation appears to contain some signal of oceanic warming. Additionally, our method of weighted 3D EOF computation and our 3D visualization Python code may be useful tools for both climate professionals and students.


INTRODUCTION
Empirical orthogonal functions (EOFs), or spatial modes, are often used to analyze general climate model (GCM) output to explore spatial patterns of climate variability or for reconstructing historical climate data (Zhang and Moore (2015), Houser et al (2020)).EOFs show which spatial patterns explain the most variance in the data, and can consequently give insight into how our data behaves.Although GCMs and climate dynamics are 3-dimensional in space and 1-dimensional in time, EOF analyses have often been performed in 2D layers in space (Shen et al. (2017), Meyers et al. (1999), Mu et al. (2018)).The primary purposes of our paper are twofold: (a) to show a computational algorithm of 3D EOFs, and (b)to interpret the climate dynamics of the first two EOFs.We do so by using the output data of a NASA JPL ocean general circulation model (OGCM).
Numerically, EOFs can be computed by finding the eigenvectors of the spatial covariance matrix from the weighted anomalies.This spatial covariance method often has convenient physical interpretations for leading modes.Most EOF calculations in recent years use singular value decomposition (SVD) with packages available in different computer languages, e.g., Python, R, and MatLab, when the size of the SVD matrix is not too big, say less than half of your computer memory (Shen and Somerville (2019)).Our 1-degree spatial resolution dataset is not very big and SVD computing can be done.However, we also use temporal covariance in this paper which can address the issue of memory limitations for future datasets of higher resolutions.This memory limitation is perhaps one of many reasons that EOFs were used to be computed layer-by-layer.Our method details in this paper are presented in the way of covariance.
Papers such as Shen et al. (2017), Meyers et al. (1999), and Mu et al. (2018) calculated 2D EOFs, and indicated that the strongest signal was El Niño Southern Oscillation (ENSO).ENSO is an ocean-atmosphere climate phenomenon characterized by the warming of sea surface temperatures in the central and/or east-central equatorial Pacific (NOAA CPC (2005)).The pattern is often described as a warm tongue in the eastern pacific.Despite its aforementioned clear spatial pattern in sea surface temperature (SST), ENSO is a phenomenon involving not only the surface but also the deep ocean, up to 150 meters, to be discussed later in this paper.Therefore, 2D EOFs lack the ability to describe the interactions between surface and deeper ocean layers.The need to consider such significant deep ocean dynamics motivated the idea of computing EOFs for all layers at once.The resulting 3D EOFs show that ENSO accounts for a smaller fraction of variability than equatorial upwelling.This validates the significance of calculating 3D EOFs.
Our paper will demonstrate 3D EOF calculations based on the NASA JPL ocean general circulation model (OGCM) from surface to 5,500 meters depth, with 32 depth layers, 1-degree latitude and longitude spatial resolution, and monthly temporal resolution.Our mode visualization and physical interpretation will focus on the first two EOFs.Our 3D EOF results demonstrate that (i) the 3D spatial pattern of equatorial ocean upwelling is mainly reflected in the first EOF mode and has its most variabilities within the depth layer between 100 and 400 meters, (ii) the 3D El Niño Southern Oscillation (ENSO) dynamic pattern is mainly reflected in the second EOF mode and is mostly confined from surface to the depth of 150 meters, and (iii) the lead eigenvalue from the 3D EOF calculation appears to contain some signal of oceanic warming.
Relevant literature already exists on (i) multiway table data analysis theory, e.g, Structuration des Tableaux A Trois Indices de la Statistique (STATIS) (Acar and Yener (2009), Vallejo-Arboleda et al. (2007), Sabatier and Traissac (1994)), and (ii) computational algorithms for big data EOFs (Gittens et al. (2016)).A special case of multiway table is a 3-way table that includes locations, parameters, and time (Stanimirova et al. (2004)).In our case, we only look at only one parameter: temperature; our locations are three-dimensional; and our time is indexed for monthly data from the NASA JPL OGCM.
In principle, multiway data analysis algorithms such as STATIS should be able to deliver the results of PCs and EOFs when the multiway data table is properly nondimensionalized.In practice, the problem is much more complex, such as the large size of the multiway data table, and the physical meanings of the PCs and EOFs.For example Gittens et al. (2016) contributed to the development of innovative computational algorithms to handle multi-terabytes of data for 3D EOF computing.Their paper considered the cases: (a) Climate Forecast System Reanalysis (CFSR) data with a 1/2-degree, 41 layers, and 6-hour resolution; (b) Community Atmosphere Model 5.0 (CAM5) data with 1/4-degree, 30 layers, and 6-hour resolution.Their CFSR dataset has 2.2 TB and their CAM5 dataset has 16 TB.They concluded that most of the run time was spent on In/Out process.Their EOFs from their multi-hour data by default have no clear physical interpretations.Accordingly, their paper included neither EOF figures nor physics interpretations.Our paper shows both the successful application of our 3D EOF algorithm that utilizes a temporal covariance matrix designed for future applications on multi-terabytes of data, and physics interpretations of the 3D EOFs from the NASA JPL OGCM.Therefore, our study has a modest purpose to show the feasibility of (i) computing 3D EOFs from a climate model output, and (ii) making physics interpretations from the point of view of ocean dynamics.
The remainder of this paper is organized as follows: Section 2 describes the data set, Section 3 describes the proposed method, section 4 includes results of performing the method on OGCM data, and section 5 contains the conclusion.Our method of weighted 3D EOF Lafarga et al. Tellus A: Dynamic Meteorology and Oceanography DOI: 10.16993/tellusa.3223computation and our 3D graphics Python code may be useful tools for both climate professionals and students.The comprehensive code can be accessed at https:// github.com/dlafarga/calc_3D_EOFs.

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.
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.

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 1950 ) and first layer (l 1 ) 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.
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 spacetime 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, 2,073,600 64,800 32 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).

METHODOLOGY
There are two methods described in this section: SVD and temporal covariance.We will first describe the SVD method to better understand the relationship between each method.We start by calculating our climatology.We use the entire data history from 1950 to 2003 as our climatology period.Let us use the following simple spacetime data matrix notation [T it ] to explain our calculation of climatology and anomalies: where i is a grid box simplified from the previous latitude-longitude notation, t is time, N is the total number of grid boxes, and Y is the total number of years.Specifically, N = 2,073,600, and Y = 54.This X can be the space-time data matrix for a particular month, January, February, or another month.
The climatology is computed as the time means of the data matrix, i.e., the mean for each row: We compute the climatology for each month: January, February, ……, December.Thus, 12 sets of climatology are computed.Figure 2 shows the climatology of January and that of August.
Our anomalies are computed as follows: These anomalies are multiplied by their spatial weights to, as mentioned previously, better reflect the geometric sizes of their corresponding grid boxes.A 3D grid box in the equator region is much larger than one in the polar region.In 2D the weight is related to the area of a grid box (see equation 4).Naturally, in our 3D case the weight for a 3D grid box i should be related to the volume of the grid box: where Δd i is the thickness of the layer for the grid box i and ϕ i is the latitude of the grid box i in radians.The square root is needed due to the symmetry of the formulation of the eigenvalue problem of a covariance matrix (Shen et al.(2017)).
Often, anomalies are standardized by dividing by the standard deviation.However, in the deeper ocean layers, temperature has a very small standard deviation.Dividing by the very small standard deviation would only amplify noise and cause large errors for the EOFs in the top layers.To avoid such errors, we simply use non-standardized volume-weighted anomalies.Up to now, the NaN missing data are still in the anomaly data matrix.We remove the rows of NaN and produce a volume-weighted space-time anomaly matrix without missing data.We denote this weighted anomaly matrix without missing data by where N′ is the number of grid boxes with data.Then, we can feed this space-time data matrix to SVD to generate EOFs to demonstrate spatial patterns and principal components (PCs) to show temporal patterns.The SVD decomposition of B is as follows: where the superscript T indicates the transpose of a matrix.In this SVD formula, U is an N′ × Y (N′ > Y) spatial matrix.Each column vector of U is an EOF.Matrix V has dimension Y × Y and is the temporal matrix.Each column is a PC.Matrix D is a diagonal matrix, whose elements d k (k = 1, 2, …, Y) correspond to the standard deviations of the EOFs and PCs.The next sections will explain these EOF and PC patterns in U and V, as well as the standard deviations in D.
Although SVD is a convenient way to compute EOFs and PCs, the covariance matrix method is often helpful in interpreting the meaning of EOFs and PCs from the point of view of variance and covariance in climate dynamics.Further, this method was traditionally used in the EOF calculations before using SVD in recent years.The traditional method often uses the eigenvalue problem for a spatial covariance matrix to compute EOFs.However, in our case, the temporal dimension Y is much smaller than the spatial dimension N′, hence the temporal covariance matrix is much smaller than the spatial covariance matrix.It is preferred to solve the eigenvalue problem for the temporal covariance matrix which is much smaller than the spatial covariance matrix.The resulting eigenvectors of the temporal covariance matrix are PCs: The EOFs can be obtained by the projection of the weighted anomaly matrix B to the PCs: where Σ is defined as a vector whose orientation is not changed after the action of [ ] tt′ Σ .Namely, the matrix action vector [ ]  is parallel to the original vector.Thus, there exists a scalar, denoted by λ so that one vector is the scalar multiplication of the other, i.e.,

[ ]
This eigenvalue problem has Y solutions whose Y eigenvectors ( , 1,2, , ) are the same PCs computed from SVD (see the rigorous mathematical proof in Shen and Somerville ( 2019)).The eigenvalues (λ Eigenvalues represent the amount of variance explained by each respective eigenvector (National Center for Atmospheric Research Staff (Eds) ( 2013)) and are sorted from the greatest variance to the smallest.Variance is a measure of how far a set of numbers is spread out from their average value (Shen and Somerville (2019)).Thus, the first eigenmode will show us the most significant pattern, as it has the most variance in climate dynamics.The last modes are usually less accurate and less important as they contain much noise and smallscale variations (Liang et al. (2012)).The details of the EOF and PC derivations can be found in Shen et al. (2017).
We have computed the EOFs, PCs, standard deviations, and variances for each mode using both the methods of SVD and temporal covariance.We have compared their results as a cross-check and explored the physical interpretation of their meaning.These results are described in the next section.

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 Lafarga et al. Tellus A: Dynamic Meteorology and Oceanography DOI: 10.16993/tellusa.3223and 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.
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 λ 1 appears.
Figure 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., 1 100.
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 1 = 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 λ 4 and λ 5 are close to each other and have a very small difference p 4 -p 5 .EOF errors are large because they are proportional to the inverse of the difference of two neighborhood eigenvalues, i.e., 1/(λ 4λ 5 ), 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 2 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.

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 λ 1 .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.
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)).
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-coolwarm 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).
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), andShen 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  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.
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).

EOF2 AND THE SPATIAL PATTERN OF EL NIÑO
Figure 8 shows the eigenvalue λ 2 of the temporal covariance matrix for every month.In contrast to λ 1 in Figure 3, there is no longer a large jump from December to January in λ 2 .As previously mentioned the smooth transition is due to an annual cycle.
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 λ 2 in winter months.Particularly December, November, and January have the three largest λ 2 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 λ 2 values are August, July, and September indicating there is little chance for an ENSO signal to occur in these months.
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 northsouth 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).
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.

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.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, 1990Niña year, -1991Niña year, winter was neutral too, and 1991Niña year, -1992 winter was an El Niño year.
The second minimum value of January PC2 (shown in Figure 12(c 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 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 λ 1 between December and January shown in Figure 3. Figure 13    The shape of ( ) 1 d λ for the de-trended anomalies is almost the same as that of λ 2 before the de-trending shown in Figure 8, although the values of ( ) 1 d λ and λ 2 are not exactly the same.The large gap between December to January in λ 1 is gone in ( ) 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 λ 1 , 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)).

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.
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), andMu 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.

CONCLUSIONS AND DISCUSSION
The main purpose of this paper is to show the method of computing 3-dimensional EOFs and interpret the physical meanings of the first two modes.To compute these EOFs, we used the temporal covariance of the volume-weighted anomalies based on the NASA JPL OGCM's monthly output of 1.0 latitude-longitude degree, 33 layers from surface to 5,500 meters depth, and a time span from 1950 to 2003.The eigenvectors of the temporal covariance are PCs.The normalized projections of the volume-weighted anomalies for the PCs form the geometric EOFs.Removing the volume factors we obtain the physical EOFs, which are displayed in our EOF figures.
Our method to compute the January PCs starts from the formulation of the N × Y space-time data matrix with the first dimension being the total number of grid boxes and the latter being the total number of years.Anomalies are multiplied by the volume weights of each depth, resulting in the weighted anomalies.After removing land/NaN data, the new weighted anomaly matrix transpose multiplied by itself forms the temporal covariance matrix.
We also computed the EOFs and PCs using SVD for the volume-weighted anomalies.Although SVD is capable of handling our data at the 1.0-degree resolution, it cannot handle the higher resolution at 1/4 degree due to computer memory resource limitations.Thus, our temporal covariance approach is useful and will be used later on to deal with 1/4-degree data or higher resolutions.Our application of SVD in this paper is for the purpose of verifying our temporal covariance approach.
NASA JPLs OGCM was used before in Shen et al. (2017) where EOFs were computed layer by layer, i.e., the twodimensional EOF, or 2D EOF.The paper shows the ENSO signal in the first one or two EOFs, and it did not reveal the equatorial upwelling mode.This was due to the nature of a 2D calculation, which does not take into account the deep ocean dynamics interconnected among different layers and instead treats each layer in an isolated fashion.Our 3D EOFs show covariances across depth layers.Consequently, equatorial upwelling appears as 3D EOF1 and explains a third of the data's variance.In contrast, the ENSO mode, or 3D EOF2, only explains 10% of the variance.
By observing the covariance among the temperatures of different layers, we were not only able to identify different ocean dynamics than the previous 2D calculations, but also we were able to pinpoint how deep each characteristic spans.Equatorial upwelling can be seen from 100 to 400 meters.ENSO is observed mostly in the upper layers until 150 meters, and at 100 meters we can observe a large cool water patch in the western tropical Pacific ocean.
Although through the 3D EOFs and PCs, we were able to identify some meaningful ocean dynamic patterns both in the deep ocean and surface, there is still much to explore in the resulting dynamics and the effect of de-trending.Further exploration of the ocean dynamics may be made to show how the 3D EOFs are related to climate oscillation patterns, such as AMO and NAO.Our EOFs can serve not only as a means to understand dynamics but also help in reconstructing historical climate data (Shen et al. (2017)).The 3D nature of the EOFs allows the covariance properties between layers to be used in the data reconstruction.The covariance between points of different depths are expected to add skills to reconstruction.For example, the 3D insitu observational data can likely support more 3D EOF modes than the 2D EOF modes in the sub-surface of the ocean where in-situ observations are scarce.Can we extend the reliable reconstruction of historical data to a longer history compared with the conventional 2D reconstruction?Can the historical in situ data support a set of 3D EOFs with a finer spatial resolution?More questions like these can be asked about the 3D EOFs with regard to their applications, computations, climate dynamics, and mathematical theory.
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 2D Jan X, looks like the following:

Figure 1
Figure 1 Ocean water temperature from the NASA JPL OGCM output for January 1950.A movable version of the figure can be found at https://ogcm-3d-visualization.herokuapp.com/where the figure can be rotated and resized.

10 )Figure 2
Figure 2 Climatology for January and August computed using the equation above from the NASA JPL OGCM data.Movable figures can be viewed at the website https://climatology-3d-vis.herokuapp.com/apps/clim_Jan_plot.

Figure 3
Figure 3 (a) The first eigenvalue λ 1 of the temporal covariance matrix for the weighted temperature anomalies for each month from January to December.Units of the first eigenvalue are the same as the covariance matrix in equation 14.(b) The scree plot for January based on data p k and q k .(c) The scree plot for August is based on data p k and q k .

Figure 4
Figure4EOF1 for January and August showing the equatorial upwelling.Depths are chosen to show that the strong upwelling occurs at 150 meters and that the upwelling stops at a deeper level around 400 meters.

Figure 5
Figure 5 A north-south cross-section of the ocean taken at 160° E to show at what layers upwelling occurs.Panels (a) and (b) are the entire scope of all ocean layers.Panels (c) and (d) show the top layers up to 400 meters as this will better show where upwelling starts and stops.
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.

Figure 6
Figure 6 NASA JPL OGCM temperature for January and August at 200 meters showing equatorial upwelling.Years are chosen to reflect time with strong upwelling.

Figure 7
Figure 7 Cross-sectional map based on the zonal mean from 0 to 360° longitude degrees to show ocean ventilation at the high latitude regions.

Figure 8
Figure 8The second eigenvalue λ 2 of the temporal covariance matrix for the weighted temperature anomalies.Units of the first eigenvalue are the same as the covariance matrix in equation 14.January has the largest value, while August has the smallest value.

Figure 9
Figure 9 EOF2 for January and August showing the ENSO pattern.Depths are chosen to show the ENSO regions in different depths from surface to 100 meters.The ENSO pattern dies out approximately at levels deeper than 150 meters.
PC1 shown in Figure12(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

Figure 10
Figure 10 Cross section at latitude of 10°S that shows the separation of the warm and cold anomaly regions in the Pacific.

Figure 11 A
Figure 11 A zonal cross-sectional map based on the average from latitude 10°N to 10°S.The average and figure are to better show a robust 3D ENSO pattern with not only the surface warm region but also a depth structure.

Figure 12
Figure 12 Principal components PC1, PC2, and PC3 of January and August from the weighted anomalies of the NASA JPL OGCM output from 1950 to 2003.
) 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.
displays the resulting first eigenvalues ( ) 1 d λ of the temporal covariance matrix from the de-trended anomalies for each month from January to December.

Figure 13
Figure 13 The first eigenvalue ( ) 1 d λ for the each month computed from the de-trended anomalies.

Figure 14
Figure14Principal components of January and August based on the de-trended anomalies.

Figure 15
Figure 15 2D EOF1 for the top layer based on the area-weighted anomalies from the NASA JPL OGCM.ENSO signal is shown in this EOF1 for the surface layer.The 2D EOFs are computed layer by layer.

Table 1
Depths of 32 layers of the NASA JPL OGCM.