The Use of MTM-SVD Technique to Explore the Joint Spatiotemporal Modes of Wind and Sea Surface Variability in the North Indian Ocean during 1993–2005

1 Department of Mathematics, Faculty of Science, Burapha University, Chonburi 20131, Thailand 2 Center of Excellence in Mathematics, Burapha University, Chonburi 20131, Thailand 3 Department of Civil, Environmental, and Architectural Engineering, University of Colorado, Boulder, CO 80309, USA 4 Co-Operative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO 80309, USA 5 Department of Aerospace Engineering Sciences, University of Colorado, Boulder, CO 80309, USA


Introduction
North Indian Ocean is the least explored among the oceans in the world. It is forced by seasonally reversing monsoons, which play an important role in its variability. In particular, the variability of monsoonal winds in the North Indian Ocean greatly influence the variability of sea surface height (SSH) and sea surface temperature (SST) which have an impact on the regional rainfall and consequently, the socioeconomic well being of the large population in the surrounding regions. Thus, a better understanding of the nature of this variability is critical for improved resource planning and management.
The ENSO-monsoon relationship during 1997-1998 created climate and oceanic anomalies [1,2]. SST anomaly resulting from ENSO condition yielded SSH anomalies along the South Asia coast. Sea level increased and retarded the outflow in estuaries resulting in floods in the deltaic area especially in Bangladesh and India [3,4].
The SST anomaly is also an important key to climate change. One strong SST variation is an out-of-phase SST in the tropical Indian Ocean, the so-called Indian Ocean Dipole [1,2,5,6] or the Indian Ocean zonal mode (IOZM). This is a nonstanding mode similar to the meridional mode in the tropical Atlantic Ocean [7]. Warm SST anomaly in the Indian Ocean, usually of the order of 0.5 • C, may lead to a strong monsoon [8][9][10]. The monsoon affects not only the South Asia area but also East Africa. Some studies show that the SST variation in the Indian Ocean is a major contribution to the East African precipitation variation [11][12][13]. The ENSO signal has a substantial modulating influence on the climate variability in the Indian Ocean basin including the IOZM [10,[14][15][16][17].
In the past few decades, several mathematical and statistical techniques have been developed and used to identify variability from signal data. Techniques such as Singular Value Decomposition (SVD) [18,19] and related orthogonal 2 International Journal of Oceanography multivariate spatiotemporal decompositions [19], Spectral Analysis [20], Principal Oscillation Patterns [21], Principal Component Analysis [22,23], Wavelet Spectrum Analysis [24] as well as Multiple-Taper Method-Singular Value Decomposition (MTM-SVD) [25] have been applied to the analysis of atmospheric and oceanic data. Some methods are limited to identifying variability in either frequency domain or spatial domain. Some methods can only identify standing modes with no information on dynamics and some methods can only identify localized variation.
Unlike other univariate or multivariate methods, MTM-SVD is a powerful multivariate tool to simultaneously detect an oscillatory signal in both spatial and temporal data that can be used to identify frequencies, where an unusual concentration of narrowband variance occurs. It can also be used to reconstruct the time history and spatial pattern associated with a frequency of interest [26]. Mann and Park [27] applied the MTM-SVD to investigate decadal climate variation in the Northern Hemisphere. Tourre et al. [28] identified joint variability between SST and sea level pressure with MTM-SVD. Rajagopalan et al. [29] and Mann and Park [26] used the MTM-SVD to show how reliable the method was, by testing through a synthetic data set with colored noise added, and then applied the method to climate data.
In this study the joint variability and dynamics of wind and two surface parameters, SSH and SST, will be investigated using MTM-SVD. The paper is organized as follows. Data used in this study is first discussed, followed by a brief description of the MTM-SVD on joint variability. Results of MTM-SVD analysis are presented next. Discussions of the results conclude the paper.

Data and Method
2.1. Numerical Ocean Model. Data used in this study is obtained from a data-assimilating numerical ocean model applied to the North Indian Ocean in a hindcast mode for the period 1993-2005 [30]. The ocean model used is the University of Colorado version of the Princeton Ocean Model (CUPOM), using topographically conformal coordinate in the vertical and orthogonal curvilinear coordinates in the horizontal. Details of the basic features of CUPOM can be found in [31,32].
The North Indian Ocean hindcast model has 1/4 • resolution in the horizontal and 38 sigma levels in the vertical, with the levels closely spaced in the upper 300 m. The model is forced by 6-hourly ECMWF winds with Smith [33] formulation for the drag coefficient C D to convert the ECMWF 10-meter wind speed to wind stress at the surface. It assimilates altimetric sea surface height anomalies and weekly composite Multichannel Sea Surface Temperature using a simple optimal interpolation-based assimilation technique. The model and the assimilation methodology based on conversion of SSH anomalies into pseudotemperature anomalies for adjusting the model temperature field via optimal interpolation, have been described in detail in [34,35]. Details of CUPOM applied to the North Indian Ocean along with validation of model's results can be found in [30,36]. The SSH is calculated explicitly using the splitmode technique.
The hindcast produced a continuous space-time data of wind and all the surface parameters (i.e., SST, SSH) for the 1993-2005 period, in addition to 3-D fields of water mass properties and currents. Weekly values are computed from the model data for use in this study. To keep the data size reasonable and maintain all the spatial information, data at every 1.

Frequency Domain Multiple-Taper Method-Singular
Value Decomposition (MTM-SVD) Approach. Robust diagnosis of the key low-frequency modes of large-scale climate entails capturing the coherent space-time variations across multiple climate state variables. Traditional time-domain decomposition approaches for univariate and multivariate data provide useful details on the broad-scale patterns of variability. However, these approaches lack the ability to isolate narrowband frequency domain structure [25,27]. Detailed methodology development and examples of the MTM-SVD methodology can be found in [25,27,37,38]. The method relies on the assumption that climate and/or oceanic modes are narrowband and evolve in a noise background that varies smoothly across the frequencies. Subsequently, spectral domain equivalents of each grid point are computed based on the multiple-taper spectral analysis [37,39].
In this study, the MTM-SVD approach is applied to examine the individual and joint spatiotemporal modes of variability of wind and sea surface parameters in the North Indian Ocean. The spatial time series of each parameter is standardized by removing the long-time mean at each grid point and normalized by dividing the long-term standard deviations by weekly (or monthly) basis so that the variability has a unit variance and also the weekly (seasonal) cycle is removed. This normalization eliminates the disparity in the units between the variables. The details of the technique can be found in the aforementioned references but below we provide a brief description of the method abstracted from the references for the benefit of the readers. n /σ (m) , where prime notation represents the anomalies and σ (m) is the standard deviation at the mth site. Note that the location of each variable needs not be at the same site but both spatiotemporal time series must have the same length in temporal space. The multiple-taper transform is next applied to the time series and the "eigenspectra" are obtained at each frequency: where Δt is the sampling interval and {w (k) n } N n=1 is the kth member in an orthogonal sequence of K data tapers (also called Slepian tapers), k = 1, . . . , K. The set of K tapered eigenspectra have energy peaks within a narrow frequency bandwidth of f ± p f R , where f is a given frequency and f R = 1/(NΔt) is the Rayleigh frequency (the minimum resolvable frequency range for the time series). Note that the number of tapers K represents a compromise between the variance and frequency resolution of the Fourier transforms [27]. Also, note that the Slepian tapers can be considered as a sequence of weighting function. An example of the first three Slepain tapers can be found in [26].
With the set of eigenspectra, one can form an (M +L)×K matrix: where each row is calculated from a different grid point time series. Each column uses a different data taper from (1) and w i represents specific weighting at each grid point. Note that each row is computed from a different series (grid point), and each column corresponds to using a different taper. Subsequently, a complex singular value decomposition (CSVD) is performed through where u k is complex spatial pattern and v k is spectral domain containing information of both variables (i.e., the first M components contain information on the first variable and the next L − M components contain information on the second variable). These eigenvectors can be inverted to obtain the smoothly varying envelope of the kth mode of variability at frequency f [27]. The localized fractional variance (LFV) provides a measure of the distribution of variance by frequency, and above a select confidence level threshold (e.g., 90%, 95%), represents a dominant narrow band mode. The confidence levels are computed based on the local white noise assumption and are constant outside the secular band. Mann and Park [27] describe a bootstrap method used to obtain the confidence bands for this study.
In general, the computed principal eigenspectrum (described above) yields a number of narrowband peaks. The MTM-SVD technique has been effectively applied to the analysis of global SSTs and SLPs [25,27], identification of dominant modes of variability in the Atlantic basin [28], and also for forecasting [29]. The LFV spectrum is used to identify significant frequencies, and temporal and spatial reconstructions are carried out to understand the joint variability of the climate fields in the North Indian Ocean. In this study, the choice of bandwidth parameter p = 2 and K = 3 tapers is used. With these choices, they provide enough spectral degrees of freedom for signal-noise decomposition and allow reasonably good frequency resolution as well as stability of spectral estimates according to [25,40]. The monthly data allows for a maximum frequency of 6 cycle yr −1 (2-month period) to be resolved efficiently.
In the next section, the dominant frequencies of each individual surface parameters SSH and SST are isolated. Next the individual dynamics of SST and SSH will be discussed followed by the analysis of spatiotemporal variability jointly between a pair of zonal and meridional wind stresses, and SSH and SST which provide insight into possible dynamical processes governing such signals.

Results
The MTM-SVD method is applied to the fields individually and also jointly. In this section results of the individual analysis of SSH and SST variability are first presented followed by the joint analysis of winds and SSH, and winds and SST.

Sea Surface Height.
The LFV spectrum of SSH based on the analysis of monthly data is shown in Figure 2 ( f ∼ 4.0 cycle yr −1 ). All these peaks appear to breach the 99% confidence level. The LFV spectrum from the weekly data (not shown) is similar to the one from monthly data.

Joint Wind and SSH Variability.
North Indian Ocean, as mentioned earlier, is most influenced by the Indian summer monsoon, thus, the variability of SSH and SST is likely to be tied firmly to the monsoonal winds. To identify the joint variability, the joint analysis of the winds (zonal, WSX and meridional, WSY winds separately) with SSH is performed. Figure 3 shows the LFV of this joint analysis. The LFV of joint variability of both pairs (WSX-SSH and WSY-SSH) yield significant peaks (over 99% confidence limit) at interannual, annual, and semiannual periods and at seasonal and intraseasonal (at 90% confidence level). The significant peaks observed in this joint analysis are the same as the frequencies identified in the individual analysis of SSH. Here too, the LFV spectra from the weekly and monthly data are similar.
To understand the joint variability in space and their evolution at these dominant frequencies, the spatial reconstructions for the two fields at the frequencies identified above are performed. Spatial patterns of the zonal and meridional wind fields and the corresponding SSH at the annual cycle ( f = 1.0 cycle yr −1 ) frequency, from the respective joint analyses are shown in Figure 4. The vector lengths in this figure (and subsequent spatial reconstructions) are indicative of the magnitude of the signal, and the phases (i.e., directions) are with respect to a reference location (grid 529th at 7.25 • N, 76.75 • E)-the vector at the reference location is always horizontal pointing right (i.e., at the 3 o'clock position). The direction of vectors at other locations corresponds to temporal difference (time lead/lag) relative to the reference vector. In the spatial pattern, clockwise vector rotation represents negative relative lag (or "lead"), and counter-clockwise rotation represents positive lag (or "lag"). A complete rotation represents a periodicity of the mode; for example, 1 year for f = 1.0 cycle yr −1 , and a grid point vector at 12 o'clock (90 • ) experiences peaks at 4 months lag relative to the reference grid point.

Annual Frequency.
The zonal winds (Figure 4(a)) exhibit an out-of-phase relationship between the North Indian Ocean region and the Southern Indian Ocean region-this is indicative of the annual cycle in that the winds in the northern hemisphere are always out of phase with the southern. The meridional winds (Figure 4   Ocean, the Bay of Bengal region, and the South China Seaall active regions of the Asian monsoon, and Indian monsoon in particular. In fact, these wind patterns are almost identical to the monsoon wind features [42][43][44]. The vectors are all in the same direction, indicative of the feature happening in the same period (i.e., the summer period), the slight changes in vector directions between the South China Sea region and North Indian Ocean shows the delay in the South China monsoon relative to the Indian. The spatial pattern of SSH variability at this mode (Figures 4(b) and 4(d)) from the joint analysis of the two wind components is similar. The vectors at the same area only differ in their magnitudes and phase due to the difference in reference vectors. The SSH signal is strong over most of the basin except in the southern hemisphere where the winds are quite a bit weaker as seen above. Variability is dominant in most areas of the North Indian Ocean except the area 50 • -80 • E and 5 • S-0 • N. From these plots, one can notice that the strong magnitude of SSH in the southwest region is not locally affected by the zonal wind since the zonal wind in that region is weak during that time. It could be remotely affected from strong easterly wind in the south leading by about 2 months.
The propagation of the SSH signal at the annual cycle frequency is apparent and can be explained as follows. Starting in the black-vector region in the southern Indian Ocean, it takes about 1-2 months (∼30 • -45 • ) to propagate westward into the southwest (green-vector) area. The large magnitude along the west coast of India propagates westward into the red-vector area in the southern Arabian Sea. Similarly, the signal in the southwest propagates northward along the Somali coast into red/blue-vector areas. Then the red area near the Socotra Island propagates eastward into the middle of Arabian Sea, while the blue area moves southward into the Equatorial region, which next propagates eastward along the equatorial waveguide into the Sumatra Island area and completes the cycle with a stronger magnitude along the east coast of Bay of Bengal. The propagation of the SSH signal lagged by a few months to the wind signal-which is consistent in that the SSH anomalies are driven in large part by the strong wind forcing in the basin [45][46][47].

Semiannual.
The spatial patterns of the wind and SSH fields at the semiannual frequency ( f = 2.0 cycle yr −1 ) are shown in Figure 5. At this period, the strongest signal in both the wind components is in the Arabian Sea regionwhich is the active Indian monsoon region. This is more so in the meridional component where the amplitudes are strongest in the Western Arabian Sea region and also a bit in the South China Sea. The main aspect at the annual and semiannual frequencies is similar in terms of the winds, in that the variability is strongest in the regions of monsoonal winds and the phase lags are consistent with spatial and temporal variability of monsoons in the region.
The SSH signal is stronger with the meridional winds ( Figure 5(d)). The large magnitude of SSH between 2 • -10 • N and 50 • -77 • E (referred as area A1 thereafter) is affected by strong WSY in the western Arabian Sea, which leads this by about a month or so. The strong magnitude of SSH in the southern region on the other hand is forced by zonal wind stress locally and is also affected by the strong magnitude of SSH in the area A1. at a reference grid point (grid 529th). The initial snapshots (at 0 • ) correspond to peak WSX anomalies in the east-central part of the North Indian Ocean and the corresponding SSH anomalies. Next snapshot reveals the evolution of WSX and SSH anomalies. The WSX anomalies are weakening in the next few months and reverse signs in about 6-7 months (at about 90 • ); though the SSH anomalies are also weakening but at a slower rate in about 7-9 months which agrees with SSH anomaly discussed in [48]. While the final snapshots correspond to the opposite conditions revealed about 15 months later. At the phase between 112.5 • and 135 • , large magnitude of SSH anomalies occur within 10 • S-10 • N across the ocean with one sign covering most of the western side of equatorial region and the other sign occurs next to the Sumatra Island. The anomalies in the eastern side intrude into the middle of this region. This is reminiscent of the "dipole" feature [1,2,5,6] identified in the wind fields.
Clarke and Liu [47] suggested that SSH variability is caused remotely by equatorial zonal wind. To analyze this, the temporal variability is examined by reconstructing the time components of these fields at locations (see Figure 7) from regions exhibiting strong magnitudes in Figure 6. The temporal reconstructions are shown in Figure 8. Relationship between WSX and SSH (b) at Sumatra is almost inphase, while the WSX in this region leads SSH (a) at the western equatorial region by 10-12 weeks. While the WSX in this region reveals 180 • out-of-phase to the SSH west of India (bottom). This suggests that the SSH variability at the eastern part of the equatorial region is directly affected by WSX. In the signal, there is a rather sudden change in WSX during the late 1997 that can be seen with an anomalous high in SSH at that time. This potentially could be driven by the strongest ENSO event in 1997-1998.

Wind and SST.
The LFV spectra of the joint analysis of WSX-SST and WSY-SST are shown in Figure 9. The spectra are very similar to that of joint wind and SSH analysis from before. Significant peaks are found at the annual, semiannual (over 99% confidence limit), and interannual bands (over 90% confidence limit).
The spatial reconstructions of the SST fields at the annual frequency periods (not shown) show large anomalies with   06  05  04  03  02  01  00  99  98  97  96  95   a phase difference of 4-5 months between the Northern Arabian Sea, Northwestern Bay of Bengal and South China Sea, and the Southwest region near the Tanzanian coast. The wind fields lead the SST variations in the Northern areas by a few months and that the strong SST variation in the southwest region of the domain is affected by both local WSX and WSY variation. At the semiannual ( f = 2.0 cycle yr −1 ) frequency, the SST anomalies are large almost everywhere with largest in magnitude found in the Arabian Sea and smallest found in the eastern side of the equatorial region. At the interannual frequency ( f = 0.3 cycle yr −1 ), the SST anomalies do not show much difference in magnitude throughout the Indian Ocean, which is consistent with [48], since the SST variations in the Indian Ocean during the ENSO periods are about 2 to 3 times smaller than those in the Pacific. From the spatial reconstruction the WSX anomalies in the equatorial region lead the SST anomalies in the Southwestern of the Sumatra by 3 to 5 months.

Conclusion
In this study, the MTM-SVD technique was applied to decimated weekly and monthly data of individual and joint fields of winds and surface variables. The MTM-SVD technique provides an attractive and complementary alternative to the traditional time domain analysis methods. Especially, its capability in isolating narrowband propagating features in the joint fields is unique unlike the traditional methods. This analysis was able to identify significant frequencies bands that are shared by a majority of spatial locations in the joint fields. It showed significant variability at semiannual, annual, and interannual frequencies in these fields individually and jointly. At the semiannual and annual frequency bands, the variability is largely linked to the seasonal reversing monsoons in the basin, which is a strong factor. The variability at the interannual frequency band is related to ENSO and the evolution of the SST spatial pattern resembles at time to a nonstanding so-called the Indian Ocean dipole. Furthermore, it appears that the winds are the main drivers of variability in both SSH and SST, which propagate around the basin. The spatial and temporal reconstructions offer prospects for long lead simulation and forecast of the climate in the North Indian Ocean which can be of great help to societies in the region that support a large part of the world's population.
One can extend the analysis of MTM-SVD technique to more than two fields but it requires more computer resources. However, a disadvantage of MTM-SVD technique is that the dominant variability at any frequency requires strong variance in many spatial locations. If there are only a few spatial locations with strong variability participating at that frequency, that area will not be dominant.