ABSTRACT Title of Document: BREEDING ANALYSIS OF GROWTH AND DECAY IN NONLINEAR WAVES AND DATA ASSIMILATION AND PREDICTABILITY IN THE MARTIAN ATMOSPHERE Yongjing Zhao, Doctor of Philosophy, 2014 Directed By: Professor Eugenia Kalnay, Department of Atmospheric and Oceanic Science, The University of Maryland, College Park Professor Sumant Nigam, Department of Atmospheric and Oceanic Science, The University of Maryland, College Park The effectiveness of the breeding method in determining growth and decay characteristics of certain solutions to the Kortweg-de Vries (KdV) equation and the Nonlinear Schrödinger equation is investigated. Bred vectors are a finite amplitude, finite time generalization of Leading Lyapunov Vectors (LLV), and breeding has been used to predict large impending fluctuations in many systems, including chaotic systems. Here, the focus is on predicting fluctuations associated with extreme waves. The bred vector analysis is applied to the KdV equation with two types of initial conditions: soliton collisions, and a Gaussian distribution which decays into a group of solitons. The soliton solutions are stable, and the breeding analysis enables tracking of the growth and decay during the interactions. Furthermore, this study with a known stable system helps validate the use of breeding method for waves. This analysis is also applied to characterize rogue wave type solutions of the NLSE, which have been used to describe extreme ocean waves. In the results obtained, the growth rate maxima and the peaks of the bred vector always precede the rogue wave peaks. This suggests that the growth rate and bred vectors may serve as precursors for predicting energy localization due to rogue waves. Finally, the results reveal that the breeding method can be used to identify numerical instabilities. Effective simulation of diurnal variability is an important aspect of many geophysical data assimilation systems. For the Martian atmosphere, thermal tides are particularly prominent and contribute much to the Martian atmospheric circulation, dynamics and dust transport. To study the Mars diurnal variability (or thermal tides), the GFDL Mars Global Climate Model (MGCM) with the 4D-Local Ensemble Transform Kalman Filter (4D-LETKF) is used to perform a reanalysis of spacecraft temperature retrievals. We find that the use of a “traditional” 6-hr assimilation cycle induces spurious forcing of a resonantly-enhanced semi-diurnal Kelvin waves represented in both surface pressure and mid-level temperature by forming a wave 4 pattern in the diurnal averaged analysis increment that acts as a "topographic" stationary forcing. Different assimilation window lengths in the 4D-LETKF are introduced to remove the artificially induced resonance. It is found that short assimilation window lengths not only remove the spurious resonance, but also push the migrating semi-diurnal temperature variation at 50 Pa closer to the estimated "true" tides even in the absence of a radiatively active water ice cloud parameterization. In order to compare the performance of different assimilation window lengths, short-term to long-term forecasts based on the hour 00 and 12 assimilation are evaluated and compared. Results show that during NH summer, it is not the assimilation window length, but the radiatively active water ice cloud that influences the model prediction. A “diurnal bias correction” that includes bias correction fields dependent on the local time is shown to effectively reduce the forecast root mean square differences (RMSD) between forecasts and observations, compensate for the absence of water ice cloud parameterization, and enhance Martian atmosphere prediction. The implications of these results for data assimilation in the Earth’s atmosphere are also discussed. BREEDING ANALYSIS OF GROWTH AND DECAY IN NONLINEAR WAVES AND DATA ASSIMILATION AND PREDICTABILITY IN THE MARTIAN ATMOSPHERE By Yongjing Zhao Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2014 Advisory Committee: Professor Eugenia Kalnay, Chair Professor Sumant Nigam, Co-Chair Professor Kayo Ide Professor Massimo Ricotti Professor Steven J. Greybush Professor Balakumar Balachandran, Dean’s Representative © Copyright by Yongjing Zhao 2014 ii Dedication To my husband Wang, Chenxi Parents Zhao, Wenxing Hou, Jianxin iii Acknowledgements First, I would like to thank my advisors, Professor Eugenia Kalnay and Prof. Sumant Nigam, who have offered me great opportunities to study on a lot of interesting projects. Prof. Kalnay guided me throughout the last three years and Prof. Nigam guided me throughout the first two years of my Ph.D. study with their patience and insights. I am deeply thankful to have them as my advisors and friends. I also would like to thank the other members of my committee: Professor Bala Balachandran, Professor Massimo Ricotti, Professor Kayo Ide, and Professor Steven J. Greybush. They have been very helpful throughout the time when I studied in Department of Atmospheric and Oceanic Science and very supportive during the graduation process. I would like to thank all my group members and colleagues John Wilson, Dr. Ross Hoffman, Dr. Alfredo Ruiz-Barradas, Dr. Bin Guan, Dr. Argyro Kavvada, Ayan Moitra, Arseny Zakharov, Dr. Ji-Sun Kang, Ying Zhang, Yan Zhou, Dr. Guo-Yuan Lien, Dr. Daisuke Hotta, Catherine Thomas, Jing Wang, and Yanni Ding, for a lot of helpful discussions. I would like to especially thank Prof. Steven Greybush for helping me getting familiar with the MGCM-LETKF system and working with me very closely in the Mars project (Chapter 3 in the dissertation). I would also like to thank Dr. Christopher Chabalko for contributing in carefully revising the breeding paper (Chapter 2 in the dissertation) and helping me getting familiar with the MATLAB codes. I would like to thank Dr. Meng Cheng for contributing to the numerical solution to the Nonlinear Schrödinger Equation in the breeding paper (Chapter 2 in the dissertation). I would like to thank everyone in Department of Atmospheric and Oceanic Science, the UMD iv WEATHER-CHAOS group, CDI group, EMARS group and other facilities on campus for providing knowledge or service. My deepest thanks are to my husband, Chenxi Wang, and my parents, Jianxin Hou, and Wenxing Zhao, for providing endless love and support during my Ph.D. study. v Table of Contents Dedication ..................................................................................................................... ii Acknowledgements ...................................................................................................... iii Table of Contents .......................................................................................................... v List of Tables .............................................................................................................. vii List of Figures ............................................................................................................ viii Chapter 1: Introduction ................................................................................................. 1 Chapter 2: Breeding Analysis of Growth and Decay in Nonlinear Waves................... 5 2.1 Introduction ......................................................................................................... 5 2.2 Introduction to the Breeding Method .................................................................. 9 2.3 Breeding Analysis of Solutions of the Korteweg-de Vries (KdV) Equation .... 12 2.3.1 Solitons ...................................................................................................... 12 2.3.2 Brief Introduction to the Korteweg-de Vries (KdV) Equation .................. 13 2.3.3 Breeding on the Korteweg-de Vries (KdV) Equation................................ 14 2.4 Breeding Analysis of Solutions of the Nonlinear Schrödinger Equation (NLSE) ................................................................................................................................. 21 2.4.1 Brief introduction to rogue waves.............................................................. 21 2.4.2 Brief Introduction to the Nonlinear Schrödinger Equation (NLSE) and the Numerical Solution ............................................................................................. 22 2.4.3 Breeding Analysis ...................................................................................... 23 2.4.3.1 Breeding Application to Modulational Instability .................................. 24 2.4.3.2 Breeding Analysis and Detection of Numerical Instability .................... 29 2.5 Conclusions ....................................................................................................... 33 Chapter 3: Data Assimilation and Predictability in the Martian Atmosphere .... 35 3.1 Introduction ....................................................................................................... 35 3.2 Background Information ................................................................................... 39 3.2.1 Thermal Emission Spectrometer (TES) Profiles........................................ 39 3.2.2 Mars Climate Sounder (MCS) Profiles ...................................................... 40 3.2.3 GFDL Mars Global Climate Model (MGCM)........................................... 41 3.2.4 Local Ensemble Transform Kalman Filter (LETKF) ................................ 44 3.2.5 Brief Review of Tides ................................................................................ 48 3.3 Experiment Design............................................................................................ 52 3.4 Results ............................................................................................................... 55 3.4.1 6-hour assimilation: spurious resonance .................................................... 55 2.4.2 Short window assimilation: resonance removal ........................................ 63 3.4.3 Impact of Water Ice Cloud Parameterization in Analyses and Forecasts .. 66 3.4.4 Impact of window length on the Martian atmosphere forecast skill .......... 71 3.4.5 Bias Correction .......................................................................................... 79 3.4.6 TES vs. MCS ............................................................................................. 87 3.5 Summary and Discussion .................................................................................. 89 Chapter 4: Summary and Lessons Learned ................................................................ 97 Appendix A: The Indian Ocean Dipole: A Monopole in SST .................................. 100 A.1 Introduction .................................................................................................... 100 vi A.2 Introduction to SST Principal Components ................................................... 105 A.3 Datasets .......................................................................................................... 117 A.4 Variability in ENSO-filtered Indian Ocean SSTs: No Evidence for a Zonal- dipole..................................................................................................................... 119 A.5 The 1994 and 1997 Fall SST Anomalies: IOD Events? ................................ 123 A.6 Seasonally Evolving Subsurface Variability in the Indian Ocean: ENSO’s Influence ............................................................................................................... 128 A.7 Recurrent Modes of Variability in Filtered Indian Ocean SSTs ................... 132 A.8 Conclusion .................................................................................................... 142 Bibliography ............................................................................................................. 146 vii List of Tables Table 3.1 Surface pressure semi-diurnal tide amplitudes at Viking Lander sites, all are shown in percentages ........................................................................................ 56 Table A.1 Correlations between Pacific and Atlantic PCs ................................... 113 Table A.2 Sensitivity Analysis: T0 is the primary analysis, which was perturbed as listed below to assess robustness of the extracted variability modes. In all cases, SST was detrended using the nonstationary SST Secular Trend mode (see footnote 1), and is referred as ‘DeTrend’ here . ‘0’ in the ‘Rotated’ column indicates that loading vectors were not rotated ........................................................................... 133 Table A.3 Sensitivity Results: Column 1 shows the leading modes identified in the primary analysis (T0); numbers following the name indicate the percentage variance explained by the mode and its rank, respectively. Columns 2–6 list attributes of the leading modes in the five sensitivity tests (T1–T5), with the three- slash numbers indicating correlation between the test case and primary analysis PC, the percentage variance explained by that mode, and its rank in the test analysis. Asterisks indicate ‘no match’. All correlation coefficients are significant at the 0.001 level ............................................................................................................. 135 Table A.4 Analog Counts: Number of “analogs” in six extended-EOF analyses of Indian Ocean SST variability during 1902-2007. An analog is identified when the absolute value of the PC of any one mode is larger than that of all others by at least one unit; note PCs are orthonormal. Mathematically, if | |PCi| - |PCj| | > 1.00 for all j not equal than i, an analog is counted. For consistent evaluation, only analogs of the two leading modes are counted ....................................................................... 135 viii List of Figures Figure 2.1 Illustration of the Breeding method ....................................................... 11 Figure 2.2 The control run (left) and bred vectors (right) of the first collision for the two-soliton case. Time increases from top to bottom, showing the solitons before, during, and after collision ....................................................................................... 15 Figure 2.3 The spatial global growth rate for the two-soliton solution case for KdV equation. Initial perturbation is set as [0.01,0.01,......]T ........................................... 17 Figure 2.4 Similar as Figure 2.3, but with different Gaussian initial perturbation. (a): N(0.01,10-4) (b): N(-0.01,10-4).. ....................................................................... 17 Figure 2.5 The control run (left) and bred vectors (right) for the Gaussian case. The solution time increases from top to bottom. Solitons emerge from the Gaussian distribution and then propagate along one direction ............................................... 19 Figure 2.6 The control run (left) and the bred vectors (right) for the Gaussian initial condition. Isolated intervals of soliton collision causing large transient instabilities, identified by trough and crest pairs, are shown ...................................................... 20 Figure 2.7 The spatial global evolution growth rate of the Gaussian initial condition case for KdV equation ........................................................................................... 20 Figure 2.8 First case: (a) analytical and (b) ETDRK4 numerical rogue wave solution ( 0 / 2iA  ) of the nonlinear Schrödinger equation ............................. 24 Figure 2.9 First case: (a) ETDRK4-based space-time bred vector solution and (b) spatial global growth rate of unstable mode of the nonlinear Schrödinger equation ................................................................................................................................. 28 Figure 2.10 First case: spatial slice at (a) x = 0; crest, and (b) x = π/2; trough of rogue wave solution (Figure 2.7b): control run (blue dashed line), the bred vectors (magenta dotted line), and the global growth rate (red solid line) .......................... 29 Figure 2.11 Second case: (a) analytical and (b) ETDRK4 numerical rogue wave solution ( 0 2iA  ) of the nonlinear Schrödinger equation ................................ 30 Figure 2.12 Second case: (a) bred vectors and (b) spatial global growth rate of rogue wave solution ( 0 2iA  ) to the nonlinear Schrödinger equation ............. 31 ix Figure 2.13 Second case: spatial slice at (a) x = 0; crest, and (b) x = -π/2; trough of rogue wave solution (Figure 2.10b): control run (blue dashed line), the bred vectors (magenta dotted line), and the global growth rate (red solid line) .......................... 31 Figure 2.14 Space-time solution based on ETDRK4 method (a), (c) and Strang- splitting method (b), (d) of the first and second type of unstable mode of the nonlinear Schrödinger equation .............................................................................. 32 Figure 2.15 (a) Strang-splitting-based bred vector solution of the first rogue wave solution and (b) the spatial global growth rate ........................................................ 32 Figure 3.1 Surface pressure semi-diurnal tide amplitude spatial distribution (a)~(c) and wave 2 amplitude (d)~(f). (a) and (d) are from free runs, (b) and (e) are from 6- hr TES assimilation and (c) and (f) are from 6-hr MCS assimilation: all are with ice cloud parameterization. Viking Lander sites location and their corresponding tides amplitudes from model simulations (in percentage, after normalization) are denoted in (a)~(c) ................................................................................................................. 56 Figure 3.2 Equatorial daily average temperature analysis increment for different assimilations. (a) represents the original 6-hr TES assimilation assimilating observation at hour 00, 06,12 and 18 with localization length scale = 600; (b) is the same as (a), but with doubled localization length scale; (c) is the same as (a), but assimilates observation at hour 03, 09, 15 and 21; (d) is the 12-hr TES assimilation assimilating observation at hour 00 and 12 with localization length scale = 600. All are with ice cloud parameterization ........................................................................ 58 Figure 3.3 Surface pressure semi-diurnal tide wave 2 amplitude for TES assimilations. (a) is the 6-hr TES assimilation updating observation at hour 00, 06, 12 and 18 with a doubled localization length scale=1200, (b) is the 6-hr TES assimilation updating observation at hour 03, 09, 15 and 21 with a localization length scale=600. All are with ice cloud parameterization ..................................... 59 Figure 3.4 Surface pressure equatorial and 23°N semi-diurnal tide amplitude (a)~(c) and wave 2 phase (d)~(f). (a) and (d) are from the original 6-hr TES assimilation assimilating observation at hour 00, 06,12 and 18 with localization length scale = 600, (b) and (e) are the same as (a) and (d), but with doubled localization length scale; (c) and (f) are the same as (a) and (d), but assimilates observation at hour 03, 09, 15 and 21: all are with ice cloud parameterization ..... 60 Figure 3.5 Surface pressure diurnal tide amplitude spatial distribution (a)~(c) and wave 1 amplitude (d)~(f). (a) and (d) are from free runs, (b) and (e) are from 6-hr TES assimilation and (c) and (f) are from 12-hr TES assimilation: all are with ice cloud parameterization ............................................................................................ 61 Figure 3.6 The Schematic illustration of the 6-hr 3D-LETKF, 6-hr 4D-LETKF, 2- hr 4D-LETKF and 1-hr 3D-LETKF Panel (a) and (b) are coming from Greybush et x al. (2011). In 4D-LETKF, observations (combined together as red/pink line segments to discern between different time slots) are compared to forecasts (solid blue lines) at their correct hours rather than just to the 6-hr/1-hr forecast as in the 3D-LETKF .............................................................................................................. 64 Figure 3.7 Surface pressure semi-diurnal tide wave 2 amplitude for TES and MCS assimilations. (a) ~ (c) are from TES assimilations, (d) ~ (f) are from MCS assimilation. (a) and (d) are 6-hr LETKF assimilation, (b) and (e) are 6-hr LETKF assimilation, and (c) and (f) are 1-hr LETKF assimilation. All are with ice cloud parameterization ...................................................................................................... 65 Figure 3.8 7-day composite daily-average surface pressure bias (a)~(c) and level 10 (50Pa) temperature bias (d)~(e) over sol 397-403 for 6-hr (a), (d) 2-hr (b) (e) and 1- hr (c) (f) of TES assimilation. All are with ice cloud parameterization ................ 66 Figure 3.9 Zonal average temperature semi-diurnal tide wave amplitude for free run and MCS 1-hr assimilation. (a) and (c) are from free run, w/o and with ice cloud parameterization, respectively; (b) and (d) are from TES assimilation, w/o and with ice cloud parameterization, respectively. Regions where ice cloud should be present (100~10Pa) are denoted between dashed lines ...................................... 70 Figure 3.10 50Pa temperature semi-diurnal tide wave 2 amplitude for free run and TES assimilations. (a) and (e) are from free run; (b)~(d) and (f)~(h) are from TES assimilations, 6-hr, 2-hr and 1-hr respectively. (a)~(d) are without ice cloud parameterization, and (e)~(f) are with ice cloud parameterization ........................ 71 Figure 3.11 Performance of the MGCM-LETKF and the free run (without assimilating observation) during NH summer (~Ls 100), evaluated by comparing RMS differences of forecasts from TES observations. (a) is without ice cloud parameterization and (b) is with ice cloud parameterization. Time slots without observation are shaded in grey. Time-averaged RMS difference is denoted in the legend ...................................................................................................................... 74 Figure 3.12 Performance of the MGCM-LETKF during NH summer (~Ls 100), evaluated by comparing RMSE differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o and with ice cloud parameterization ...................................................................................................... 76 Figure 3.13 Performance of the MGCM-LETKF and the free run (without assimilating observation) during NH fall (~Ls 200), evaluated by comparing RMS differences of forecasts from TES observations. (a) is without ice cloud parameterization and (b) is with ice cloud parameterization. Time slots without observation are shaded in grey. Time-averaged RMS difference is denoted in the legend ...................................................................................................................... 78 xi Figure 3.14 Performance of the MGCM-LETKF during NH fall (~Ls 200), evaluated by comparing RMSE differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o and with ice cloud parameterization ...................................................................................................... 79 Figure 3.15 7-day composite bias pattern (shaded colors) and observation track (dots) over sol 397-403 at level 10 for 1-hr TES assimilation (a)~(d), and 6-hr TES assimilation (e)~(h), at different times a day: hour 00 (a), (e); hour 06 (b), (f); hour 12 (c), (g); and hour 18 (d), (h). Different colors of dots represent observation path at different day. All are without ice cloud parameterization ................................. 82 Figure 3.16 7-day composite daily-average bias over sol 397-403 at level 10 for 1- hr (a) and 6-hr (b) TES assimilation. Both are without ice cloud parameterization ................................................................................................................................. 82 Figure 3.17 Performance of the 1-hr MGCM-LETKF TES reanalysis during NH summer (~Ls 100), w/o and with bias correction (both sliding window daily average bias correction and sliding window diurnal bias correction) in a RMSE sense. Time slots without observation are shaded in grey. All are without ice cloud parameterization. Time-averaged RMS is denoted in the legend ........................... 83 Figure 3.18 Performance of the MGCM-LETKF during NH summer (~Ls 100), evaluated by comparing RMSE differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o bias correction, with sliding window bias correction and with sliding window diurnal bias correction. All are without ice cloud parameterization. (a) is TES assimilation, (b) is TES assimilation only considering 30S-30N, (c) is MCS assimilation .............. 84 Figure 3.19 Bias pattern of 120 hour forecast in 1-hr TES LETKF (a)-(d) and 6-hr TES LETKF (e)-(h) runs. Left panels are original forecast without bias correction; second column panels correct the sliding window daily-average bias every 1hr/6hr; third column panels correct the sliding window diurnal bias every 1hr/6hr; the above-mentioned panels are all without RAC parameterization; right panels are original forecast without bias correction, but with RAC parameterization ............ 86 Figure 3.20 The same as Figure 3.19, but for MCS assimilations ......................... 87 Figure 3.21 MCS (a) and TES (b) reanalysis of equatorial temperature difference (3pm-3am). Both are with ice cloud parameterization .......................................... 88 Figure 3.22 Free run (a), TES (b) and MCS (c) reanalysis of zonal mean temperature diurnal tide amplitude. All are with ice cloud parameterization ........ 89 Figure A.1 Principal components (PCs) of the seven leading modes of recurrent spatiotemporal SST variability in the Pacific (shaded) during 1900-2009. Tick xii marks on the vertical axis denote three standard deviations. Smoothed PCs after 50 applications of a 1-2-1 smoother are represented in solid black lines .................. 108 Figure A.2 Spatiotemporal evolution of recurrent SST variability in the Pacific: Seasonal lead-lag regressions of the leading principal components (PC 1-7) ...... 110 Figure A.3 Principal components (PCs) of the four leading modes of recurrent spatiotemporal SST variability in the Atlantic (shaded) during 1901-2008. Tick marks on the vertical axis denote three standard deviations. Smoothed PCs after 50 applications of a 1-2-1 smoother are represented in solid black lines .................. 114 Figure A.4 Spatiotemporal evolution of recurrent SST variability in the Atlantic: Seasonal lead-lag regressions of the leading principal components (PC 1-4) ...... 116 Figure A.5 Evaluation of the basis for Indian Ocean Dipole (IOD) mode of SST variability from both recent (1958-1998) and century-long (1900-2007) SST records; the former period is identical to that analyzed by Saji et al. (1999): (a) Boreal fall SST regressions of the western-IOD box (W-IOD, marked) based SST index; (b) likewise for the eastern-IOD box (E-IOD, marked) index; note, the dipole mode index (DMI, Saji et al.) is obtained by subtracting the E-IOD index from the W-IOD one. Fall regressions of the W-IOD and E-IOD indices on ENSO- filtered SSTs are shown in (c)-(d); ENSO-filtered SSTs are generated, as noted in the text, by removing the spatiotemporally-evolving ENSO influence. Corresponding regressions for the century-long period are shown in the right column. In all cases, regressions are computed on detrended SST (following footnote 1). Contour interval and shading threshold is 0.05K, and positive values are shown in warm colors; the zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black. ............................................. 119 Figure A.6 Two additional constructs of ENSO variability are filtered from Indian Ocean SSTs prior to 716 computation of the regressions of the W-IOD and E-IOD based SST indices: ENSO variability is constructed 717 from the canonical and non-canonical ENSO components only (i.e., without contribution from Pacific biennial 718 variability) in the left panels, and from the widely used index of ENSO’s mature-phase variability (Nino3.4 SST 719 index) in the right panels. Boreal fall regressions are for the century-long period (1900-2007); rest as in Figure A.5. Regressions significant at the 0.05 level are contoured in black ....... 120 Figure A.7a-j Origin of the Fall 1994 and 1997 Indian Ocean SST anomalies: Observed anomalies are in the top panels (a,f); ones reconstructed from ENSO’s full influence (canonical, non-canonical, biennial modes) in the next two panels (b,g); those reconstructed from ENSO’s partial influence (canonical, non-canonical modes) in the middle panels (c,h); those reconstructed using ENSO’s mature-phase index (Nino3.4 SST index) in the next from the bottom panels (d,i), while the fully ENSO-filtered anomalies are in the bottom panels (e,j). Detrended SSTs are the starting point in all cases. The western and eastern IOD box is displayed in all xiii panels to facilitate IOD recognition. Contour interval and shading threshold is 0.15K, and positive values are shown in warm colors; the zero-contour is omitted in all panels ........................................................................................................... 126 Figure A.7k-n Origin of the Fall 1994 and 1997 Indian Ocean upper-ocean (0- 381m) heat content anomalies: Those from SODA 2.1.6 reanalysis (Carton et al. 2000) are in the top panels (k,m), while the fully ENSO-filtered ocean heat content anomalies are in the bottom panels (l,n). The anomalies are displayed after 4 applications of the 9-point smoother (smth9 in GrADS) with a shading threshold/interval of 2x108 J/m2; see color bar. Detrended ocean heat-content is the starting point in all cases. The western and eastern IOD box is displayed in all panels to facilitate IOD recognition ...................................................................... 128 Figure A.8 Fall-centered seasonal lead-lag regressions of 3 versions of the E-IOD SST index on upper-ocean (0-381m) heat content (OHC, from SODA 2.1.6 ocean reanalysis) and 850 hPa winds (from ERA-40 atmospheric reanalysis): Regressions of the ‘raw’ (i.e., unfiltered) index are in the top row; of an index constructed from just ENSO related IO variability in the middle; while those of the ENSO-filtered E- IOD index are in the bottom row. Note, the index and regressions are obtained from the same data set in all cases. The summer (JJA) anomalies (first column) are the 1-season lead regressions while the winter (DJF) ones (last column) are the 1- season lag regressions; contemporaneous (fall, SON) ones are in the middle column. The OHC regressions are shaded/contoured in warm (positive) and cool (negative) colors with an interval/ threshold of 108 J/m2: see color bar. Wind vector scales are shown below each panel. All data sets were detrended (following footnote 1) and the zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black .................................................................... 131 Figure A.9 Principal components (PCs) of the two leading modes of recurrent spatiotemporal SST variability in the Indian Ocean are displayed in the upper panel. The PCs are obtained from rotated, extended-EOF analysis of Pacific- uninfluenced, non-seasonal Indian Ocean SST variability during 1902-2007; PC-1 and PC-2 explains 10.9% and 5.5% of the variance, respectively. Tick marks on the vertical axis are drawn every three standard deviations. The autocorrelation structure of the twos PCs and the IOD and Subtropical-IOD (SIOD) indices is shown in the lower panel ...................................................................................... 138 Figure A.10 Spatiotemporal evolution of recurrent SST variability in the Indian Ocean: Seasonal lead-lag regressions of the leading principal components (PC-1 and PC-2) are shown in the left two columns, with time increasing downward. Fall- centered regressions (i.e., of the fall PC values) are shown to facilitate comparison with the IOD structure which is robust in boreal fall (SON). Indian Ocean SSTs are filtered for the Pacific’s influence (as in T0 analysis) prior to computation of regressions. Fall-centered regressions of the IOD, and winter-centered (DJF) regressions of the Subtropical IOD on unfiltered Indian Ocean SSTs are shown in the far right columns. The western and eastern IOD box is displayed in all panels to xiv facilitate IOD recognition. Regressions are for the 1902-2007 period, and detrended SSTs are the starting point in all cases. Contour interval and shading threshold is 0.05K, and positive values are shown in warm colors; zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black ...................................................................................................................... 139 Figure A.11 Seasonal regressions/correlations of the Indian Ocean SST principal components on the Pacific-uninfluenced Indian Ocean SST and adjoining continental precipitation: JJA (top), SON (middle), and DJF (bottom). SST regressions are shaded and contoured in warm (positive) and cool (negative) colors with an interval/threshold of 0.05K. Correlations with CRU TS3.1’s half-degree continental precipitation are shown with a shading threshold of 0.10 and an interval of 0.05; see color bar. All data sets were detrended (following footnote 1) and the zero-contour is omitted in all panels. Regressions/correlations significant at the 0.05 level 798 are contoured in black ................................................................... 141 1 Chapter 1: Introduction Ships sinking into the ocean with little warning, primarily attributed to rogue or freak waves have brought attention to extreme ocean waves. Extreme waves are actually effects of energy focusing, with large concentration of energy gathered from surrounding smaller perturbations. Modulational instability, or Benjamin Feir (BF) instability, suggests that a rogue wave is a kind of unstable wave which “sucks” energy from surrounding waves and grows into a large wave through nonlinear processes. Therefore, as part of the work of my Ph.D. dissertation, I worked on a project focused on studying the instability and predictability of rogue waves (expressed by the nonlinear Schrödinger equation (NLSE)), as well as the shallow water waves (expressed by Kortweg-de Vries equation, or KdV equation) using breeding (Toth and Kalnay, 1993, 1997). Bred Vectors (BVs) are particularly suited to detect instabilities in a nonlinear system as well as model the behavior of the growing errors. Breeding results show that the soliton solutions of KdV equation have periods of growth and decay only during soliton collisions, but are overall stable. For the rogue wave prediction in the NLSE, the peaks of bred vector always precede the rogue wave peaks, indicating that bred vectors may serve as precursors to predict rogue wave energy localization. The breeding method was found to be successful in indicating stability for the stable solitons of the KdV equation, impending growth for the unstable solutions of the NLSE equation, and detecting numerical instability. This part of work is discussed in Chapter 2. 2 The effectiveness of bred vectors in characterizing chaotic behavior has also been established in weather-forecasting: bred vectors are representative of forecast errors due to inexact initial conditions. The application is not only for the Earth, but also for Mars (Greybush et al., 2013). After identifying regions and seasons of instability of the Martian atmosphere (Greybush et al., 2013), as well as studying Martian atmospheric variability using data assimilation (Greybush et al., 2012), the forecast skill of Martian atmosphere has become equally important. This is the second project I am currently working on. Following the basic configuration in Greybush et al. (2012)’s reanalysis, initially we use “traditional” 6 hour assimilation cycle like that used on the Earth atmosphere, but we found it is problematic -- spurious resonance is found in the Kelvin waves represented in both surface pressure and mid-level temperature. This is because the 6 hour assimilation cycle introduces a wave 4 pattern in the diurnal averaged analysis increment that acts as a "topographic" stationary forcing. As a result, we developed short assimilation window scheme, within which the spurious resonance has been removed. It is found that short assimilation window lengths are also found to push the migrating semi-diurnal temperature variation at 50 Pa closer to the estimated "true" tides even in the absence of a radiatively active water ice cloud parameterization. Short-term to long-term forecast skills at different seasons using different window length are evaluated and compared, suggesting that during NH summer, it is not the assimilation window length, but the radiatively active water ice cloud that influences the model prediction. A “diurnal bias correction” that includes bias correction fields dependent on the local time is shown to effectively enhance Martian atmosphere prediction. Our discovery of assimilation-induced 3 artifacts in atmospheric motions may have relevance to the Earth as well: although the role of thermal tidal waves is much larger in Mars than in the Earth, the widely used 6hr assimilation windows in operational data assimilation does have the potential of forcing resonant atmospheric modes, and should be studied in detail in future studies. This part of work is discussed in Chapter 3. Prior to these two projects, I was involved in a project under the guidance of Prof. Nigam. The primary goal of the project is to estimate the effects of the principal component (PC)s of Pacific, Atlantic and Indian Ocean (IO) SSTs on Chinese summer rainfall, as part of a broader effort to understand the reasons for the drying trend of northern China in recent decades. Following Guan and Nigam (2008, 2009)’s studies about the PCs of Pacific and Atlantic SST, we extracted the IO SST PCs, and found that the IO SST and Pacific SST are strongly correlated. Especially, a zonal- dipole structure in interannual variations of the tropical IO SSTs – the Indian Ocean Dipole (IOD), is strongly correlated with the ENSO signal. After removal of ENSO’s influence from IO SSTs, the IOD poles are essentially uncorrelated in the ENSO- filtered SSTs in both recent (1958-1998) and century-long (1900-2007) periods. As a result, we proposed a monopole SST index to account for this factor, and regressions of the eastern IOD pole on upper-ocean heat content suggest that internally-generated basin variability can have zonal-dipole structure at the subsurface, other than the surface. A spatiotemporal analysis using the extended-EOF technique, which targets spatial and temporal recurrence and extracts modes (rather than patterns) of variability, also does not support the existence of zonal-dipole variability at the surface, however, this analysis did yield a dipole-like structure in the meridional 4 direction in boreal fall/winter, when it resembles the Subtropical IOD pattern (but not evolution timescale). This part of work is discussed in Appendix A due to limited content. 5 Chapter 2: Breeding Analysis of Growth and Decay in Nonlinear Waves 2.1 Introduction Toth and Kalnay (1993, 1997) introduced breeding as a simple method to determine the characteristics of atmospheric instabilities. Breeding was initially used for the analysis of forecasting errors in atmospheric data assimilation. The usable prediction interval of detailed atmospheric models is often limited by chaos. The sensitivity of the models to small perturbations was characterized by application of the breeding method. The effectiveness of bred vectors in characterizing chaotic behavior has been established in many fields. In weather-forecasting, bred vectors are representative of forecast errors due to inexact initial conditions, as noted not only for the Lorenz models (Annan, 2004; Balci et al., 2012; Bowler, 2006; Hallersberg et al., 2010; Pazó et al., 2010, 2011, 2013; Primo et al., 2008), but also for complex models both for the Earth (Carrassi et al., 2007; Corazza et al., 2003; Hoffman et al., 2009; Keller et al., 2010; Magnusson et al., 2008) and Mars (Greybush et al., 2013; Newman et al., 2004). The application of breeding has also been extended to the ocean. Ocean instabilities of different types and time scales have been identified through the application of bred vector analysis (Hoffman et al., 2009). In this case, the work was based on a global ocean model driven by wind stress and surface fluxes from the NCEP/NCAR reanalysis (Kalnay et al., 1996). Toth and Kalnay (1997) also suggested that the breeding method could be applied to a coupled ocean-atmosphere system to isolate the El Niño -Southern Oscillation (ENSO) coupled instabilities from fast instabilities, which are relevant to weather systems. In this case, the breeding method is employed to advantage of the early nonlinear saturation of convective 6 instabilities compared to baroclinic instabilities. In coupled systems (Cai et al., 2003; Tang and Deng, 2010), breeding was shown to be effective in forecasting ENSO events in a perfect modeling scenario through atmospheric data assimilation and ensemble forecasting. This provided a basis for applying the breeding method to a fully coupled atmosphere-ocean general circulation model (GCM) with the ultimate goal of improving operational climate predictions through ensemble forecasting and data assimilation (Yang et al., 2006; Yang et al., 2009). Breeding has also been used as an undergraduate education tool to study the Lorenz system (Evans et al., 2004). The occurrence and duration of regime changes of the model were predicted by measuring bred vector growth. Motivated by prior work, this effort will focus on the efficacy of the breeding method in identifying the stability characteristics of solutions to the Kortweg-de Vries (KdV) equation and the nonlinear Schrödinger equation (NLSE). These applications provide new insight into the growth and decay of nonlinear waves and demonstrate additional application areas for the breeding method. Furthermore, this work helps gain a new perspective in analyzing the stability of KdV solutions, as well as gaining insight into energy localization associated with rogue waves. By construction, bred vectors are a nonlinear, finite time generalization of the construction used for leading Lyapunov vectors (LLV). Their logarithmic growth rate is equivalent to the LLV exponent. In the context of atmospheric modeling, two states that are initially similar will diverge as they evolve in a chaotic system. During this process, stable perturbations will decay while unstable perturbations will grow until they dominate the state. The differences between the perturbed and unperturbed states 7 are bred vectors. They can capture the most unstable modes manifesting physical instabilities. In breeding, the differences between the two solutions are periodically rescaled in order to maintain them as small. Fast-growing perturbations are thus captured in a "breeding cycle", i.e. the fastest-growing modes are naturally selected (Toth and Kalnay, 1993, 1997). If the locations, shapes, and time periods of the bred vectors are known, the dynamics or physical interactions that govern these phenomena can be further explored. As a result, breeding can serve as a useful tool to explore and characterize chaotic systems (Evans et al., 2004). Shallow water waves can be modeled by the KdV equation. This equation admits solitary-wave solutions, also referred to as solitons. The KdV equation was derived to describe unidirectional shallow water wave propagation. It has been widely used as a prototype of an exactly integrable nonlinear system in different fields, including fluid dynamics, theoretical physics, pure mathematics, and related ones. The KdV equation serves as a model to describe weakly nonlinear surface waves in shallow water, providing insights to the study of surface gravity waves (Boussinesq, 1877; Korteweg et al., 1895). This application is seen as one of the most complete and convincing conformations of the Inverse Scattering Transform (IST) theory of KdV (Hammack and segur, 1974, 1978). The KdV equation has also been used to model internal solitons in the ocean through experimental observations (Osborne and Burch, 1980; Segur and Hammack, 1982). Recent advances in sensor data and computational analysis have provided data related to extreme waves in oceanic environments. The definition of a "rogue wave" has been well established as a wave whose height is more than twice the significant 8 wave height (SWH). On the basis of this definition and current theories, rogue waves occur more often than predicted by traditional statistical analysis. The new findings suggest that loss of ocean vessels, which have been anecdotally attributed to extreme waves, could possibly be mitigated through the use of predictive technologies (Carrassi et al., 2007; Ansmann et al., 2013; Rosenthal et al., 2005; Baschek et al., 2011; Kharif et al., 2003). Due to the intermittent nature of oceanic rogue waves, and the difficulty in gathering data about them, several governing mechanisms have been proposed (Forristall, 2005). These include diffractive focusing and the interaction of current effects and wind waves (http://www.opc.ncep.noaa.gov/perfectstorm/mpc_ps_rogue.shtml). Recently, it has been proposed that rogue waves can be generated by a combination of several physical factors such as strong winds and strong currents as well, rather than by a single distinct cause. The Benjamin-Feir instability mechanism is commonly associated with rogue wave formation in deep water (Calini and Schober, 2002; Henderson et al., 1999; Kharif et al., 2001; Onorato et al., 2001; Pelinovsky et al., 2000). This mechanism has been tested in fiber optical systems, and similar results have been obtained in controlled wave tank experiments. Recently, a laboratory water tank experiment successfully demonstrated a rogue wave consistent with the Benjamin-Feir instability observed in a NLSE model (Cho, 2011). The Benjamin-Feir instability is used to describe energy transfer from smaller adjacent perturbations to a larger unstable wave. This unstable wave reaches a maximum, and then other nonlinear effects, such 9 as breaking in ocean waves, lead to energy dissipation. This modulational instability in deep water conditions can be modeled by the NLSE and is often the basis for a theoretical study of rogue waves (Henderson, 1999; Dudley et al., 2008; Kibler et al., 2010; Ablowitz et al., 2001; Carter and Contreras, 2008; Schober, 2011; Islas and Schober, 2011; McLaughlin and Schober, 1992; Osborne et al., 2000; White and Fornberg, 1998). While the breeding method can be applied without prior knowledge of the instability mechanism, the application to a NLSE wave solution can help gain insights into the growth rate of theoretically predicted rogue waves. The investigation of solution stability for the KdV equation is still ongoing. It has long been recognized that, in most cases, KdV solutions are soliton solutions. In this work, the effectiveness of the breeding method is evaluated in the identification of the stability of soliton solutions to the KdV equation, and rogue wave solutions to the NLSE. Breeding is also applied to identify instabilities in a numerically unstable system. The breeding method is described in detail in Section 2.2, and the results of the application of breeding to the KdV, NLSE, and numerically unstable system are presented in Sections 2.3 and 2.4. Section 2.5 contains a summary and a discussion of the results. 2.2 Introduction to the Breeding Method Lyapunov vectors (LVs) are widely used to characterize the directions of highest expansion and contraction in dynamical systems. A Lyapunov exponent can be obtained by tracking a perturbation in a system. The perturbation is periodically renormalized to avoid strongly nonlinear interactions and computational overflow. If a system is unstable, all perturbations will converge towards the fastest growing 10 perturbation, namely the leading Lyapunov vector (Nayfeh and Balachandran, 1995). Bred vectors are obtained by a nonlinear, finite time generalization of the method that is used to obtain the leading Lyapunov vector, and therefore, “breed" the fastest nonlinear perturbations (Yang et al., 2006; Kalnay, 2003; Yang et al., 2008). Like the first leading Lyapunov vector (LLV), bred vectors are independent of the definition of norm, and could represent the direction of maximum sustainable growth (or minimum decay) within a system without the presence of external forcing (Kalnay et al., 2003). The breeding perturbations or bred vectors (BVs) are obtained through the following steps: (a) add a small arbitrary perturbation to the initial condition at the initial time t0, (b) integrate the nonlinear model for a period of time (t1-t0) with the unperturbed (control) initial condition and the perturbed initial condition, (c) subtract the unperturbed (control) forecast from the perturbed forecast, and (d) scale down the difference between the perturbed and unperturbed runs to keep them the same size (using any norm |dx|) as the initial perturbation, (e) add the rescaled perturbation to the analysis at t1, and finally, (f) repeat steps (b)-(e) forward in time. For small perturbation sizes, the BV becomes the same as the local leading Lyapunov vector. In the presence of fast and slow growing instabilities, the amplitude of the initial perturbation and the rescaling period can be used to separate the leading bred vectors corresponding to the instabilities (Peña and Kalnay, 2004). An initial perturbation is illustrated as dx0 in Figure 2.1. The bred vector is the difference between the perturbed and unperturbed (control) forecasts integrated in time using the nonlinear model M , tctfndttcndttft xxxMxMdx   )()( . (2.1) 11 The breeding window length, or rescaling interval is Tw = ndt , where n is the number of time steps before each rescaling and dt is the model time step. The bred vector amplification factor is then defined as the size of the bred vector after n steps divided by its original size ||/|| 0dxdx t , and the local growth rate per unit time g at time t can be written as, g = 1ndt  ln dxt dx0     . (2.2) It is recalled that a chaotic system is indicated by a Lyapunov exponent > 0. In a similar fashion, a non-positive long term logarithmic bred vector growth rate indicates a stable system, and a positive long term logarithm bred vector growth rate indicates the system is unstable. Figure 2.1 Illustration of the breeding method. 12 2.3 Breeding Analysis of Solutions of the Korteweg-de Vries (KdV) Equation 2.3.1 Solitons Dispersive partial differential equations (PDEs) have long been used to model the behavior of idealized vibrating media in the absence of frictional or other dissipative forces, for example, waves on a string, on the surface of water, or in air. The Korteweg-de Vries (KdV) equation and the non-linear Schrödinger equation both belong to this category. Solutions to these PDEs have two seemingly contradicting properties. The first property is to be conservative with time, which means after a considerable long time t, the state of the solution remains "similar" to the initial state; which also indicates that energy is conserved. The second property is dispersion, which refers to the fact that waves move at different speeds based on their frequency. Although energy is conserved, different components of the solution may convect at different speeds. Therefore, the amplitudes of the solution after a long period of time may appear to be attenuated as compared to the initial state. Solitons, in which the dispersive tendency is perfectly balanced with nonlinear effects, are also possible. As a self-reinforcing wave (wave packet or pulse), a soliton's shape remains invariant as it travels. Solitons are found as the solutions to a class of nonlinear dispersive partial differential equations, including the KdV equation. An important property of the soliton solutions studied is that the shape remains invariant before and after interaction with another soliton. For a non-dissipative system, the amplitude, shape, and velocity of the soliton are conserved after soliton 13 collision. The breeding method will be used to demonstrate the persistence of the soliton properties in the following subsection. 2.3.2 Brief Introduction to the Korteweg-de Vries (KdV) Equation In 1895, Korteweg and de Vries derived a nonlinear PDE that bears their names, which serves as a mathematical model of shallow water waves, explaining the soliton phenomena mathematically (Korteweg et al., 1895): 06  xxxxt  . (2.3) They argued that this equation could explain Russell’s experiments. This equation stated that the rate of change of the wave height in time is determined by two factors, a nonlinear one (6x) and a dispersive one (xxx). As a result, the solution of this equation can be seen as a balance between nonlinear and dispersive effects. It was not until 1965 when Zabusky, and Kruskal published the numerical solution of the KdV equation (Zabusky and Kruskal, 1965) that Russell’s observations and the KdV equation became widely studied by mathematicians, physicists and engineers interested in water waves. Gary Deem, Zabusky, and Kruskal also introduced a modified KdV equation ( 0 xxxxt  ) and showed the solitary wave interaction process (see (Zabusky, 1984)). The following Section is aimed at applying breeding method to this modified KdV equation, and discussing the results. 14 2.3.3 Breeding on the Korteweg-de Vries (KdV) Equation In this work, numerical solutions to the modified KdV equation ( 0 xxxxt  ) have been determined through the Fourier spectral/ exponential time-differencing fourth-order Runge–Kutta (ETDRK4) scheme due to Kassam and Trefethen (Kassam and Trefethen, 2005). Spatial periodic boundary conditions have been modeled over [-,]. Since the modified KdV equation is a function of time and space, the breeding method takes the following definition. Let )0,(xbv be the initial perturbation, and   bv(x,t) x ft  xct be the bred vector after time t. As described above, the breeding window length is Tw = ndt . The bred vector amplification factor is calculated as, amplificationfactor = bv(x,Tw )bv(x, 0) , (2.4) and the growth rate is computed as 1 ln( )GrowthRate amplificationfactorn dt  . (2.5) The breeding method is then numerically applied to analyze the stability of the solutions. Stable solutions are indicated by a non-increasing bred vector. Unstable solutions are indicated by a bred vector of increasing amplitude in time. In this case, a small difference in the initial condition will result in large and unpredictable discrepancies in the forecasts. In this work, the perturbation is rescaled every eight time steps. The initial perturbation is set as [0.01,0.01,......]Tfor simplicity. Other initial perturbations, such as uniform stochastic distributions, have also been tested and found to yield similar results, as discussed later. Since the bred vectors resemble the leading Lyapunov vectors localized in both space and time, they do not depend on 15 the norm and the scaling period (Toth and Kalnay, 1993, 1997). In the current study, the independence of bred vectors on the rescaling interval and norm has been verified in numerical experiments with different rescaling periods and different scaling factors multiplied by the initial perturbations. Figure 2.2 The control run (left) and bred vectors (right) of the first collision for the two-soliton case. Time increases from top to bottom, showing the solitons before, during, and after collision. The first case examined is the two-soliton solution. In this case, the integration time interval is set as 10-6, and the total integration time is 0.025. Initially, two solitons with different amplitude exist. As the solution evolves, dispersion causes the larger amplitude soliton to convect faster than the smaller amplitude as expected (Russell, 1844). The first soliton collision process occurs at time (t = 0.001-0.004) as 16 seen in the left portion of Figure 2.2. The bred vectors have non-zero values during the soliton collision as observed in the right portion of Figure 2.2. In order to characterize possible instabilities throughout the spatial domain, the global growth rate is considered. The global growth rate is defined as ( , )1 ln ( ,0) bv x tGlobalGrowthRate n dt bv x        , (2.6) this can be defined using a quadratic norm as ( , ) ( , )1 ln ( ,0) ( ,0) T T bv x t bv x tGlobalGrowthRate n dt bv x bv x         . (2.7) The norm of the bred vector shows that instabilities only occur during the collision process as seen in Figure 2.3. During the collision, the perturbation growth rate decreases rapidly followed by an even faster increase associated with a transient instability. After the collision, the perturbation growth rate quickly attenuates. After an initial spin-up time, two characteristic “trough” and “crest” signatures emerge as seen in Figure 2.3. An adjacent trough and crest pair is characteristic of growth and decay. The positive growth rate, that is, larger amplitude of the bred vectors, is indicative of a transient instability. The global growth rate demonstrates an initial transient instability while the perturbed initial condition transitions to a soliton before and during the first collision. The second (and following) collisions show an initial decrease in perturbation energy followed by a compensating growth. The long term net BV energy growth remains zero indicating that the soliton solutions are stable except during the initial transition from the perturbed state and during collisions. 17 Figure 2.3 The spatial global growth rate for the two-soliton solution case for KdV equation. Initia perturbation is set as [0.01,0.01,......]T. Figure 2.4 Similar as Figure 2.3, but with different Gaussian initial perturbation. (a): N(0.01,10-4) (b): N(-0.01,10-4). 18 Note that the initial perturbation is set as [0.01,0.01,......]Tin Figure 2.3. When the solitons collide, the perturbation growth rate decreases quickly first, then grows even more quickly, and decays again, keeping a smooth shape. Both perturbed bred vectors in Figure 2.4 show a very similar shape in the bred vector growth rate compared to Figure 2.3 but with initial high-frequency small-amplitude perturbations, indicating that the bred vectors are very robust to changes in the initial perturbations. The second type of initial condition used was a Gaussian pulse. In this case, the integration time interval was 10-6, and the total integration time was 0.05. A Gaussian initial condition decays into solitons that convect in the positive direction as shown in the left portion of Figure 2.5. These solitons mutually interact and also cause transient periods of stability and instability as shown by the amplitude of the bred vectors in the right portion of Figure 2.6. After the Gaussian initial condition breaks into five solitons with different amplitudes (and speeds), the larger solitons collide with smaller ones. A detailed view of the system is shown at a later time in the left portion of Figure 2.6. In this case, time intervals of isolated soliton interactions are demonstrated. The corresponding bred vectors, indicating transient instability and overall long term stability, are shown in the right portion of Figure 2.6. The global growth rate exhibits many trough and crest pairs. Each one is associated with a soliton collision, causing transient instability, as shown in Figure 2.7. The transient instabilities are absent in the absence of soliton collisions. 19 Figure 2.5: The control run (left) and bred vectors (right) for the Gaussian case. The solution time increases from top to bottom. Solitons emerge from the Gaussian distribution and then propagate along one direction. The global growth rate of the Gaussian initial condition contains many trough and crest pairs with varying amplitude. Large and medium signatures occur when two groups of solitons collide simultaneously while smaller signatures correspond to the collision of a single pair of solitons. The growth rate curve decreases to zero after collision, demonstrating long term stability. 20 Figure 2.6 The control run (left) and the bred vectors (right) for the Gaussian initial condition. Isolated intervals of soliton collision causing large transient instabilities, identified by trough and crest pairs, are shown. Figure 2.7 The spatial global evolution growth rate of the Gaussian initial condition case for KdV equation. 21 2.4 Breeding Analysis of Solutions of the Nonlinear Schrödinger Equation (NLSE) 2.4.1 Brief introduction to rogue waves Rogue waves are often referred to relatively large and spontaneous ocean surface waves that represent a huge threat to vessels in the ocean. Waves whose height is more than twice the significant wave height (Hs or SWH) are usually defined as rogue waves. Traditionally, the SWH is defined as the mean wave height, the vertical distance between trough and crest, of the highest third of the waves. Another definition is four times the standard deviation of the surface elevation. Linear ocean wave theory suggested that the rogue wave occurrence was extremely low since the rogue wave probability is so small. Only until decades ago, people started realizing that rogue waves could occur much more often than they thought when solid evidence became available. Rogue waves could also be explained by linear random theory, in which the wave can be seen as the superposition of a series of sinusoidal wave components consisting of different propagation directions and frequencies. The surface therefore follows Gaussian distribution suggested by the central limit theorem. Linear random wave theory also suggests that for narrow-banded unidirectional wave, wave height and thus wave occurrence follow a Rayleigh distribution (Longuet-Higgins, 1952). However, it has been suggested by observations that the occurrence of large waves whose heights are greater than twice the significant wave height may not be properly described by the classical linear theories (Skourup et al., 1997; Stansell, 2005). Considering nonlinearity, for general wave with broad directional spreading, it is hard to describe rogue wave occurrence statistically, and as a result, relevant studies are 22 mainly based on numerical simulations. The NLSE, which provides model equations based on the narrow-band wave envelope approximation, is introduced in Section 2.4.2. For broad directional waves, the probability of rogue wave occurrence is found to approximately follow Gaussian statistics (Gramstad and Trulsen, 2007). Nonlinear probability of rogue waves is found to be dependent of water depth and directional spreading angle (Xiao, 2013). 2.4.2 Brief Introduction to the Nonlinear Schrödinger Equation (NLSE) and the Numerical Solution The nonlinear Schrödinger equation is used to model many physical nonlinear systems including waves in the deep ocean, fiber optics, nonlinear acoustics, hydrodynamics, quantum condensates and many other nonlinear media and systems (Cho, 2011; Dudley et al., 2008; Kibler et al., 2010; Ablowitz et al., 2001; Malomed et al., 2005; Nore et al., 1993; Som et al., 1979). In this Section, two different initial conditions corresponding to rogue wave solutions to the NLSE are analyzed. The breeding method is used to characterize the stability of each one. Furthermore, the predictive nature of the breeding method will be addressed. The nonlinear Schrödinger equation for deep-water waves can be written as: 2( ) 0t g x xi c v       , (2.8) where  is a complex field. For deep water waves, the absolute value of  can be interpreted as the amplitude of the wave envelope. The constant, real coefficients are given by 2/ ,8 ,2 2 002 0 0 0 0 kvkkcg    , (2.9) 23 and the associated linear, deep water dispersion relation is given by 020 gk . (2.10) where 0 is the frequency of the carrier wave, 0k is the wave number, and gc is the group velocity of the wave. The solution to the nonlinear Schrödinger equation is a complex narrow-banded envelope. The amplitude of the wave train (x,t) is given by:   (x,t)(x,t)eik0i0t const, (2.11) where k0 is the wavenumber of the carrier wave, 0 is the frequency of the carrier wave and cg is the linear group speed. The solution to the nonlinear Schrödinger equation is represented in terms of its modulus ( , ) ( , )A x t x t and its phase, ( , )x t as:   (x,t) A(x,t)ei(x,t ). (2.12) Transforming Eq. (8) in a frame moving with the group velocity of the wavetrain gives the nonlinear Schrödinger equation for deep water used in this work (Vitanov et al., 2013); that is, 2 2 2 2 0i t x         . (2.13) Numerical solutions are constructed by using the ETDRK4 method (Kassam and Trefethen, 2005) similar to the numerical solutions of KdV equation. A Fourier differentiation matrix is used to solve the spatial derivatives. 2.4.3 Breeding Analysis The breeding method is applied to the solutions presented above for the nonlinear Schrödinger equation. Since the solutions are complex, the control run solutions are 24 2cx , and the corresponding bred vectors are 2f cx x . For both cases, the spatial domain has been discretized into 512 points, and the bred perturbation has been rescaled every eight time steps similar to the analysis presented for the KdV equation. The initial perturbation is also chosen as [0.01, 0.01, ...]T. The ETDRK4 scheme is found to exhibit better numerical stability compared to the Strang-splitting scheme (Press et al., 1992), and this is discussed in Section 2.4.3.2. 2.4.3.1 Breeding Application to Modulational Instability Figure 2.8 First case: (a) analytical and (b) ETDRK4 numerical rogue wave solution ( 0 / 2iA  ) of the nonlinear Schrödinger equation. The first solution examined, which is based on a hyper elliptic functions and an eigenvalue 0 / 2iA  , is expressed as 2 0 2 2 20 0 0 0 2 0 0 cos( 2 )sec (2 ) 2 tanh(2 )( , ) 2 cos( 2 )sec (2 ) iA tA x h A t i A tx t A e A x h A t       (2.14) 25 This rogue wave is periodic in space, but not in time. The requirement of periodic boundary conditions in the ETDRK4 method is satisfied by choosing a spatial domain of [0, 2] and setting 0 1/ 2A  . The numerical time integration begins at t = 10, and the initial condition takes the form of the corresponding analytical solution. The numerical solution matches the analytical solution for the first localization, as shown in Figure 2.8a. Due to temporal periodicity, a wrap around effect is exhibited in the numerical solution as shown in Figure 2.8b (Press et al., 1992). The breeding method was applied to the solution. Since the solution is periodic in time due to the wrap around effect, only the first trough and crest pair is analyzed. The global growth rate in the EDTRK4 scheme (Kassam and Trefethen, 2005), shown in Figure 2.9b exhibits a modulational instability. The first crest and trough pair occurs during t ~ [-3, 3]. After t = 7, the growth rate takes on an exponential trend. Since the growth rate of the norm of the unperturbed solution (control run) remains zero, the growth rate for the norm of the bred vector indicates instabilities and possible energy transfer. Details regarding the growth rate of the control run are as follows: growth rate for the norm of the NLSE solution: (2.15) where * is the complex conjugate of . 26 At t = 0 (time step = 500), the spatial solution can be characterized in terms of crests and troughs. The crests are identified at x ~ (0, )&(7) and the troughs are identified at x ~ (7) as seen in Figure 2.8b. As the solutions are symmetric about x =  only the first half, x ~ (0, ) is the focus of the current study. For t = 0, the solution reaches its maximum at x = 0, 2 and minimum at x = /2 and 3 /2. At /4 and 7 /4, the solution is close to the carrier wave amplitude A0. In order to focus on the solution crest, a spatial slice through x = 0 is shown in Figure 2.10a. The temporal solution, global growth rate, and bred vector in the crest region area are plotted together. A similar analysis in the trough region x =  /2 is shown in Figure 2.10b. Correspondingly, only the first trough and crest pair is analyzed since the solution is periodic in time. In the crest region, there is a mild trough and crest pair in the global growth rate indicating strong growth. The trough and crest pair is not as significant for the spatial slice through the trough region. The second solution analyzed is also a rogue wave, with similar features to the first. Based on hyperelliptic functions with an eigenvalue of 0 2iA  , A0, the solution is expressed as 2 0 2 2 0 0 2 0 2 0 0 2 cos(4 2 ) 2 sin(4 2 ) ( , ) 1 cos(4 2 ) 2 cosh(2 ) iA t A t i A t x t A e A t A x           ( .16) where A0 = 1/2. This solution is periodic in time, but not in space. However, it is symmetric about x = 0. Therefore the solution can be treated as spatially periodic over a domain that is symmetric about x = 0. The application of periodic boundary condition necessarily requires that the solution is 27 identically valued at the domain boundary. For a spatial domain of [-, ], the periodic boundary conditions required by Fourier differentiation/ETDRK4 are satisfied. The initial condition corresponds to the analytical solution at t = 0. The solution remains numerically stable for long periods of integration (Figure 2.11). Again, the analytical and numerical solutions are found to match each other for long periods of integration. Since the solution is periodic in time, an interval of t = [0, 10] is identified for analysis. Large modulations are apparent in the second initial condition. This case is highly modulated and contains very few regions of pure carrier wave amplitude A0. Troughs can be identified by the light blue color (-~/8, 3/8~) and crests are identified above the green color regions (/8~/8) as seen in Figure 2.11b. This analysis is only focused on the first half (-~) as the solution is symmetric about x = 0. For the first trough and crest pair the solution reaches its maximum at x = 0. Minimum values occur close to x = -/2 and /2. At x = -38 and 3/8, the amplitude remains at an intermediate value between maximum and minimum. The breeding method is applied to the solution, as shown in Figure 2.12a, and the global growth rate is computed as shown in Figure 2.12b. 28 Figure 2.9 First case: (a) ETDRK4-based space-time bred vector solution and (b) spatial global growth rate of unstable mode of the nonlinear Schrödinger equation. Several features in the global growth rate can be better understood when compared to the solution behavior. A spatial slice at x = 0 isolates the crest region. The temporal solution, bred vector amplitude, and global growth rate are shown for this spatial location in Figure 2.13a. The global growth rate maximum appears several steps before the first bred vector peak; coincident with the solution crest observed in Figure 2.11a, just after t = 0. In this case, the growth rate serves as a powerful indicator of the impending rogue wave instability. Similar solution behavior is observed during troughs. A spatial slice at x = -/2 isolates the trough region. The temporal solution, bred vector amplitude, and global growth rate are shown for this spatial location in Figure 2.13b. The bred vector minimum appears simultaneously with the growth rate maximum, and the growth attains a maximum slightly before the solution trough. However, for practical considerations, the first crest is most significant for the prediction of oceanic rogue waves. The bred vector growth rates are always ahead of, and therefore serve as a precursor to, the crests. 29 Figure 2.10 First case: spatial slice at (a) x = 0; crest, and (b) x = π/2; trough of rogue wave solution (Figure 2.7b): control run (blue dashed line), the bred vectors (magenta dotted line), and the global growth rate (red solid line) 2.4.3.2 Breeding Analysis and Detection of Numerical Instability In this section, the use of breeding analysis to identify instabilities in numerically unstable systems is evaluated. The first type of rogue wave solution ( 0 / 2iA  ) to the NLSE is solved using both the ETDRK4 scheme and the Strang-splitting scheme (Strang, 1968) as shown in the left and right portions of Figure 2.14. In this case, the ETDRK4 solution is numerically stable for long periods of integration. The Strang- splitting technique shows similar results until t ~ 35, after which a numerical instability develops and corrupts the solution. Breeding was applied to the solution computed with the Strang-splitting technique as shown in Figure 2.15a. The bred vector shows typical behavior until around t ~ 35, after which the behavior changes dramatically due to the numerical instability. A qualitative change can also be identified in the global growth rate, as shown in Figure 2.15b. In the period from t ~ [0, 35], the global growth rate exhibits low frequency oscillations associated with the 30 previously analyzed rogue wave solutions. After t ~ [0, 35], the oscillations in the global growth rate increase in frequency. The change is less dramatic than that exhibited by the bred vectors. This can be attributed to the spatial averaging which is incorporated in the global growth rate. The observed qualitative changes in the analysis suggest that the numerical instability in the Strang-splitting solutions can be detected by the bred vectors (Figure 2.14c). The Strang-splitting scheme for the second rogue wave solution is also numerically unstable after t ~ 20 (Figure 2.14d). Bred vector analysis also detects the numerical instability in the second rogue wave solution. However, in this case, the numerical instability is not as severe. Figure 2.11 Second case: (a) analytical and (b) ETDRK4 numerical rogue wave solution ( 0 2iA  ) of the nonlinear Schrödinger equation. 31 Figure 2.12 Second case: (a) bred vectors and (b) spatial global growth rate of rogue wave solution ( 0 2iA  ) to the nonlinear Schrödinger equation. Figure 2.13 Second case: spatial slice at (a) x = 0; crest, and (b) x = -π/2; trough of rogue wave solution (Figure 2.10b): control run (blue dashed line), the bred vectors (magenta dotted line), and the global growth rate (red solid line). 32 Figure 2.14 Space-time solution based on ETDRK4 method (a), (c) and Strang-splitting method (b), (d) of the first and second type of unstable mode of the nonlinear Schrödinger equation. Figure 2.15 (a) Strang-splitting-based bred vector solution of the first rogue wave solution and (b) the spatial global growth rate. 33 2.5 Conclusions In this study, bred vectors have been used to study the stability of soliton solutions to the KdV equation, rogue wave solutions to the NLSE equation, and numerically unstable solutions. This novel analysis brings forth new perspectives for demonstrating the stability of soliton solutions, as well as identifying a precursor of energy localization for rogue waves. In one of the cases studied, a Gaussian initial condition breaks into a train of solitons. Solitons originate from the initial center peak and propagate along one direction. The solitons mutually interact due to dispersion in the system. The bred vector analysis shows transient energy growth and decay during soliton collisions. The same analysis indicates a stable system in the absence of soliton collisions. This serves as a validation case for the ability of the breeding method to detect unstable growth. Breeding was also applied towards the identification of rogue wave solutions of the nonlinear Schrödinger equation. These rogue wave solutions, induced by modulational instabilities, are characterized by trough and crest pairs in the bred vector growth rates. Since the growth rate for the global norm of the NLSE solution remains zero (see Eq. 2.15), the growth rate for the norm of bred vector may indicate the possible energy transfers and transient instabilities. In the cases studied, the growth rate maxima always preceded the solution crest. This suggests the global growth rate and bred vector could serve as indicators of a rogue wave peak. Global growth rate troughs are not as easily associated with minima of the rogue wave envelope. Consider a rogue wave happening on Feb 1933 lifting the sea into huge 15m swells, which has a period of 14.8s, a calculated wavelength of 342m and a 34 calculated speed of 23m/s (Bascom, 1964; Mayo, 1997). We could approximately treat this case as the second type of solution to NLSE discussed in Section 2.4.3.1 since they are both identified periodic in time. As described in Figure 2.12a, the BV global growth rate maximum when the first “energy localization” occurs is approximately 18% (compared with the temporal period) ahead of the control run rogue wave peak, which is approximately 62m, or 2.7s. This may not be exciting enough since 2.7s ahead is too short for the vessels to take action, although 2.7s, or 62m, is 18% of a wave cycle. In a full numerical model including the carrier waves evolution, it should be possible to predict and understand the transfer of energy to and from the carrier waves and thus predict the appearance of rogue waves further ahead. The bred vectors and global growth rate are found to exhibit qualitative changes in the presence of a numerical instability. The bred vectors show high frequency oscillations similar to that in the unstable solution. The growth rate fluctuations transition to lower amplitudes and higher frequencies. Due to spatial averaging in the global growth rate, the high frequency spatial content is reduced to more modest fluctuations than the frequencies present in bred vectors. In summary, the breeding method was found to successful in indicating stability for the stable solitons of the KdV equation, impending growth for the unstable solutions of the NLSE equation, and detecting numerical instability. In the future, breeding will be applied to study the energy transfer in a full field rogue wave solution, rather than the envelope solutions shown here. This will be accomplished through the study of the Dysthe equation (Islas and Schober, 2011), which contains both KdV and NLSE terms in addition to damping terms. 35 Chapter 3: Data Assimilation and Predictability in the Martian Atmosphere 3.1 Introduction Atmospheric thermal tides are planetary-scale gravity waves, with periods that are harmonics of the solar day (Chapman and Lindzen, 1970; Forbes, 1995). Tides are particularly prominent in the Martian atmosphere so that temperature, wind fields and surface pressure, exhibit strong dependence on local solar time (LT). This is because in Mars, the very low atmospheric density near the surface, (two orders of magnitude smaller than on Earth) and the low surface thermal inertia (compared to that of the terrestrial oceans) lead to very large diurnal variations in surface temperature and a resulting strong atmospheric response. The rapid radiative damping time results in a strong coupling between surface and boundary layer temperature. Aerosols, particularly dust, also lead to significant solar heating within the atmosphere, particularly during dust storm episodes. (Zurek, 1976; Wilson and Hamilton, 1996; Zurek et al., 1992). Thermal tides are particularly significant for the understanding the Martian atmosphere due to their broad impacts on Martian atmosphere dynamics (Forbes et al., 2002), and they also have impacts on dust and water transport, which can in turn affect the atmospheric stability, dynamics, chemistry and thermal structure (Zurek et al., 1992). They are also of great importance for interpreting spacecraft observations. Tides can be decomposed into westward propagating migrating (sun-synchronous) waves that are strongly driven by solar heating, and non-migrating waves resulting from zonal variations in thermo-tidal forcing. Zonal modulation of forcing can arise 36 from longitudinal variations of topography, albedo and surface thermal inertia, as well as from spatial variations in radiatively active aerosols including dust and water ice clouds (Zurek, 1976; Wilson and Hamilton, 1996). The study of Martian thermal tides and stationary waves has been mostly based on temperature observations provided by the Mars Global Surveyor (MGS) Thermal Emission Spectrometer (TES), providing full latitude and longitude coverage profiles (Smith et al., 2001). TES observes in a sun-synchronous orbit that provides twice- daily temperature observations (2 am and 2 pm local solar time, LT), which is insufficient to fully characterize the basic features of the thermal tides. TES provides temperature observations, mostly in nadir viewing mode, resulting in atmospheric soundings from 0-35 km (surface to 10 Pa), with a vertical resolution of ~10 km (Smith et al., 2001). The twice-daily TES observations at 2am and 2pm LT provide some constraint on the tides, but far from enough. As a result, a Data Assimilation System (DAS), which optimally combines temporally and spatially irregular observations with a background state from a global climate model (GCM) forecast, is needed for thermal tidal studies. Mars data assimilation systems have already been developed to yield a dynamically consistent set of dynamical and thermal fields providing the basis for detailed investigations of the Martian dynamical features and circulations system (Houben, 1999; Zhang et al., 2001; Montabone et al., 2005, 2006; Lewis et al., 2007; Hoffman et al., 2010; Lee et al., 2001; Greybush et al., 2012). The resulting TES reanalyses give an estimation of the evolving state of the Martian atmosphere during the MGS mission. Thus, the thermal and dynamical atmospheric variables obtained 37 through data assimilation provide a basis for analyzing Martian atmospheric variability, circulation features, dust transport, etc. Lewis et al. (2007) performed the first complete set of TES-based reanalysis, named “UK Reanalysis”, made many contribution to the understanding of Martian atmospheric properties such as atmospheric tides (Lewis and Barker, 2005), dust storms variability (Montabone et al., 2005), tropical water ice clouds (Wilson et al., 2008), and predictability (Rogberg et al., 2010). Since then, ensemble-based methods have been developed and matured (see Kalnay, 2010), and have been identified as a key advancement in meteorological data assimilation. Ensemble-based methods estimate the flow-dependent background error covariances, i.e., the background uncertainty, from the ensemble spread, and this has led to significant improvements in Numerical Weather Prediction (NWP) skill (e.g., Kleist, 2012). An example of an ensemble-based data assimilation systems used in Martian atmosphere reanalysis is the Data Assimilation Research Testbed (DART) (Anderson et al., 2009) used by Lee et al. (2011), where TES radiances were assimilated into the MarsWRF GCM (Richardson et al., 2007). Some atmospheric features such as the circulation, the CO2 cycle, and surface properties were shown to be improved through this ensemble-based assimilation. Montabone et al. (2011) assimilated MGS TES temperature retrievals and column-integrated dust opacities into the UK-LMD-MGCM in use at the University of Oxford and at The Open University in the UK. Like the UK Reanalysis, this assimilation scheme uses the Analysis Correction scheme (Lorenc et al., 1991), and the successive corrections method that demonstrated its robustness in Mars reanalysis, even during the imperfect MGS aerobraking period (Lewis et al., 2007). 38 The resulting Mars Analysis Correction Data Assimilation (MACDA) provides a multi-year reanalysis spanning almost three complete Martian seasonal cycles. An advanced version assimilating both MGS TES and Mars Reconnaissance Orbiter (MRO) Mars Climate Sounder (MCS) (McCleese et al., 2007) profiles into the UK- LMD-MGCM provided the first results of a complete assimilation of both datasets (Montabone et al. 2012), covering almost six complete Mars years, and allowing for the study of the inter-annual variability of the Martian weather and climate with respect to its components including the thermal tides. Hunt et al. (2007) developed the Local Ensemble Transform Kalman Filter (LETKF) system, an advanced data assimilation system that has been successfully tested in the terrestrial data assimilation (Miyoshi and Yamane, 2007; Szunyogh et al., 2008; Miyoshi et al. 2010). The LETKF optimally combines the ensemble forecasts in order to optimally fit the observations. It has several advanced features such as estimating observation errors, which can aid in the improvement of covariance inflation estimation in assimilation studies (Li et al., 2009). Greybush et al. (2012) used a 4D-LETKF scheme that assimilates observations at their time of observation, rather than in 6 hours windows centered at the assimilation time (as had been used in Earth studies), coupled with the Geophysical Fluid Dynamics Laboratory (GFDL) Mars Global Climate Model (MGCM; Wilson, 2011), and generated a new reanalysis of Martian weather and climate including observed and unobserved model variables (temperature, surface pressure, and winds). The 6hr forecast root mean square differences (RMSD) with respect to observations were found to be much smaller in the 4D-LETKF assimilation than in a free-running 39 MGCM, suggesting a successful assimilation. However, we found that the assimilation resulted in eastward propagating semidiurnal tides with unrealistically large amplitudes. One of the main results of this paper is that these unrealistically large tides are due to a spurious resonant forcing associated with the use of 6-hour assimilation windows. In this study, the 4D-LETKF coupled with the MGCM is used to assimilate TES and MCS observations, providing time series of complete Mars fields that reveal the thermal tides features and exploring the impact of the use of different assimilation windows. Background information on the TES/MCS observations, MGCM and LETKF, and a basic description of tides are presented in Section 3.2, the configuration of the system with different assimilation windows is discussed in Section 3.3, results are described in Section 3.4, and Section 3.5 is a summary and discussion. 3.2 Background Information 3.2.1 Thermal Emission Spectrometer (TES) Profiles The Mars Global Surveyor (MGS) spacecraft was launched in November 1996, reached Martian orbit in September 1997, started its primary mission in February 1999, and failed in November 2006, during its third extended mission. The MGS spacecraft followed a sun-synchronous orbit during its primary science mission, providing 12 daily observations each at 2 am and 2 pm local solar time (LT), drifting ~30 degrees in longitude with each subsequent orbit. In the current study, we assimilate temperature retrievals into the MGCM, obtained from TES radiance measurements from 1997 to 2004 (Mars Year (MY) 24 to early MY 27) (Christensen 40 et al., 1992; Christensen et al., 2001). Atmospheric temperature soundings extending from the surface to ~40km (0.1 hPa) with a vertical resolution ranging from ~10km near the surface to ~15km at higher altitudes are retrieved from nadir radiances in the 15 micron CO2 absorption band (Conrath, 2000). The TES nadir radiances are also used to retrieve column dust and ice opacity, but these are not used in this study. Temperature errors are estimated to be around 3K, with higher errors in polar regions and near the surface (Eluszkiewicz et al., 2008). Systematic errors can also vary from orbit-to-orbit. During dust storm seasons, retrieval uncertainties of both temperature profile and dust properties are relatively high, partially due to dust being assumed to be composed of non-scattering particles (Smith et al., 2001). 3.2.2 Mars Climate Sounder (MCS) Profiles As the first science investigation at Mars that is capable of performing a "4- dimensional" systematic study of the key features of Martian weather and climate. Mars Climate Sounder (MCS) on the Mars Reconnaissance Orbiter (MRO), began its primary mission in November 2006, and has been providing vertical profiles of temperature, humidity, pressure, dust content, and clouds of the lower 80km of Martian atmosphere (McCleese et al., 2007). Modern infrared remote sensing techniques are utilized in the MCS mission, including limb sounding and filter radiometers, which enable high vertical resolutions. The Mars Reconnaissance Orbiter is now in its extended mission, which allows the study of interannual variability or climate variations. It did successfully "watch" the planet-scale dust storms that occasionally occur in late southern spring. 41 The Mars Climate Sounder typically observes horizontally from the limb of Martian atmosphere at 3am and 3pm LT, achieving high vertical resolution (i.e., 5km, approximately one half the atmospheric scale-height on Mars, compared to one atmospheric scale-height, or 10km in the TES profiles), providing daily, three- dimensional global weather maps by measuring infrared emitted radiance from the limb of the atmosphere. Besides that, it also uses its azimuth actuator to view the limb 90° “cross track” to the left and right with respect to the orbit and thus provides additional local time observations (Kleinboehl et al., 2013). It can “see” in nine channels spanning the visible and infrared bands. One channel is in the visible and near infrared range (0.3-3.0 microns), which aids in the study of solar energy- atmosphere/surface interaction; the other eight channels are in the thermal infrared range (12-50 microns), measuring temperature, water vapor, pressure, and dust. 3.2.3 GFDL Mars Global Climate Model (MGCM) Originally based on the GFDL SKYHI terrestrial GCM and with an early version of the model described in Wilson and Hamilton (1996), the Mars Global Climate Model (MGCM) is now adapted to the GFDL Flexible Modeling System (FMS), using a finite volume dynamical core (Lin, 2004). It has already been used in Mars atmospheric tides and planetary wave analysis (Wilson and Hamilton, 1996; Hinson and Wilson, 2002; Wilson et al., 2002; Kleinbohl et al. 2013). The MGCM includes surface and subsurface physics essential for the surface temperature calculation; a budget containing gaseous and condensed CO2 that is needed in the average surface pressure annual cycle estimation; inventories of water (regolith, surface frost, water 42 vapor and condensate) and dust aerosol that maintain mass conservation, and a Mellor-Yamada 2.5 boundary layer scheme. Wilson (2011) introduced a parameterization of subgrid-scale gravity waves forced by topography. In this study, we use essentially the same baseline model configurations as for the assimilating experiments in Hoffman et al. (2010) and Greybush et al. (2012). The basic configuration for the MGCM integrations described here is 5 × 6 degree resolution (or ~300 km) with a 36 × 60 latitude-longitude horizontal grid. A hybrid pressure-sigma vertical coordinate system is employed in this model, with 28 vertical levels. A higher resolution (2 × 2.4 degrees, or ~120 km) simulation was also used in Greybush et al. (2012) to test the model sensitivity to horizontal resolution: the features found at both resolutions in most of the free atmosphere such as the polar westerly and tropical easterly jets, and the zonal mean temperature structure are very similar. Although the model is mass-conserving, data assimilation can lead to sources and sinks. To maintain CO2 mass conservation in the assimilations, at each grid point the analysis surface pressure is scaled by the ratio of background atmospheric mass to analysis atmospheric mass after each analysis step, which allows local surface pressure to match the temperature increments while conserving atmospheric CO2 (Greybush et al., 2012). Dust aerosol is a major source of atmospheric heating the Mars atmosphere. A “seasonal dust” dust configuration in which the dust mixing ratio vertical profile is specified as in Forget et al. (2001) and Montmessin et al. (2004), is used in several assimilation experiments in this study. Simulations with an interactively evolving three-dimensional dust distribution were also carried out in this study, with the 43 opacity field represented by radiatively active dust tracers (with an effective diameter of 1.0 micron), although the multiple dust tracers should have had different effective particle sizes, which influence both optical properties and sedimentation velocity. Although model winds advect dust, dust sedimentation is independent of wind and a function of particle size, and it has a spatially and temporally varying vertical distribution that evolves dynamically, as described in Greybush et al. (2012). However, such an evolving 3-dimensional dust distribution did not yield qualitative difference in the tidal analysis at this time of year from the analysis using “seasonal dust”, so that in order to simplify our experiment we only use “seasonal dust”. Water ice clouds also make a significant radiative contribution in the Mars atmosphere (Wilson et al., 2008; Madeleine et al., 2012). The water ice cloud microphysics scheme employed in this model is as described in Montmessin et al. (2004). Cloud microphysical processes are critical in this model. The microphysical processes consider transportation of cloud ice mass, and predict a spatial mean cloud particle size by specifying the ice mass and the number (N) and size (ro) of dust-like cloud condensation nuclei. Assumptions associated with the background dust particle size distribution (PSD) impact the cloud particle size, which controls the sedimentation and cloud optical properties and eventually determines cloud vertical distribution. Another parameter involved in this study is a scaling constant (CR) that determines the cloud opacity. The introduction of CR facilitates the exploration of cloud radiative impact due to the straightforward switch between simulations with radiatively passive (CR=0) and radiatively active clouds (CR > 0). As discussed in 44 Greybush et al. 2012 and Kleinbohl et al. (2013), we can vary CR to “tune” the cloud radiative impact. 3.2.4 Local Ensemble Transform Kalman Filter (LETKF) Data assimilation methods create an analysis that optimally combines information from current observation and a prior short-time forecast, and thus represent the best estimate of the state of the atmosphere. They have two stages: the analysis stage, where the state of the atmosphere is estimated, and the forecast stage. During the analysis stage, observational information and model forecast are statistically combined to yield the optimal estimate of the state of the atmosphere, and this state estimate is then used to initialize the next short-term forecast, called the “first guess” or the background, during the “forecast” stage. In the current study, the MGCM is used to advance the analysis in time to provide a background state in the next analysis, which becomes the best estimate of the state of the atmosphere prior to the utilization of the next time observation. The background error covariance (denoted B) and the observation error covariance (denoted R), characterizing the background and the observation uncertainties determine the optimal weighting of model forecast and observations when producing the analysis. The background error covariance (B) not only quantifies the forecast uncertainty but also determines how one type of observation influences other model forecast variables through the covariance between different variables. In this study, only temperature observations are utilized to update other model variables such as surface pressure, surface stress and wind field. Hoffman et al., (2012) showed the effectiveness of assimilating temperature retrievals in the reduction of RMS error compared to the freely running model with Observing 45 System Simulation Experiments (OSSE). As in Hoffman et al. (2012) we adopted the LETKF scheme (Hunt et al., 2007) which is an efficient formulation of the Ensemble Kalman Filter approach that estimates the current model state within a local region. Each ensemble member with different parameters (dust and water ice cloud properties) is advanced in time by the forecast model to produce the ensemble of forecasts, and therefore yield different atmospheric states (temperature, wind, surface pressure) during the model forecast stage. In this study, 16 ensemble members with dust scales gradually increasing from 1.0 to 2.5 are advanced by the MGCM. Greybush et al. (2012) conducted experiments with both 16 and 32 ensemble members, but found the improvement in forecast accuracy not significant. As a result, we use 16 ensemble members in this study to save computational costs. While obtaining comparable performance as other sequential ensemble-based Kalman Filters applied using the same numerical weather prediction model (Whitaker et al., 2008), the LETKF is computationally efficient because the LETKF assimilation is performed independently at each grid in the center of a local region, making the LETKF computations completely parallel. The LETKF algorithm is now briefly described. Starting from a background ensemble of the model states, where the kth member of the total K ensemble members is represented as a vector kbx , and the ensemble mean is denoted as bx . Their relationship can be expressed as: kb b bx x X  , (3.1) where bX is a matrix, the columns of which represents the background ensemble perturbation departure from the ensemble mean. The column k of bX denotes the 46 difference between model state kbx of the kth ensemble member and the ensemble mean bx . Since K represents the ensemble size, the background error covariance matrix B is expressed as 1 1 T b bK B X X . (3.2) Similarly, the local analysis error covariance is written as Pa = 1 K -1XaXa T = XbPaXbT , (3.3) where aX represents the analysis ensemble perturbation matrix departure from the ensemble mean, and aP is defined as the transform matrix that transfers the background ensemble perturbations into the analysis ensemble perturbations. It can be interpreted as the analysis error covariance in ensemble space (dimension K x K), and derived as Pa = (K -1)I + YbTR-1Yb -1 , (3.4) where bY refers to the background ensemble perturbation matrix transformed into the observation space with the observation operator H; R is the observation error covariance, assumed to be diagonal for simplicity (i.e., the observation error covariance between different variables is neglected). I is the identity matrix with a dimension of K x K. Generally, the ensemble size K is much smaller than the model dimension, and as a result, the calculation in ensemble space highly enhances the computational efficiency. Since we have a limited ensemble size, sampling errors due to spurious long- distance correlations are inevitable. Localization is used to reduce this kind of 47 sampling error, as well as to reduce the computational cost. The cutoff distance is chosen to be 3.65 times the horizontal localization distance L (Miyoshi et al., 2007), and an inverse Gaussian function is applied to gradually decrease the observation impact outside of the cutoff distance: 2 2exp 2 ij R df L       . (3.5) Here L is the localization length scale, and ijd the distance between the i th observation and the j th model grid point (Hunt et al., 2007). Therefore, this function increases the "observational error" gradually. Within the LETKF scheme, "R- localization" is achieved by multiplying the observation error covariance matrix R by Rf . Therefore, uncertainty assigned to the observations increases gradually until they have infinite uncertainty and therefore no influence on the analysis. Greybush et al. (2011) demonstrated that the performance of R-localization using the LETKF in examining the wind-height geostrophic balance has similar performance in terms of analysis accuracy as those of using B-localization (localization is applied to the background covariance matrix B) (Hamill et al., 2001; Houtekamer and Mitchell, 2001) in other ensemble Kalman filters. This localization scale is one of the topics that we investigate, and is discussed in Section 3.4.1. The weight for the ensemble mean can be expressed as wa = Pa Yb( )T R-1 y0 - yb( ), (3.6) so that the analysis ensemble mean is thus obtained as a b b ax x X w  , (3.7) the analysis ensemble perturbations are obtained from: 48   1/21a b a b aK    X X P X W. (3.8) Let kaw be the columns of the matrix that are the results by adding the ensemble mean aw to each of the columns of aW , the resulting "weight" vectors kaw give the linear combinations of the background ensemble perturbations, and add them to the background ensemble mean to obtain the model space analysis in model space: k ka b b ax x X w  . (3.9) The linear combination of the background ensemble perturbations kb aX w can be seen as the analysis increment, in which the relative contribution of each ensemble member to the analysis is available. If each ensemble member contains a unique configuration of model parameters (in this study, the optimal dust scale and the radiatively active water ice cloud content), equation (3.9) can be used to assess the most appropriate physical parameter that best fit the observations. Since the model deficiencies and sampling errors are underestimated in the ensemble Kalman filters, the background ensemble spread should be artificially increased to allow the LETKF to give more weight to the observations due to the probable underestimation of the uncertainties of the ensembles. This process, called "inflation", is discussed in Section 3.4.3. 3.2.5 Brief Review of Tides Tides are a dominant aspect of Mars meteorology, and are much more prominent than those on Earth. Unlike on Earth, Mars tides are strongly influenced by aerosols (dust and water ice clouds). The diurnal and semidiurnal tides have atmospheric components such as surface pressure, temperature, zonal and meridional wind, and 49 surface pressure, and can be further decomposed into westward propagating and eastward propagating waves. The longitude-time dependence of thermal tides can be represented as: , (3.10) and can be expressed in a local time reference frame as: , (3.11) where s is the zonal wavenumber,  denotes east longitude,  represents the temporal harmonic ( = 0 refers to a stationary wave, not discussed in this study.  = 1 stands for the diurnal tide,  = 2 represents the semidiurnal tide), tLT means the local solar time, and Ss,is the phase. Universal and local solar time are related by t = tLT - . Tides with s > 0 (s < 0) correspond to westward (eastward) propagation, and s = 0 means zonally symmetric tides. When s, tides are sun-synchronous without longitude dependence, and are also called "migrating" tides. Tides that do not follow the motion of the sun (s) are called ''non-migrating" tides, which typically result from thermotidal forcing zonal variations. In a local solar time reference, an observed wave number m zonal wave may represent a mix of stationary Am,0, and non-migrating tide components (As,, with (s - m)(Wilson, 2000, Forbes and Hagan, 2000). For example, an observed wave 2 variation includes contributions from the wavenumber 3, westward propagating diurnal tide (A3,1) and wavenumber 1, eastward diurnal tide (A-1,1), as well as higher temporal harmonics (A0,2, A4,2, A1,3, A5,3, ...). Classical tide theory can be used to understand the latitude and vertical structure of the atmopsheric response to imposed thermal forcing. The diurnal response can be 50 represented as a series of vertically propagating waves in the tropical region and vertically trapped waves in the extratropics (Chapman and Lindzen, 1970). The dominant tropical diurnal tide response is a vertically propagating component with a ~33km vertical wavelength (Zurek, 1976; Wilson and Richardson, 2000), while other westward nonmigrating tides are consisting of a series of shorter vertical wavelengths. Since the vertical wavelength of diurnal tide in tropical region and the vertical smoothing length of the TES nadir retrievals have the same order of magnitude, the observed tidal amplitudes are significantly reduced compared to the atmospheric response (Wilson and Richardson, 2000). By contrast, the diurnal tide response in the extratropical region is vertically trapped and so varies little over the depth of the forcing region. There is a class of Kelvin waves that are particularly prominent on Mars. As noted in Hamilton and Garcia (1986) and Zurek (1988), the external mode on Mars is such that the free mode period for Kelvin waves is roughly 24 hours, rather than 33 as on Earth. The proximity of this period to the solar day is what makes resonant enhancement much more likely on Mars than on Earth---after all there is not much forcing available with a 33 hour period. The Kelvin mode contributions to A-s,s (s = 1,2,3…) are relatively non-dispersive and are thereore all potentially nearly resonance for forcing at 24/s periods.1 The detection of the near-resonant Kelvin wave in the atmosphere including a comparison with a MGCM simulation has already been discussed in previous study (Hinson et al., 2008). They can be readily excited by the scattering of the migrating tide by planetary scale topography. Attention has been 1 Note that 1 sol = 1 Mars day = 24 hours 37 minutes. A Martian hour = 1/24 of a sol. 51 paid to the most prominent components of the eastward propagating, in particular to diurnal Kelvin waves (DK1, DK2,... corresponding to s = -1,-2, ...) which are solutions of the Laplace Tidal Equation. They show meridionally broad symmetric structures (Chapman and Lindzen, 1970). Strictly speaking, A-s,s can be represented by a number of different Hough mode with different meridional structure. The Kelvin wave is simply the dominant mode with the broadest meridional structure. The vertical structure of DK1 resembles well the equivalent barotropic Lamb wave, and may be enhanced by resonance (Zurek, 1976; 1988; Wilson and Hamilton, 1996). Both DK2 and DK3 are vertically propagating waves with ~80 and ~60 km vertical wavelength, respectively; the amplitudes of both grow exponentially with height (Wilson, 2000). The zonal wavenumber 2, Semidiurnal East wave (SE2), which could be resonantly enhanced due to different factors, i.e., assimilation window length, is a major point discussed in this study. We are also interested in the westward propagating, wavenumber 2 semidiurnal tides, the dominant contribution to A2,2. This tide has a very long vertical wave length and a broad meridional structure, so it effectively integrates tide forcing distributed over a large range of heights and latitudes (Wilson and Hamilton, 1996; Forbes and Miyahara, 2006).A difficulty in assessing the accuracy of the tides in the reanalysis fields is the sparseness of observations. A fundamental issue for data assimilation is whether the twice-daily spacecraft observations are sufficient to constrain the tides in the assimilating model. One hypothesis is that that the migrating tides are not really constrained by the observations though the non-migrating tides may be more readily 52 constrained. In the current study, the time series of temperature and surface pressure are decomposed into a sequence of harmonics of a solar day. Like in the study of Forbes and Miyahara (2006), in which the fully resolved diurnally varying surface pressure observations provided at the Viking Lander sites are used to study the semidiurnal tides in dusty season, the tidal signals at Viking Lander sites, are also compared to the observations to get quantitatively assessments. 3.3 Experiment Design Assimilation experiments of both TES and MCS profiles are conducted for a 20- sol period during the NH Martian summer, approximately Ls (solar longitude) 94.5°- 103.7°, or sols 388-408 past the perihelion, since during this time atmospheric variability is primarily governed by the thermal tides and aerosol structure, as traveling wave activity is at a minimum. Experimental runs for TES and MCS assimilation are from Mars Year (MY) 25 and MY 31, respectively. Allowing for a 9-sol spin-up period, the starting date of analysis time for Martian tides is sol 397 (Ls 98.6). As the Martian atmosphere shows a low degree of interannual variability in the aphelion season (Smith, 2004), delaying the starting year for MCS profiles to allow for spin-up does not result in substantial differences when comparing the basic tidal features between TES and MCS assimilations at the same period of different years. Recall that we use 16 ensemble members with varying dust opacities increasing linearly from 1.0 to 2.5. A dust configuration with a prescribed vertical profile of dust mixing ratio that evolves with the season and latitude according to an analytic 53 formula (Montmessin et al., 2004), which is labeled as "seasonal dust", is applied in the current study. In this study, observations that are geographically close enough are combined into “super-observations” (Alpert and Kumar, 2007), which significantly reduces the number of observation profiles assimilated in the MGCM, and efficiently reduces computational time, as well as filtering out small scale features that can not be represented by the model and reduce random observation error (Greybush et al., 2012). Only observations with good quality control flags are assimilated as indicated by the Planetary Data System (PDS), from which the TES data was acquired. Quality control is also applied within the LETKF to reject the observations whose departure from the background is larger than a QC threshold. We use the same configuration in the current study as in Greybush et al. (2012): the prescribed TES observation error is 3K. For MCS, the observation error less than 3K is hard-coded to 3K considering probably underestimated errors, and the QC threshold is set as five times the observational error. Like other EnKF systems, the LETKF system is overconfident about the quality of the ensemble forecast because it neglects model errors and the impact of nonlinearities and thus underestimates the analysis and forecast ensemble spread. A significant source of model error is the poor knowledge of radiatively active aerosol in the atmosphere, including dust and water ice clouds. This results in an underestimation of the weights of the observations, and requires “inflating” the ensemble spread at each assimilation time. We use adaptive inflation (Miyoshi, 2011) in our assimilation, in which innovation statistics are used in estimating multiplicative 54 inflation values at each grid point via a Kalman filter (Hunt et al., 2007), which serves as a temporal smoother. Since the satellite path only covers a small fraction of the whole domain, we only inflated the ensemble spread at the grid points where there are observations. Zonal and meridional winds are also inflated at the grid points where there are temperature observations to maintain dynamic consistency. Surface pressure is not inflated since there are no direct observations. For TES profiles, as no observation information available above ~0.l hpa (10 Pa) to constrain the ensemble spread, inflation is not needed in these levels, and a linear decrease in spread inflation values from the top level of the TES observations is applied. For MCS profiles, since observations exist at all model levels, this special treatment is no longer needed. In this study, we use basically the same model configuration as in Greybush et al., (2012). Horizontally 600km localization, vertically 0.4 log P localization and temporally 3h localization. We have found that the semidiurnal east wave response can be sensitive to the horizontal localization radius, which will be discussed in Section 3.4.1. Several post-processing steps are applied at each analysis step. We apply a Shapiro filter to the analysis increments to remove the high frequency noise at high latitudes, and then scale the analysis surface pressure at each grid point by the ratio of background atmospheric mass to analysis atmospheric mass to keep CO2 mass conservation. 55 3.4 Results 3.4.1 6-hour assimilation: spurious resonance Figure 3.1 shows the spatial patterns of the semi-diurnal surface pressure amplitude (in percent) for a free run simulation and TES/MCS data assimilations with a 6hr window. In all cases, the surface pressure field is first normalized by the diurnal mean surface pressure to account for topography. An examination of the tide components confirms that the semi-diurnal tide response is dominated by two modes, the migrating tide (following the Sun) and SE2 (Wilson and Hamilton, 1996), with their amplitudes as a function of latitude shown in the bottom row of Figure 3.1. Wave number 4 patterns (Figures 3.1a-c) can be seen as the interference between these two dominant semidiurnal modes (Figures 3.1d-f), the westward propagating and eastward propagating tides, A2,2 and A-2,2 as described in Section 3.2.4. Compared to the free run spatial patterns, the 6-hr TES assimilation shows stronger wavenumber 4 patterns. At Ls ~100, the surface pressure oscillation at Viking Lander 1 (22.7°N 48.2°W) has an amplitude of about 0.3%~0.4% (Wilson and Hamilton, 1996). The 6-hr TES assimilation shows a much too strong signal compared to the Viking Lander observation data (Table 3.1). The enhanced wave number 4 patterns in the TES and MCS assimilation (Figures 3.1b,c) are due to the significantly enhanced amplitude of A-2,2 (SE2) (Figures 3.1e,f). 56 Figure 3.1 Surface pressure semi-diurnal tide amplitude spatial distribution (a)~(c) and wave 2 amplitude (d)~(f). (a) and (d) are from free runs, (b) and (e) are from 6-hr TES assimilation and (c) and (f) are from 6-hr MCS assimilation: all are with ice cloud parameterization. Viking Lander sites location and their corresponding tides amplitudes from model simulations (in percentage, after normalization) are denoted in (a)~(c). Table 3.1: Surface pressure semi-diurnal tide amplitudes at Viking Lander sites, all are shown in percentages. free run TES 6hr LETKF MCS 6hr LETKF VL1 0.05% 2.5% 0.6% VL2 0.03% 0.6% 0.1% It is clear that the 6hr assimilations have a much larger amplitude of the semidiurnal Kelvin wave than the free run and the Viking observations. This suggests that there is a resonant forcing associated with the analysis increments calculated with 6hour assimilation windows. In order to test this hypothesis we recall that Wilson et 57 al. (1996) showed that zonal modulation of the sun-synchronous wave can be excited by stationary topographic forcing, as introduced in Section 3.2.5. An external forcing with wave number m could scatter the migrating tide (s ) into spatial harmonics  - m and  + m: . (3.12) Equation (2.13) shows that m = 4, corresponding to a wavenumber 4 topographic forcing, can indeed excite a wave of the form cos(-2+2t), or eastward propagating semidiurnal wave. Thus m=2 can yield a scattering of the diurnal westward (migrating) tides into A-1,1 (the Diurnal East wave 1, or DE1). We are now interested in determining whether it is possible that the 6-hr assimilation cycle could introduce a geographically fixed forcing with a wave number 4 pattern that would also excite the non-migrating semidiurnal tide, or SE2 . The equatorial temperature analysis increments in the original 6-hr TES assimilation do show a significant wave- 4 pattern prevalent in several vertical levels, such as surface-800Pa, 100-30Pa and 20-10Pa (Figure 3.2a). This introduces a forcing similar to a wave 4 topographic feature that excites the semidiurnal east wave (Equation 3.13). The use of 6hr assimilation windows and a localization length of 600km produces a maximum amplitude in wave 4 (Figure 2.2a), and thus maximum Kelvin wave forcing. This is because the analysis increments due to observations are computed 4 times a day, and the impact of the observations is constrained by the localization length L=600km, which keeps the analysis increments separated. The amplitude of the wave 4 in the analysis increments is significantly decreased by broadening the localization to 1200km (Figures 3.2b and 3.3a). The wave 4 forcing in the analysis increments can cos( )cos( ) cos[( ) ] cos[( ) ]t m m t m t               58 also be modified by changing the updating hours, thus shifting the longitude phase of the pattern. Figure 3.2c shows the equatorial temperature analysis increment pattern after changing the assimilating hours to 03, 09, 15 and 21 (and subsequently, which spacecraft orbits are emphasized). It is clear that this wave-4 pattern has an eastward shifting trend compared to updating observation at hour 00, 06, 12 and 18, and the resulting resonance in Kelvin wave is also slightly diminished (Figure 3.3b), in this case, the SE2 response is slightly diminished (Figure 3.3b). Actually, the resulting SE2 will be a result of forcing by the “real” physics and the spurious forcing. These two contributions may constructively or destructively interfere. Figure 3.2 Equatorial daily average temperature analysis increment for different assimilations. (a) represents the original 6-hr TES assimilation 59 assimilating observation at hour 00, 06,12 and 18 with localization length scale = 600; (b) is the same as (a), but with doubled localization length scale; (c) is the same as (a), but assimilates observation at hour 03, 09, 15 and 21; (d) is the 12-hr TES assimilation assimilating observation at hour 00 and 12 with localization length scale = 600. All are with ice cloud parameterization. Figure 3.3 Surface pressure semi-diurnal tide wave 2 amplitude for TES assimilations. (a) is the 6-hr TES assimilation updating observation at hour 00, 06, 12 and 18 with a doubled localization length scale=1200, (b) is the 6-hr TES assimilation updating observation at hour 03, 09, 15 and 21 with a localization length scale=600. All are with ice cloud parameterization. The large sensitivity at VL1 site (23°N) tide amplitudes on the semidiurnal Kelvin wave suggests the importance of SDK2, which has a big influence on the resulting interference pattern between the migrating semidiurnal tide and the Kelvin wave, which dominates the spatial structure of the local semidurnal tide (Figures 3.4 a-c). Modeling also suggests that the decline in S2 amplitude at Ls ~ 90 is consistent with the increase in SDK2 amplitude, leading to a greater destructive interference (Wilson and Hamilton, 1996). Obviously, this effect is dependent on the longitude phase of the Kelvin wave, and so may be useful for constraining the assimilation. 60 Changing the update times (0,6,12,18 to 3,9,15,21) significantly changed the phase of the simulated SDK2 (Figure 3.4f); this may have a clear influence on the surface pressure field (Figure 3.4c). On the contrary, doubling the localization radius changed the equatorial SDK2 phase only by ~10° in equatorial to mid-latitude regions (Figure 3.4e), and the phase of the migrating semidiurnal tide is almost stable in equatorial to mid-latitude regions, despite of sharp fluctuations in high-latitude (Figure 3.4d-f), suggesting that large sensitivity of equatorial to mid-latitude phase change is mainly related with the change of the assimilating time. Figure 3.4 Surface pressure equatorial and 23°N semi-diurnal tide amplitude (a)~(c) and wave 2 phase (d)~(f). (a) and (d) are from the original 6-hr TES assimilation assimilating observation at hour 00, 06,12 and 18 with localization length scale = 600, (b) and (e) are the same as (a) and (d), but with doubled localization length scale; (c) and (f) are the same as (a) and (d), but assimilates observation at hour 03, 09, 15 and 21: all are with ice cloud parameterization. 61 Figure 3.5 Surface pressure diurnal tide amplitude spatial distribution (a)~(c) and wave 1 amplitude (d)~(f). (a) and (d) are from free runs, (b) and (e) are from 6-hr TES assimilation and (c) and (f) are from 12-hr TES assimilation: all are with ice cloud parameterization. A 12hr assimilation window with a short localization of 600km produces a wave 2 forcing (Figure 3.2d), as could be expected from the previous arguments. Thus a 12-hr assimilation window induces a wave-2 pattern in the temperature analysis increment (Figure 3.2d), serving as a spurious forcing to modify the wave number 1, eastward propagating diurnal tide (DE1) (Figure 3.5f). The combination of DE1 and t the diurnal migrating tide yield an interference pattern having 2 maxima (Figure 3.5c). It is obvious that the near-resonant Kelvin waves are sensitive to the DA implementation, while the migrating tides do not show this sensitivity. The non- uniqueness of the tide result has a highly dependent on aspects of the of the DA system implementation (the assimilation phase and the localization length scale), the 62 former of which plays a more significant role by altering the forced and resonantly enhanced Kelvin wave considerably. It is important to note that the normal mode Kelvin waves on Mars have periods of ~24/s hours (s is the zonal wave number) as demonstrated in the current section, which also happens to be close to the period forcing associated with harmonics of the diurnal tide. A sol on Mars is 24 hours 37 minutes, and can be divided into 24 "hours", as on Earth. This is why resonant excitement by forcing (legitimate or spurious) can have a significant impact. On Earth, the normal mode period is more like 33/s hours, which does not have any strong source of forcing. Thus this resonance issue is much more significant for Mars than for Earth. As noted previously in Section 3.2.5, the resonance condition is an intrinsic feature of Mars. In a resting atmosphere, the resonance will depend on atmospheric temperature, and the atmosphere is closest to resonance when the atmosphere is coldest, during the aphelion season (~Ls 100 in this study) (Wilson and Hamilton, 1996). Topography yields a way for forcing by the westward propagating (migrating) tides to project onto eastward propagating modes. The assimilation is able to create a spurious forcing that can be amplified by the resonance. Consequently, the amplitude and/or phase of the new response can be notably different from the "natural" forcing, which is what we have seen, both for the 6 hour and the 12 hour updates. Changing the phase of the assimilation updates lead to changes in phase of the eastward- propagating components of the tide response, as we would expect. The migrating components are quite insensitive (Figure is not shown). 63 2.4.2 Short window assimilation: resonance removal We have tested the hypothesis that large amplitude semidiurnal and diurnal Kelvin waves in the reanalysis are spuriously forced by the use of 6 and 12 hour assimilation windows with short localization scales (600km), which concentrate the analysis increments into zonal wave number 4 and 2 respectively. In order to assess whether this spurious forcing can be eliminated (or at least reduced) we replaced the 6 hour assimilation window used in the LETKF with either a 2 hour or a 1 hour assimilation window. This reduction should actually improve the accuracy of the LETKF since shorter assimilation windows maintain a more linear growth of the forecast perturbations, as assumed in the development of the EnKF schemes (Kalnay et al., 2007). Compared to the 3D-LETKF (Figure 3.6a), in which observations are gathered together in the center of a 6-hr window, 4D-LETKF, in which observations are assimilated at their correct time of occurrence within a one-hour window, yields more realistic observation increments. This is particularly important for the Mars reanalysis, since the Martian boundary layer is dominated by large diurnal temperature variations and thus strong thermal tides (Greybush et al., 2012). For a forecast starting at hour 0, and using a traditional 6-hr assimilation (Figure 3.6b), a 9- hour forecast is launched first; observations are then combined together within 1-hr interval: t = 3–3.5 hr, t = 3.5–4.5 hr, t = 4.5–5.5 hr, and then t = 8.5–9 hr. These observation files are then compared to forecasts at t = 3, 4, 5, 6, 7, 8, and 9 hr to create an analysis at t = 6 hr. Figure 3.6c and Figure 3.6d show the basic structures of the 2-hour and 1-hour assimilation configurations. For the 2-hr assimilation (Figure 3.5c), also for a forecast starting at t = 0 hr, observations are combined into 1-hr 64 interval, i.e., t = 1-3 hr. These observations are compared to forecasts at t = 1-3 hr, providing an analysis at t = 2 hr (each model run produces 3-hr forecasts). For the 1- hr assimilation (Figure 3.6d) starting at t = 0 hr, observations within t = 0.5 - 1.5 are simply combined together to be compared to the forecast at t = 1 hr, yielding an analysis at t = 1hr (each model run only generates a 1-hr forecast). Figure 3.6 Schematic illustration of the 6-hr 3D-LETKF, 6-hr 4D-LETKF, 2- hr 4D-LETKF and 1-hr 3D-LETKF Panel (a) and (b) are coming from Greybush et al. (2011). In 4D-LETKF, observations (combined together as red/pink line segments to discern between different time slots) are compared to forecasts (solid blue lines) at their correct hours rather than just to the 6-hr/1-hr forecast as in the 3D-LETKF. Figure 3.7 reveals the tide amplitudes for the 6-hr, 2-hr and 1-hr TES and MCS assimilation results, all using the short localization of 600km. We can easily see that the amplitudes of A-2,2 (SE2), are significantly reduced compared to that if the 6-hr 65 assimilation, and they even drop to a value below the free run. This could be further demonstrated in the analysis increments, both demonstrated in surface pressure (Figures 3.8a-c) and 10Pa temperature (Figures 3.8d-f): weak wave 4 patterns exist in both surface pressure analysis increment (Figure 3.8a) and 10Pa temperature analysis increment (Figure 3.8d) in the 6hr assimilation, but not in 2hr (Figures 3.8b,e) and 1hr (Figures 3.8c,f) assimilation. As a result, we could confirm that the 6-hr (4 times/sol) assimilation leads to an erroneous geographically fixed forcing that excites the near- resonant Kelvin wave. This problem does not exist in short window assimilations like 2-hr and 1-hr windows. Because the artificial separation of the analysis increments into 4 equally spaced groups disappears with short windows. Figure 3.7 Surface pressure semi-diurnal tide wave 2 amplitude for TES and MCS assimilations. (a) ~ (c) are from TES assimilations, (d) ~ (f) are from MCS assimilation. (a) and (d) are 6-hr LETKF assimilation, (b) and (e) are 6-hr 66 LETKF assimilation, and (c) and (f) are 1-hr LETKF assimilation. All are with ice cloud parameterization. Figure 3.8 7-day composite daily-average surface pressure bias (a)~(c) and level 7 (10Pa) temperature bias (d)~(e) over sol 397-403 for 6-hr (a), (d) 2-hr (b) (e) and 1-hr (c) (f) of TES assimilation. All are with ice cloud parameterization. 3.4.3 Impact of Water Ice Cloud Parameterization in Analyses and Forecasts Aerosol radiative forcing (dust and/or water ice cloud) has a strong influence on the depth-weighted 50-10 Pa semi-diurnal tide amplitude. There is evidence suggesting that clouds, not dust, provide the dominant radiative forcing in non-dusty seasons, such as the aphelion season studied in this paper (Wilson et al. 2008; Kleinböhl et al. 2013; Steele et al. 2014). Tropical clouds contribute to strong daytime heating by absorbing upwelling IR radiation from the relatively hot surface. 67 While clouds contribute to nighttime cooling by IR emission at cloud base, there is also significant adiabatic heating in the 100-10 Pa region by the tide-induced downward vertical velocity, resulting in a net heating effect in the tropics (Wilson et al. 2008; Madeleine et al. 2012; Wilson and Guzewich, 2014). Figure 3.9 shows both the free run and 1-hr LETKF assimilating MCS profiles, with and without the presence of Radiatively Active water/ice Clouds (RAC). It can be seen that for free runs, in the presence of RAC (100~10Pa), zonal average temperature semidiurnal tides amplitude show a significant signal in tropical regions (Figure 3.9c). However, in the absence of RAC, the free run shows a cold bias within this region (Figure 3.9a), suggesting that adding ice cloud yields a substantial difference in semidiurnal tide amplitude through the region where the RAC is present (100~10Pa). However, in the MCS LETKF 1-hr assimilation, the difference in the RAC region (100~10Pa) between with (Figure 3.9d) and without RAC (Figure 3.9b) is much smaller, which suggests that LETKF data assimilation can compensate for the absence of ice cloud parameterization. Therefore, the assimilation run does "correct" the cold tropical temperature bias in zonal mean temperature, as observed in the MACDA assimilation (Wilson et al., 2008; Montabone et al., 2011). At upper levels, since there are no direct temperature observations to constrain the assimilation, the forecast may contain large errors and therefore depart from the "true" state. Strong signals in upper level SH mid-latitude (Figures 3.9b-d) may be due to the vertical propagation of waves that transfers the mid-level RAC signal to upper levels, or may be the upper level response to mid-level RAC signals. 68 Figure 3.10 presents the semidiurnal tide amplitude of the temperature at 50Pa. As discussed in the previous section and shown in Figure 3.2a, the presence of a spurious wave 4 in the analysis increments with a 6hr assimilation window leads to the excitation of a mid-level temperature semidiurnal Kelvin wave (Figure 3.10b,f), and is removed by using short window assimilations (Figures 3.10c,d,g,h). Without the presence of RAC (Figures 3.10a-d), the westward propagating wavenumber 2, or A2,2, is reasonably enhanced in the 2-hr and 1-hr LETKF assimilating MCS profiles, while the free run and 6-hr assimilation may underestimate them (Figure 3 in Kleinboehl et al., 2013, indicates that at 50Pa A2,2 should be around 2K in the equator). The deficiency of the amplitude of A2,2 is largely due to the absence of RAC parameterization in the free run and 6-hr assimilation, since A2,2 in the free run and 6- hr assimilation are largely enhanced with RAC parameterization (Figures 3.10e,f). The amplitude differences of A2,2 with and without RAC parameterization in 2-hr and 1-hr assimilation (Figures 3.10c,d,g,h) are much smaller compared to the difference in the free run (Figures 3.10a,e). The difference in the 6-hr assimilation with and without RAC is smaller than that in the free run, but larger than that in the short window assimilations (Figures 3.10e,f). This may further indicate that besides removing spurious resonance, short assimilation windows further "push" the analysis closer to the true state even without the presence of necessary RAC parameterization. An amplitude increase in the short window length semidiurnal migrating tide at around 60°S is significant in scenarios both with and without RAC (Figures 3.10 c,d,g,h). The recent study (Figure 3 in Kleinboehl et al., 2013) suggests that the semidiurnal amplitude is ~2K in ~60°S. It is apparent that for both free run and 6-hr 69 assimilations, either with or without RAC, the amplitudes around 60°S are underestimated (Figures 3.10a,b,e,f) while short assimilation windows yield a more reasonable results that is closer to the estimated state in (Kleinboehl et al., 2013). (Figures 3.10c,d,g,h). This is not surprising since, as discussed before, shorter assimilation window lengths are favorable for EnKF schemes like LETKF, while long window lengths are optimal in 4D-Var assimilations. Since the twice-daily observations are insufficient to provide constraints on the semidiurnal tides, it may not be possible to state that the short window assimilations are closer to the true state, or determine whether 1-hr assimilation indeed produces closer results to the true state than 2-hr assimilation. In fact, the observations are currently not very constraining (Wilson and Guzewich, 2014), which makes the identification of the optimal assimilation window length harder. However, we did show that the spurious resonance, apparent both in surface pressure and mid-level temperatures induced by the 6-hr assimilation cycle, could be removed in the short window assimilations. Moreover, the results with both 2hr and 1hr assimilation windows are quite similar, suggesting that 2hr windows may be short enough to ensure linear growth. 70 Figure 3.9 Zonal average temperature semi-diurnal tide wave amplitude for free run and MCS 1-hr assimilation. (a) and (c) are from free run, w/o and with ice cloud parameterization, respectively; (b) and (d) are from TES assimilation, w/o and with ice cloud parameterization, respectively. Regions where ice cloud should be present (100~10Pa) are denoted between dashed lines. 71 Figure 3.10 50Pa temperature semi-diurnal tide wave 2 amplitude for free run and TES assimilations. (a) and (e) are from free run; (b)~(d) and (f)~(h) are from TES assimilations, 6-hr, 2-hr and 1-hr respectively. (a)~(d) are without ice cloud parameterization, and (e)~(f) are with ice cloud parameterization. 3.4.4 Impact of window length on the Martian atmosphere forecast skill We next examine the impact of assimilation window length on the skill of forecasts (compared to future observations) at various lead times. As mentioned before, adaptive inflation (Miyoshi, 2011) was used in the previous 6-hr, 2-hr and 1- hr experiments. In order to conduct “fair” comparison of forecast skill among experiments with different assimilation window length, initial adaptive inflation values should be scaled proportional to the assimilation window length (adaptive inflation reads inflation values from an initial file, and update the values during the assimilation process). In this study, for 6-hr LETKF assimilating TES profiles, covariance inflation value is initially 77% at model level 7 (the top level of TES observation), and decreases linearly to 0% at model level 4 and the levels above 72 where there are no observations. This configuration of inflation values is applied in the initial covariance inflation file. Since the background spread is inflated three times within 6 hours for 2-hr window assimilation, covariance inflation value is adjusted to 21% ( 3 1.77 =1.21). Similarly, covariance inflation value is adjusted to 10% ( 6 1.77 =1.1) for 1-hr assimilation window. In the current study, we only inflate the grid points within the localization radius of the observation, which means that if a model grid point is far away from any observations, there is no need to inflate this grid point. In order to estimate the effect of assimilation window length on the Martian atmosphere prediction, as well as study the impact of RAC on the prediction, two different periods are analyzed, one is Northern Hemisphere summer, near the Aphelion (~ Ls100) as in our previous tidal analysis, the other is NH autumn (~ Ls200). Although these two time periods are not symmetric about the sun (the Perihelion is ~ Ls 250), they do represent two typical time periods of an annual cycle. During the first time period (NH summer), the dominant variability is from the tides. There is relatively low dust loading during this season and model errors are dominated by the aerosol vertical distribution, in which water ice clouds likely play the dominant role (Wilson et al., 2008, 2014). Little to no traveling waves and dynamical instabilities are present during this time period. Figure 3.11 shows forecast root mean square differences with respect to temperature observation of 6-hr, 2-hr and 1-hr assimilation under different scenario, with and without RAC parameterization during NH summer. RAC demonstrates again its significance in the MGCM: adding RAC yields a substantial reduction in the RMS difference in the 73 freely running models. However, for the LETKF assimilations, RAC only slightly reduces the RMS difference (2%~5%). With the presence of the RAC, the free run RMS difference with observations is much closer to the assimilations, about 20% higher. This would further suggest the importance of the RAC parameterization, as indicated by Wilson et al. (2008) and Rogberg et al. (2010), that during the aphelion season the absence of a realistic representation of the radiative effects of the low latitudes water ice clouds may cause a significant systematic error in the Martian atmospheric prediction. This is confirmed by the tendency to have increasing errors right after the missing observations, especially in the absence of RAC, because during those times the "training" of the model stops and the model is “pushed back” towards the free run). 74 Figure 3.11 Performance of the MGCM-LETKF and the free run (without assimilating observation) during NH summer (~Ls 100), evaluated by comparing RMS differences of forecasts from TES observations. (a) is without ice cloud parameterization and (b) is with ice cloud parameterization. Time slots without observation are shaded in grey. Time-averaged RMS difference is denoted in the legend. Although 1-hr assimilation shows higher frequency fluctuations, its forecast RMSD with observations is the smallest. However, considering that 1-hr/2-hr assimilation RMS difference is comparing observations with 1-hr/2-hr forecast and 6- hr assimilation RMS difference is comparing observations with 6-hr forecast, longer forecast time always results in reduced prediction. 6-hr assimilation indeed shows larger RMSD than 1-hr and 2-hr assimilation. This is reasonable, but does not allow to judge which assimilation window length yields the best results. 75 In order to choose the most appropriate assimilation window, a comparison of prediction that requires the same forecast length is needed. 6-hr to 336-hr (9-day) forecast based on hour 00 and hour 12 forecast results during sol 393-407 are conducted, and compared to observations within ±3 hours (e.g., 6-hr forecast is compared to observations within hour 3-9). Since there are fluctuations in representativeness error depending on where the orbit passes, averaged RMSDs over the referred time period (both with and without RAC) are used, and shown in Figure 3.12. For the model runs with RAC parameterization, it is clear that although 6-hr assimilation benefits from more observations (at hour 00, 6-hr LETKF assimilates observations from hour 21 (previous day) to hour 03; 2-hr assimilation uses observations from hour 21 (previous day) to hour 01, and 1-hr assimilation utilizes observations from hour 21 (previous day) to hour 0.5), the prediction between different window lengths does not show large discrepancies prior to the “saturation” (~2 days). After having reached “saturation”, free run and assimilations should have RMSDs that are very close to each other, which means losing all "predictability". Repeating the experiments using different observation windows of ±0.5 and ±1 hours for verification showed similar results, and is not discussed in this study. 76 Figure 3.12 Performance of the MGCM-LETKF during NH summer (~Ls 100), evaluated by comparing RMS differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o and with ice cloud parameterization. For the model runs without RAC parameterization, it seems hard to judge which assimilation window performs best prior to the “saturation” since 6-hr assimilation performs best for short-range forecast (1-3 days) and short assimilation windows yield the smallest RMSDs for mid-range forecast (7-9 days). In the absence of the RAC, “saturation” occurs at day 9, which is postponed by a large extent from the RAC scenario. The low "predictability" (short time until error saturation) in the RAC scenario might be due to the near-perfect nature of the freely running model, suggesting that the RAC configuration in the model may be realistic, but not perfect considering the overshooting phenomenon. 77 During the NH fall (~Ls 200) when baroclinic instabilities associated with the synoptic waves are active, dust plays a more important role in the atmosphere while ice cloud role somewhat diminished. Figure 3.13 shows the same content as in Figure 10, but during NH fall. In this particular season, RAC plays little role in the model; adding RAC yields almost no difference in the RMS difference in both the freely running models and the LETKF assimilations. In the 1-9 days predictability shown in Figure 3.14, the RMS differences between with and w/o RAC parameterization are small, indicating that water ice cloud is not significant in this time period and its role could be neglected. During this time period, 6-hr assimilations yield the smallest RMSD while 1-hr assimilations show the worst, but the predictability difference between different assimilation window lengths is small, and we could think the benefit from utilizing more observations may make the 6-hr assimilation yield the optimal predictability considering all the assimilation window lengths, but assimilation window length does not play a significant role in affecting the Martian atmospheric predictability. From the above discussion, we may conclude that the short assimilation window lengths do not play a significant role in improving the Martian atmospheric predictability although they can remove the spurious resonance in the Kelvin wave component of the tides as well as improve the semidiurnal migrating tides in mid- level temperature simulation. Water ice cloud is an important factor that could affect the predictability in NH summer when aerosol loading is dominated by water ice clouds, but has little to no impact in NH fall when dust loading is high and there is less water ice cloud prevailing (synoptic waves and dynamical instabilities may 78 dominate during this period). When looking at the Martian atmospheric prediction, time to error saturation (RMSD from forecast matching free run) is an important index determined by both the radiative time scale, and the dynamical error doubling time. Figure 3.13 Performance of the MGCM-LETKF and the free run (without assimilating observation) during NH fall (~Ls 200), evaluated by comparing RMS differences of forecasts from TES observations. (a) is without ice cloud parameterization and (b) is with ice cloud parameterization. Time slots without observation are shaded in grey. Time-averaged RMS difference is denoted in the legend. 79 Figure 3.14 Performance of the MGCM-LETKF during NH fall (~Ls 200), evaluated by comparing RMS differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o and with ice cloud parameterization. 3.4.5 Bias Correction Since every model has its own deficiencies, including both “external” errors (e.g., sub-grid scale parameterization could never be perfect since it is hard to parameterize physical processes happening within the model grid) and “internal” errors (imperfect initial condition, uncertainties and model instabilities), bias always exists within the model. Model bias itself accounts much for the temperature analysis error; for example, a Mars model requires parameterization assumption for water ice cloud and dust properties, and inaccuracies within these parameterization processes will result in model bias form the true state. As a result, an empirical bias correction 80 method based on the idea of Danforth et al., (2007) is applied to improve the model prediction. One of the advantages of data assimilation is to help compensate for imperfect model parameterization by approximating the missing model forcings, treating the model as a black box and only considering forecast error statistics. These bias statistics can then be used to provide insights for improving model parameterizations and aerosol distributions, which are key goals for Mars science. The basic idea of empirical bias correction is to estimate model bias (in temperature, wind, surface pressure, and surface stress) as the time-mean analysis increment over a period of time, and then add the calculated bias to each of the ensemble member forecast before conducting the analysis. For short period assimilation, Greybush et al. (2012) used time-averaged analysis increments over sol 540-550 past perihelion and did the bias-corrected run over that period of time, and for a one-year reanalysis conducted a sliding-window bias correction using a 10-sol time window centering at the analysis date. In our experiments, we assume that there is no previously existing assimilation run, and the bias correction fields created from the analysis increments time averaged within the previous 7 sols are used to correct the forecast at each time step - we call this as "real-time" bias correction, which can save considerable computational time and resources. These results were compared to those using the bias correction method that relies upon a previous run, and their difference is only 1% in RMSD (figure not shown), which is small enough to justify using real-time bias correction since the speed is doubled. Considering 2-sol spin-up time from sol 388, we utilized the biases starting from sol 390, and started bias correction after day 391. If the date that is 7 sols prior to the current analysis time is prior to sol 390, then the 81 starting date to calculate the 7-sol average is hard-coded to sol 390, e.g., sol 391 uses daily-averaged bias of sol 390, and sol 392 uses daily-averaged bias averaged from sol 390 to 391. Since the model without RAC parameterization is known to have large deficiency in representing mid-level temperature, bias correction in the model that is far from perfect ought to have larger improvement. Greybush et al. (2012) improved empirical bias correction by treating separate bias correction fields for different times of day (hour 00, 06, 12, and 18, for traditional 6-hr reanalysis), which is termed as "diurnal" bias correction since different times of a day have different bias structures, and daily bias structures show an obvious diurnal cycle. Bias patterns at level 10 (near the top of TES observation levels) at different times of a day (hour 00, 06, 12 and 18) of 1-hr and 6-hr TES assimilation is shown in Figure 3.15, indicating that bias structure varies at different times of a day. Note that in the absence of RAC, the daily average bias (a - b) show significant positive patterns in 60°S-40°N region of both 1-hr and 6-hr analysis increments., and 6-hr daily average bias pattern shows a significant wave-4 pattern, further demonstrating its resemblance to the wave-4 topography forcing that scatters the wave 2, migrating semidiurnal tide to wave 2, eastward propagating Kelvin wave (Figure 3.16). This may suggest that diurnal bias correction - correcting bias separately at each correct hour should be needed. Figure 3.15 also shows the corresponding observation path at different times of a day (hour 00, 06, 12 and 18) over 7 sols (sol 397-403). It is also clear that, within 6-hr window there are more observation tracks, and the observation path transition at the same hour (say hour 00) between different sols just shows a smooth position shift (Figures 3.15e-h). However, 82 since the 1-hr observation paths are more sporadic, the observations at the same hour (say hour 00) of different sols contain not only a significant position shifting, but also a significant contrast in observation numbers (Figures 3.15a-d). Besides, the number of observations of different hours in the 1-hr observation path may vary much more than that in the 6-hr observation path. Due to the reasons presented above, larger sampling (representativeness) errors are expected in the 1-hr bias correction. Figure 3.15 7-day composite bias pattern (shaded colors) and observation track (dots) over sol 397-403 at level 10 for 1-hr TES assimilation (a)~(d), and 6- hr TES assimilation (e)~(h), at different times a day: hour 00 (a), (e); hour 06 (b), (f); hour 12 (c), (g); and hour 18 (d), (h). Different colors of dots represent observation path at different day. All are without ice cloud parameterization. 83 Figure 3.16 7-day composite daily-average bias over sol 397-403 at level 10 for 1-hr (a) and 6-hr (b) TES assimilation. Both are without ice cloud parameterization. Figure 3.17 Performance of the 1-hr MGCM-LETKF TES reanalysis during NH summer (~Ls 100), w/o and with bias correction (both sliding window daily average bias correction and sliding window diurnal bias correction) in a RMSD sense. Time slots without observation are shaded in grey. All are without ice cloud parameterization. Time-averaged RMS is denoted in the legend. Figure 3.17 clearly shows that in the absence of the RAC parameterization, the real-time sliding window bias correction runs have reduced the forecast RMSD compared to both no bias-correction run and the run with sliding window daily average bias correction. We choose the forecasts at hour 00 and 12 during sol 397- 407 from both 1-hr and 6-hr real-time sliding window diurnal bias correction runs, 84 conduct the prediction runs with bias correction at their assimilation times and compare forecast (1-9 days forecast) with observations within ±3 hours, as described in Section 3.4.4. The bias used in the prediction bias-correction experiment is the sliding window 7-sol composite prior to the forecast time. Figure 3.18 Performance of the MGCM-LETKF during NH summer (~Ls 100), evaluated by comparing RMS differences of forecasts based on hour 00 and hour 12 forecasts from TES assimilations with TES observations, w/o bias correction, with sliding window bias correction and with sliding window diurnal bias correction. All are without ice cloud parameterization. (a) is TES assimilation, (b) is TES assimilation only considering 30S-30N, (c) is MCS assimilation. The bias-corrected predictions show a big part of RMSD is persistent temperature bias in this season (Figure 3.18). Figure 3.18a shows the prediction of the TES 85 assimilation. It is clear that for 6-hr bias correction, using diurnal bias yields better results than using daily-averaged bias. However, for 1-hr bias correction, situation is reversed. The reason why 1-hr diurnal bias correction does not work well could be explained by the observation increment in the 120-hr forecast with bias correction (Figure 3.19c), where unrealistic signal is amplified in upper level, SH polar front region. This could be explained by the sensitivity of baroclinic waves to large sampling errors in 1-hr diurnal bias. We then limit the analysis region in the tropics (30°S-30°N) where baroclinic waves are not prevailing, and then the diurnal bias works better for both 1-hr and 6-hr corrections (Figure 3.18b). It is also apparent that 1-hr diurnal bias correction reduces the bias in the observation increment in the upper and middle atmosphere in the tropics (Figure 3.19c) compared to the 1-hr daily- average bias correction (Figure 3.19b). Another possible reason that could explain the failure of 1-hr diurnal bias correction in testing the prediction is due to error in the TES upper boundary (~10Pa). Note that the observation increments in the bias- corrected runs without the presence of RAC (Figures 3.19b,c,f,g) are also smaller than in the original runs without bias correction but with RAC parameterization (Figures 3.19d,h), suggesting that bias correction could also well compensate the absence of water ice cloud. Forecast skill test with bias correction is also conducted in MCS assimilations to overcome this problem as the coverage of the MCS data is the whole atmosphere. It is evident that for MCS experiments, diurnal bias correction in the forecast always works better than the daily-average bias correction (Figure 3.18c, Figures 3.20c,f). This may further demonstrate the advantage of using diurnal bias correction. The 86 prediction in the MCS runs without bias correction without the presence of RAC (3-5 sol) is different from that in the TES assimilation runs (7-9 sol), suggesting that in a highly imperfect model without RAC, the prediction may depend on the chosen of the observing system, again demonstrating the importance of the RAC parameterization. Figure 3.19 Bias pattern of 120 hour forecast in 1-hr TES LETKF (a)-(d) and 6-hr TES LETKF (e)-(h) runs. Left panels are original forecast without bias correction; second column panels correct the sliding window daily-average bias every 1hr/6hr; third column panels correct the sliding window diurnal bias every 1hr/6hr; the above-mentioned panels are all without RAC parameterization; right panels are original forecast without bias correction, but with RAC parameterization. What must also be noticed here is correcting the bias every 6 hour shows some advantages over correcting every hour (Figures 3.18a-c). Both TES and MCS 120-hr forecast with bias correction show that 1-hr bias corrections (Figures 3.19b,c; Figures 3.20b,c) do not work as well as 6-hr bias corrections (Figures 3.19f,g; Figures 3.20e,f) in SH polar front regions, suggesting that frontal dynamics may not work well for the 1-hr bias correction. Besides, 1-hr bias correction may induce larger 87 sampling errors than the 6-hr bias correction as shown in Figure 3.14, which may also decrease the forecast skill. Figure 3.20 The same as Figure 18, but for MCS assimilations. 3.4.6 TES vs. MCS Our current work also gives an overview of assimilating MCS profiles using MGCM with the 4D-LETKF framework. Since MCS temperature profiles have greater vertical resolution than TES which smooths out the temperature structure, they are expected to result in a more detailed longitudinal and vertical tide signal (Figure 3.21, Figure 3.22). Although MCS has a better representation in equatorial region vertical structures than TES does (TES data yields a fairly uniform lapse rate through much of the atmosphere), it may be problematic in high latitude when dealing with strong temperature gradient, especially in polar vortex regions. Although MCS 88 profiles extend to 80km, upper level observations contain large error, which significantly influence the assimilation results. Model forecast in MCS assimilation yields even much larger error in upper levels compared to observation error, producing huge spread and large covariance inflation values in upper levels. This makes forecast RMSD in MCS assimilations much larger than those of TES assimilations (Figure 3.18). Even if only considering TES levels (~40km), MCS assimilations still have larger RMSDs than TES assimilations (Figure is not shown), this might be attributed to MCS large errors in polar front regions. Another possible explanation is because the model finds it more challenging to match the more complex vertical structure in MCS observations. Figure 3.21 MCS (a) and TES (b) reanalysis of equatorial temperature difference (3pm-3am). Both are with ice cloud parameterization. 89 Figure 3.22 Free run (a), TES (b) and MCS (c) reanalysis of zonal mean temperature diurnal tide amplitude. All are with ice cloud parameterization. 3.5 Summary and Discussion We examine the diurnal cycle of Mars atmospheric analyses created using the LETKF, MGCM, and TES/MCS spacecraft observations. Features show a large sensitivity to assimilation phase and moderate sensitivity to localization radius, representing the non-uniqueness of the 6-hr data assimilation system, which is unwanted. The use of a “traditional” 6-hr assimilation window creates wave 4 patterns in the daily averaged analysis increment, serving as an external forcing to trigger spurious resonance, which is similar to how a wave 4 pattern in topography modulation excites an eastward propagating semidiurnal east wave. Harlim and Majda (2008) also noted a lack of observability when assimilating regularly-spaced, sparse observations with a resonant period in a simplified chaotic model; the data assimilation system cannot discern the true signal from among the various Fourier harmonics in the aliasing set based on observations alone. 90 In order to test whether this resonance mechanism is responsible for the large amplitude semidiurnal tidal wave observed in our reanalysis, we recomputed the reanalysis using different assimilation window lengths. As expected, a 12-hr assimilation window led to wave 2 pattern in the analysis increments. Short assimilation windows such as 2-hr and 1-hr not only eliminated this artificial resonance induced by the traditional 6-hr assimilation, but also reduced the RMS difference between the forecast and the observations, and better compensated for the absence of the radiatively active water ice cloud parameterization. It is not surprising that the short assimilation windows show an advantage over the 6-hr window within LETKF. Ensemble Kalman Filter assumes that the ensemble forecast perturbations with respect to the mean and the observational errors are both Gaussian, so that when they are combined using the Kalman Filter formulation, the analysis ensemble perturbations are also Gaussian (Hunt et al., 2007). This implies that short windows are advantageous for EnKF. In our experiments with the assimilation of TES and MCS temperature profiles on the Mars atmosphere we have found that assimilation windows of 1 hour and 2 hours succeed in eliminating the resonance observed when using a 6 hour window, and their results are similar except that one hour analyses tend to have slightly larger RMS differences with the observations. This suggests that for very short windows, since there are fewer observations per cycle, the analysis increments may be more sensitive to noise from observation error and sampling error in the background ensemble error covariance matrix, and that the optimal EnKF should be short enough to ensure linear growth, but not too short in order to ensure enough observations are processed 91 simultaneously. By contrast, 4D-Var data assimilation is optimal for very long windows, in which the assimilation is exposed to the maximum number possible of noisy observations (Pires et al., 1996). They found, however, that for very long windows, due to the presence of multiple minima in the variational cost function, it was necessary to introduce a Quasi-static Variational Assimilation (QVA) in which the assimilation windows were started small, in order to reach the true single minimum, and then increased in small increments to reach the optimal long window. This was observed by Singleton (2011), who showed that for a toy coupled ocean- atmosphere model, the optimal coupled variational assimilation was reached using very long windows and the QVA algorithm. In contrast, the optimal coupled LETKF analyses were obtained with very short windows, but for optimal windows, the analyses errors were similar for the EnKF and the 4D-Var approaches. Although all the assimilations showed significant improvements in forecast RMSD with respect to independent observations compared to a freely running model, short window assimilations further reduced the RMSD compared to the 6-hr assimilation. Considering that short window assimilations use short-term forecast to compare to nearby observations, this provides an advantage to short windows so that it may not be fair to compare among different assimilation forecasts with different assimilation window length. Therefore, we also launched 14 day forecasts during both the NH summer (~Ls100) and the NH fall (~Ls200), starting from hour 00 and 12 initial conditions, and calculate the averaged forecast RMSD compared to all the observations within ±0.5, 1, and 3 hour of the forecast time. In this longer forecasts, the assimilation window length appears to have only a minor impact on the forecast 92 skill. Instead, the model forecast skill highly depends on the role of water ice clouds (RAC) during NH summer when water ice cloud is prevalent. During NH fall, the forecast skill may be dominated by other effects such as synoptic waves and dynamical instabilities. The ultimate goal of this study is to provide the basis for a multi-year, short assimilation window reanalysis that contains diurnal, seasonal and interannual variability, removes spurious resonance, as well as gets closer to the possible “true” state of the Martian atmosphere. We are interested in further improving the model forecast, and therefore develop a real-time sliding window bias correction method, which does not rely on previous reanalyses. Since bias structure is different at different times of a day due to the changing satellite track, “diurnal” bias correction, which treats bias correction fields differently for different times of a day, is developed to further reduce forecast RMSD as well as enhance prediction. It has been demonstrated that conducting diurnal bias correction in the model forecast improves the prediction compared to conducting the daily-average bias correction, and correcting the model bias every hour in the forecast may induce larger sampling error than correcting it using longer window length. The bias corrections are also demonstrated to have the ability of compensating for the absence of water ice clouds. Although our ultimate goal is to develop a multi-year reanalysis that could best describe the Martian atmospheric state, it is still hard to determine the whether 1-hr or 2-hr is the optimal window length although we have already determined that short window lengths have a better performance than the traditional 6-hr assimilation. The sparseness of the observations makes the determination of the optimal assimilation 93 window length a hard task. Although the Viking program provides surface pressure time series that are relevant to the TES and MCS years in this particular season (NH summer), the observation data are missing at Viking Lander 2 around ~Ls100, only data at Viking Lander 1 are available for comparison. The wave number 2/4 pattern is the interference between two dominant wave modes of diurnal/semidiurnal tides, respectively, and therefore are sensitive to assimilation window length. As a result, the comparison between the simulated diurnal and semidiurnal tides and the tides provided by VL1 only provides weak evidence to determine the optimal assimilation window length since the “truth” provided by the Viking Lander sites is far from enough. There are several measurements of diurnal temperature variations, such as measurements made during the MGS TES aerobraking and science phasing period when the orbit drifted with respect to local solar time during which good sampling of the diurnal varying thermal structure is allowed (Banfield et al., 2000); and during the MGS mapping period when two local times were sampled because of the sun- synchronous orbit that results in aliasing of higher order modes into the diurnal response (Banfield et al., 2003). During the aerobraking and science phasing period, semidiurnal tides were also estimated, however, there are still not enough evidence for providing a reliable “truth” since data are limited in the southern hemisphere, especially south of 20°S (Banfield et al., 2000). Harmonic contributions such as diurnal and semidiurnal wave could be calculated by tides fitting only if the data is available at least six local times per day. Hence, errors from fitting the sparse observations could be expected. One major source of the error is the departure of the 94 fitted quantity from the fit without perturbation. Besides that, there are uncertainties accompanied with the uneven data coverage in local time at low latitudes, which make the “true” tides estimation, especially the semidiurnal tides estimation an even harder task (5K in Kleinbohl et al., 2013). The sparseness of the observation, as well as the errors and uncertainties with the tides estimation, introduces large ambiguity to the “truth” of the tides, and thus makes the estimation of the optimal window length an uneasy task. Our finding of assimilation-induced artifacts in atmospheric motions may have relevance to the Earth as well. Mitchell et al. (2014) highlights the dynamical similarities between the Martian atmosphere and the Earth’s stratosphere. Tidal motions are important in the Earth’s upper atmosphere as well (Coy et al., 2011; Zakasaki et al., 2012; Hagan and Roble, 2012). Swinbeck et al. (1999) noted a contrast in tidal amplitudes between simulations and assimilations for the Earth’s stratosphere, suggesting possible improvements that should be made to the satellite retrievals. With the Earth’s upper atmosphere being an increasing area of research for assimilation studies (e.g. Polavarapu et al., 2006; Sankey et al., 2007; Hoppel et al., 2013), the impact of satellite orbits and assimilation windows should be considered. Since the diurnal period of Kelvin wave resonance is special to Mars, and the comparable Kelvin resonance period on Earth is 33 hours, it does not have the same potential of forcing resonant atmospheric modes, as discussed in Section 2.4. Note that atmospheric tides on the Earth are analogous to Martian tides in many ways. Like Martian tides exist in all levels mainly due to absorption of solar heating by aerosols, Earth tides exist in both troposphere and stratosphere, as water vapor in the 95 troposphere (~0–15 km) and ozone in the stratosphere (~30 to 60 km) absorb solar radiation during the day. The former one absorbs longwave radiation accounting for 20% of the semidiurnal tides amplitude while the latter one absorbs shortwave radiation accounting for 70% of the semidiurnal tide amplitude. The rest 10% of the semidiurnal amplitude could be explained by latent heat related with tropical convection (Lindzen and Chapman, 1969; Lindzen, 1978). They both contain migrating (sun-synchronous) and non-migrating tides originating from topography, surface interactions and latent heat release due to deep convection in the tropics. Diurnal waves originate from direct solar heating with a period of 24 hour, while semidiurnal tide and other higher frequency waves, i.e., with periods of 8hr or 6hrs, are mainly due to the decomposition of the approximate square wave profile of solar heating of the atmosphere, which is rich in harmonics, into separate frequency components using a Fourier transform. However, tides on the Earth have a much smaller amplitude. For example, the averaged surface pressure semidiurnal tide amplitude in Mars is around 1% throughout the year, and it is one order of magnitude smaller in the Earth (Lindzen and Chapman, 1969). Another major difference between Martian tides and Earth tides is that the diurnal tide dominates in the Martian tropics in a clear atmosphere, and the semidiurnal tide only becomes apparent with increasing dust loading. By contrast, semidiurnal tides prevail on Earth’s surface pressure due to the fact that vertical propagation of inertia-gravity wave becomes evanescent for diurnal tides (Lindzen, 1990). In summary, atmospheric tides in Mars and the Earth, although driven by the diurnal heating cycle, have different forcing mechanisms and dominant periods. 96 Nevertheless, the presence in Mars of a clear spurious resonance associated with the use of an assimilation cycle of 6hrs (1/4 of a day), suggests that similar (but weaker) effects may be present in the Earth analysis. 97 Chapter 4: Summary and Lessons Learned The breeding method, proposed by Prof. Kalnay, is a nonlinear, finite time generalization of the method to compute leading Lyapunov Vectors that has been effectively used in different areas such as data assimilation, weather and climate forecasting and estimation of growing errors in nonlinear systems. The study in Chapter 2 is the first time that bred vectors have been used to identify characteristics of consistent instability in ocean wave interactions and forecast them. The main result we have found, that Bred Vector growth precedes the growth of the rogue waves in the NonLinear Schrödinger Equation, is very encouraging and indicates that the use of breeding should lead to prediction of rogue waves in full nonlinear equations. As precursors of energy localization, bred vectors may bring potential societal and economical benefits, especially for global shipping and offshore energy harnessing. The results also suggest that in a complete numerical solution, breeding would allow to establish the transfers of energy leading to rogue waves. This project is a joint project incorporating efforts from different members from Department of Mechanical Engineering, Department of Astronomy and Department of Atmospheric and Oceanic Science, and is supported by NSF Grant No. CMMI-1125285. The paper “Breeding Analysis of Growth and Decay in Nonlinear Waves” authored by Yongjing Zhao, Eugenia Kalnay, Christopher Chabalko and Balakumar Balachandran was submitted for publications in Physica D on June 3rd 2014. Since the breeding method requires numerical solutions, this project involves efforts of obtaining numerical solutions of partial differential equation (PDE)s, such 98 as KdV equation and NLSE. Developing and applying numerical schemes to obtain the numerical solutions to PDEs became the most challenging part. For both KdV equation and NLSE, the Fourier spectral/exponential time-differencing fourth-order Runge–Kutta (ETDRK4) scheme is used due to its numerical stability. In addition to the ETDRK4 scheme, Strang-Splitting scheme is also used for the NLSE and the detection of numerical instability by bred vectors is demonstrated since the Strang- Splitting scheme exhibits numerical instability after a period of time. The original purpose for studying the application of the ETDRK4 scheme is for the Dysthe equation, since the direct solution to the NLSE is the envelope solution without providing any carrier waves, and it is hard to see the possible energy transfer through the envelope solution using bred vectors. The Dysthe equation (Islas and Schober, 2011) contains both KdV and NLSE terms plus damping terms, and the breeding will be applied to study the energy transfer in a full field rogue wave solution, rather than the envelope solutions shown here. The study of the Dysthe equation is part of the most important future works. Thermal tides are a prominent feature of Martian atmosphere; representing thermal tides correctly therefore becomes a very important factor when evaluating the Martian atmosphere reanalysis. Greybush et al. (2012) assimilated TES temperature profiles into the GFDL Mars Global Climate Model, created a full-year reanalysis consisting of temperature, wind, and surface pressure, extending from autumn of MY 24 to autumn of MY 25. However, after studying the semidiurnal Kelvin wave, we found that the “traditional” 6-hr reanalysis used before (Greybush et al. (2012)) introduces spurious resonance by forming a wave-4 analysis increment acting as 99 topographic forcing. Short assimilation window length such as 1-hr and 2-hr assimilation could remove this spurious resonance, as well as improve the mid-level semidiurnal migrating tides representation. The study in Chapter 3 is the continuous study following Greybush et al. (2012)’s study. The assimilation of MCS profile is introduced, and the results of assimilating MCS profiles are expected to be published in the future. Forecast skill studies assessing the short to long range forecasts (1-15 sols) initiated from this reanalysis can help evaluate the quality of the MGCM model as well as the 4D-LETKF data assimilation system, and therefore provide a better understanding of the behavior of the reanalysis during observation data gaps. This Mars project is also a team-collaborating project, including efforts of team members from University of Maryland (UMD), the Pennsylvania State University (PSU), Geophysical Fluid Dynamics Laboratory (GFDL) and Atmospheric and Environmental Research (AER). This work was supported by the NASA MDAP and PATM programs, including NASA grants NNX07AM97G and NNX11AL25G. Due to these advantages of short assimilation window length, we re-generated a 3-yr Ensemble Mars Atmosphere Reanalysis System (EMARS) reanalysis with MCD dust using 1-hr assimilation window length. Other publications regarding the 3-yr EMARS reanalysis are also anticipated. 100 Appendix A: The Indian Ocean Dipole: A Monopole in SST Published in Journal of Climate, doi: 10.1175/JCLI-D-14-00047.1 with coauthor S. Nigam. A.1 Introduction The tropical Indian Ocean (IO) basin – home to pronounced seasonal low-level wind variability including direction-reversal (monsoonal flow) – exhibits notably weak interannual variability in SST and surface-winds (e.g., Nigam and Shen, 1993), especially in comparison with the Pacific where interannual El Nino Southern Oscillation (ENSO) variability is impressive and influential. The proximity of the two basins and ENSO’s large-scale structure and near-global response provides scope for basin interaction. ENSO, in fact, does influence the Indian Ocean SST and low- level winds, especially during July-November (Nigam and Shen, 1993, cf. Figure 2). More recently, Saji et al. (1999) identified a dipole pattern of SST variability in the tropical Indian Ocean which is widely referred as the Indian Ocean Dipole (IOD). The identification was based on the EOF analysis of SST anomalies of all calendar months in the 1958-1998 period. The leading pattern, describing basin-scale anomalies of uniform polarity, was taken to represent the ENSO influence on the Indian Ocean, and the second pattern – the IOD – was thus considered temporally independent of ENSO. This basic premise of Saji et al. is questioned in this study. Saji et al.’s characterization of their leading EOF pattern as representing ENSO’s influence on the Indian Ocean is reexamined for the following reasons. First, ENSO’s influence varies with season whereas their leading EOF is seasonally invariant: The ENSO-related warming of the tropical Indian Ocean, for example, is 101 focused in the northwestern basin in summer-fall and in its southeastern sector in boreal winter (Nigam and Shen, 1993, Figures 2, 4). Second, characterizing (filtering) ENSO’s influence is challenging as its evolution is not representable by a single index or EOF due to complex spatiotemporal development (e.g., Guan and Nigam, 2008; Compo and Sardeshmukh, 2010). Saji et al.’s characterization of the leading EOFs thus requires reconsideration. Since Saji et al.’s analysis, the IOD structure and impacts have been widely analyzed using the Dipole Mode Index (DMI; Saji et al., 1999). The DMI index, interestingly, is neither tightly correlated with the second EOF’s time series (~0.7, or only 50% common variance) nor independent of ENSO (correlation with Nino3.0 SST index is ~0.35); all as reported in Saji et al. (1999). The correlation of the EOF time series itself with the Nino3.0 (or related) SST index was not reported by the authors, leaving open the question of IOD’s independence from ENSO. Not surprisingly, the IOD–ENSO link has been extensively debated: Several observational studies lend support to IOD’s independence from ENSO (e.g., Behera et al., 1999; Webster et al., 1999; Murtugudde et al., 2000; Ashok et al., 2003; Behera et al., 2003; Yamagata et al., 2003), while others question the same (e.g., Chambers et al., 1999; Reason et al., 2000; Allan et al., 2001; Nicholls and Drosdowsky, 2001; Baquero-Bernal et al., 2002; Dommenget and Latif, 2002; Hastenrath, 2002; Xie et al., 2002; Dommenget et al., 2007; Jansen et al., 2009; Dommenget, 2011). Allan et al. (2001) claim that ENSO’s spatiotemporal evolution is aliased in the DMI index but these authors did not question the basis of the index itself. Meyers et al. (2007) could 102 not infer a clear link between IOD and ENSO events, attributing the lack of clarity to the decadal variations in thermocline depth off Java-Sumatra. Modeling studies (e.g., Baquero-Bernal et al., 2002; Shinoda et al., 2004; Yu and Lau, 2005; Kug and Kang, 2006; Jansen et al., 2009) also provide insights on the contribution of ENSO and ocean dynamics in generating Indian Ocean SST variability. From analysis of the tropical Indian Ocean response to observed wind forcing, Shinoda et al. found the dipole mode to be the leading mode of thermocline depth variability, with the ENSO-independent dipole variability more strongly expressed in the upper-ocean heat content than SSTs. The successive IOD events of 2006 (El Nino year) and 2007 (La Nina year) spurred further interest in IOD genesis and predictability (e.g., Luo et al., 2008), including the possibility of IOD’s link with non-canonical ENSO variability (Ashok et al., 2007; Luo et al., 2010); the latter is referred as the Non-Canonical ENSO mode by Guan and Nigam (2008) and El Nino Modoki by Ashok et al. (2007). The evidence for ENSO’s influence on IOD is growing: Yamagata et al. (2004) found one-third of IOD events connected with ENSO using seasonally-stratified correlations. Behera et al. (2006) found a significant fraction of IOD events correlated with tropical Pacific variability, including ENSO, in their modeling study. The influence of IOD on ENSO has also been investigated (e.g., Wu and Kirtman, 2004; Dommenget et al., 2006; Jansen et al., 2009; Izumo et al., 2010; Frauen and Dommenget, 2012), with indications of inter-basin interaction on both interannual and decadal timescales (Huang and Shukla, 2007a,b; Dommenget, 2011). 103 Given the large body of literature on analysis of the IOD–ENSO link, it is, perhaps, necessary to articulate the goals of this observational study:  Most observational investigations of the IOD–ENSO link begin with IOD’s putative characterization – a zonal dipole structure – manifest in the DMI index, and seek index refinement by factoring for ENSO’s influence; that is, they seek a posteriori adjustments. This study, in contrast, questions the IOD’s canonical characterization itself, i.e., the basis for the DMI index. It thus seeks an a priori accounting of ENSO’s influence, like Meyers et al. (2007).  ENSO’s influence on Indian Ocean SSTs is characterized taking into account ENSO’s complex development: The spatiotemporally varying impact of both canonical and non-canonical ENSO variability is estimated and filtered prior to the search for recurrent modes of SST variability in the Indian Ocean. Many previous studies have estimated this influence – incompletely, in our opinion – from the tracking of ENSO’s mature phase alone and related compositing (e.g., Yamagata et al., 2004).  ENSO-filtered Indian Ocean SSTs are analyzed using the extended empirical orthogonal function technique (extended-EOF) which focuses on spatial and temporal recurrence, and as such, yields modes (rather than patterns) of variability – all rooted in the Indian Ocean basin, in this case. Unlike some earlier studies (e.g., Behera et al., 2003), SSTs are not smoothed or filtered in any manner, except for ENSO variability, avoiding potential aliasing of the 104 SST record. A century-long SST record is also analyzed here in the interest of robust findings. The ENSO characterization and the follow-on Indian Ocean (IO) SST analysis are based on the recent innovative analysis of natural variability and secular trend in the Pacific (and Atlantic) SSTs in the 20th century (Guan and Nigam, 2008, hereafter GN2008). By focusing on spatial and temporal recurrence, the extended-EOF analysis discriminates between interannual and decadal-multidecadal variability, and the nonstationary secular trend – all without any advance filtering (and potential aliasing) of the SST record. The Atlantic SSTs were similarly analyzed but after excluding the influence of Pacific SSTs and the SST secular trend on the Atlantic basin (Guan and Nigam, 2009); leading to a clarified view of the Atlantic Multidecadal Oscillation (AMO) structure and the implicit decadal time-scale exchanges of sub-Arctic and North Atlantic water (e.g., the 1980s Great Salinity Anomaly, Guan and Nigam, 2009). The Atlantic basin analysis serves as a prototype for this Indian Ocean SST analysis. The related Pacific and Atlantic PC analysis is introduced in section A.2, and the related data used in the Indian Ocean analysis is introduced in section A.3. The claim for the dipole structure of interannual SST variability in the tropical Indian Ocean is examined in section A.4. Section A.5 examines if the 1994 and 1997 SST anomalies indeed represent IOD events, as suggested by Saji et al. (1999, Figure 1). The search for dipole variability is extended to the Indian Ocean subsurface in Section A.6, which examines the structure of upper-ocean heat content variations linked with the eastern pole of the IOD (the more viable of the two poles). The IO 105 SST analysis is described in section A.7 and concluding remarks follow in section A.8. A.2 Introduction to SST Principal Components A.2.1 Introduction As Sea Surface Temperature (SST) is identified to play a significant role in global climate, large attention has been paid to large-scale SST departures, or anomalies manifest on interannual, decadal, and multidecadal to century-long time scales in the Pacific and Atlantic basin. Guan and Nigam (Guan and Nigam, 2008) explored 7 Pacific spatiotemporal variability modes by applying the extended EOF (EEOF) analysis (Weare and Nasstrom, 1982), which can effectively characterize variability evolution, including nascent phase structure. Same analysis is conducted to analyze Atlantic and 4 spatiotemporal variability modes are explored (Guan and Nigam, 2009). However, their analysis is until 2002. Since we have the latest datasets, it is necessary to expand the time period to the most recent years. A.2.2 Data and Analysis The data used in this chapter are the Met Office’s (UKMO) Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) 1.1 SST, globally available on a 1° x 1° grid for the period from 1870 onward, which were then constructed by applying a reduced-space optimal interpolation and subsequent insertion of quality- controlled gridded observations onto the reconstruction (Rayner et al., 2003). 106 HadISST is chosen as the dataset for analysis of Pacific and Atlantic low-frequency variability as its relatively long record, data continuity, and statistical homogeneity. Seasonal means are then calculated based on the monthly values, and the 1° x 1° grid are interpolated to a 5° x 2.5° longitude-latitude grid in order to meet with the computational efficiency. What must be mentioned here is that gridded SSTs were weighted by the square root of the cosine of the latitude (e.g., Chung and Nigam, 1999) in order to guarantee each grid cell’s influence in the analysis is commensurate with its area. Therefore, the SST data were not subjected to any additional normalization, especially by the local standard deviation, and then the undertaken EOF analysis is covariance based. Extended EOF, a variant of regular EOF (Guan and Nigam, 2008, 2009) is applied to conduct the following analysis. For Pacific Ocean, the analysis is conducted in the Pan-Pacific domain (20°S–60°N, 120°E–60°W), and the time period is 1900-2011. A five-member anomaly sequence at seasonal intervals is used, and 7 modes are rotated. Prior to the EEOF analysis, residual SST anomalies are calculated by removing the PC regressions of the six Pacific natural variability modes and the secular trend mode from the seasonal anomalies, which were then area weighted, but not normalized. The analysis of residual Atlantic SST anomalies is then conducted within 20°S–60°N, 80°W-20°E. A seven-member anomaly sequence at seasonal intervals is used, and 4 modes are rotated. As the Matlab EEOF program only requires whole years, so the time period for Atlantic analysis is 1901-2008. Note that the slightly different start and end time points arise from the use of extended-EOF 107 (EEOF) analysis technique which results in the principal component (PC) losing a few time points at both the beginning and end of the analysis period, on account of the temporally-extended sampling window. Filtering of the Pacific (and Atlantic) SST influence, and our analysis code quirk requiring whole years leads to the slight end- point differences. For reference, the Pacific analysis period is 1900-2009, while the Atlantic one (after removal of Pacific’s influence; losing two seasons at the beginning and end due to the 5-band analysis in Pacific) is 1901-2008, and loses three more seasons at the beginning and end due to the 7-band analysis in Atlantic. A.2.3 Pacific PC Analysis Figure A.1 lists seven leading principal components of Pacific SST variability. The first two leading modes describe the decaying phase (ENSO+) and growing phase (ENSO-) of ENSO, which explain 18.5% and 15.0% of the explained variance respectively, together they depict the canonical ENSO (cf. Rasmusson and Carpenter, 1982). The fifth model which exhibits different spatial and temporal evolution pattern as canonical ENSO is referred as the noncanonical ENSO mode, which accounts for 4.3% of the explained variance, also as El Nino Modoki by Ashok et al. (2007) and similar as the Pacific meridional mode (Chiang and Vimont, 2004). In the post-climate shift (1977-onward), noncanonical ENSO begins to show longer time scale and stronger amplitude, which coincides with previous studies on increased ENSO energy and persistence in the post–climate shift period (Gu and Philander, 1995; Zhang and Busalacchi, 2005). Biennial variability, the seventh leading mode of Pacific spatiotemporal variability (3.7% of the explained variance), described by Rasmusson et al. (1990) also serves as an important component which affects Indian 108 Ocean Dipole variability (See Section A.4 – A.6). Characterized by the near 2-yr time scale of individual warm and cold events and their phase locking with the annual cycle, Rasmusson et al. (1990) claimed for its importance in the ENSO cycle, who viewed ENSO as a superposition of biennial variability and lower-frequency variability. 109 Figure A.1: Principal components (PCs) of the seven leading modes of recurrent spatiotemporal SST variability in the Pacific (shaded) during 1900- 2009. Tick marks on the vertical axis denote three standard deviations. Smoothed PCs after 50 applications of a 1-2-1 smoother are represented in solid black lines. 110 Figure A.2: Spatiotemporal evolution of recurrent SST variability in the Pacific: Seasonal lead-lag regressions of the leading principal components (PC 1- 7). 111 The fourth and sixth mode of pacific SST variability capture decadal variability, say Pan Pacific (fourth) (5.8% of the explained variance) and North Pacific (sixth) (6.0% of the explained variance) decadal variability, show variations on long time scales (longer than ENSO-related variability), and capture the major climate shifts in the twentieth century. The 1920 shift is captured in PC4 (PDVPP) and the 1940s and mid-1970s are manifest in PC6 (PDVNP). These shifts are also discussed by Mantua et al. (1997) and Minobe (1997). The secular trend of the Pacific SST record is captured in the third mode (11.6% of the explained variance). This trend mode is much similar as both the global surface air temperature trend (e.g., National Aeronautics and Space Administration Goddard Institute for Space Studies (NASA GISS)) and the ocean heat content trend (Levitus et al., 2005). Figure A.2 depicts the spatiotemporal evolution of all the seven Pacific SST variability modes. During the growing phase of ENSO (ENSO-), positive SST anomalies appear in the eastern equatorial Pacific at the beginning and then expand westward until they reach and across the date line. In comparison, the decaying phase, say ENSO+, is characterized by the loss of amplitude across the basin. The above are the development process of canonical ENSO. For noncanonical ENSO, the spatiotemporal evolution is characterized by the eastward development of positive SST anomalies along the equator which begins at the central Pacific. Monahan (2001) found a noncanonical ENSO mode that also characterized temporal variability which is relative to the 1976/77 climate shift by using nonlinear, neural network–based EOF 112 analysis, however, spatial evolution of this mode differs from the one we talked about here. The mature phase of PDVPP is featured by small region of negative SST anomalies in the central North Pacific, and stronger positive anomalies outside which expand eastward from the Bering Sea and downward until Baja California, and finally go back southwestward toward the tropical Pacific. What must be noticed here is the lack of SST anomalies in the central and eastern equatorial Pacific. Moreover, it is also featured by the clockwise development along the North American west coast. Similar Pan-Pacific modes are also discussed in several other studies (e.g. Barlow et al., 2001). North Pacific mode (PDVNP) is featured by a zonally elongated band of positive SST anomalies which are centered along the sub-Arctic front and the date line. Spatiotemporal evolution process of PDVNP can be described as follows: cold SSTs begin to develop in the eastern basin while warm SSTs begin to develop in the central North Pacific. Similar North pacific modes are also discussed in several other studies (Deser and Blackmon, 1995; Nakamura et al., 1997; Mantua et al., 1997; Nigam et al., 1999; Barlow et al., 2001; Wu et al., 2003). PDVNP plays a strikingly significant role in the rainfall reconstruction in China, which is part of the undergoing work. Modal structure of secular trend mode shows that Northern Hemisphere is dominated by warm SST except for two small regions of cooling: an equatorial sliver in the central Pacific and another off the southern tip of Greenland, with largest warming appearing along the eastern coasts of Asia and North America. It should be 113 noted that another cooling center in the Northern Hemisphere is found in previous studies (e.g. Cane et al., 1997; Kaplan et al., 2000). A.2.4 Atlantic PC Analysis Figure A.3 presents four leading principal components of Atlantic SST variability. They exhibit a wide range of timescale ranging from interannual to multidecadal. None of the modes shows a secular trend as the Pacific influences have all been removed. These four modes should be uncorrelated with the Pacific PCs and uncorrelated with each other (Pacific PC should also be orthogonal to each other), and the near-zero cross correlation coefficients between each Pacific and Atlantic PC during 1901-2008 are manifested in Table A.1, to guarantee each PC is orthogonal to each other. Table A.1: correlations between Pacific and Atlantic PCs. ensop ensom trnd pdvpp ensonc pdvnp bien amo ninop ninom nao ensop 1 -6.27E-05 -4.01E-05 3.03E-05 -3.01E-05 2.28E-05 4.51E-06 4.25E-02 -4.66E-02 1.09E-02 4.20E-02 ensom -6.27E-05 1 -8.86E-05 6.69E-05 -6.64E-05 5.03E-05 9.95E-06 -3.89E-02 -2.30E-02 -4.85E-02 7.07E-03 trnd -4.01E-05 -8.86E-05 1 4.27E-05 -4.24E-05 3.22E-05 6.34E-06 -3.65E-02 -1.61E-02 -3.28E-02 2.55E-03 pdvpp 3.03E-05 6.69E-05 4.27E-05 1 3.20E-05 -2.42E-05 -4.80E-06 2.92E-02 8.06E-02 4.43E-03 1.95E-02 ensonc -3.01E-05 -6.64E-05 -4.24E-05 3.20E-05 1 2.42E-05 4.79E-06 4.38E-02 -1.39E-02 4.93E-02 9.94E-03 pdvnp 2.28E-05 5.03E-05 3.22E-05 -2.42E-05 2.42E-05 1 -3.59E-06 2.17E-02 -6.32E-03 9.31E-03 -4.68E-02 bien 4.51E-06 9.95E-06 6.34E-06 -4.80E-06 4.79E-06 -3.59E-06 1 -8.18E-03 5.68E-02 1.45E-02 9.22E-03 amo 4.25E-02 -3.89E-02 -3.65E-02 2.92E-02 4.38E-02 2.17E-02 -8.18E-03 1 -9.39E-05 -2.61E-06 5.17E-05 ninop -4.66E-02 -2.30E-02 -1.61E-02 8.06E-02 -1.39E-02 -6.32E-03 5.68E-02 -9.39E-05 1 -1.98E-06 4.07E-05 ninom 1.09E-02 -4.85E-02 -3.28E-02 4.43E-03 4.93E-02 9.31E-03 1.45E-02 -2.61E-06 -1.98E-06 1 1.18E-06 nao 4.20E-02 7.07E-03 2.55E-03 1.95E-02 9.94E-03 -4.68E-02 9.22E-03 5.17E-05 -1.98E-06 1.18E-06 1 114 Figure A.3: Principal components (PCs) of the four leading modes of recurrent spatiotemporal SST variability in the Atlantic (shaded) during 1901- 2008. Tick marks on the vertical axis denote three standard deviations. Smoothed PCs after 50 applications of a 1-2-1 smoother are represented in solid black lines. Atlantic Multidecadal Oscillation (AMO) is identified as the leading mode of residual Atlantic SST PCs, which accounts for 9.0% of the explained variance. Based on the presence of a complete cycle of AMO PC time series, it approximately has a 70-year period (Figure A.3). Recent warming is suggested not to be unique as same warming phenomenon also occurred during 1920-30s. The decaying and growing phase of interannual Atlantic SST variations in the eastern tropical basin are described in the following two leading modes, say AtlNino- and Atlnino+, which account for 7.3% and 6.8% of the explained variance, respectively. The fourth leading mode describes decadal variation of a north-south SST triple in the North Atlantic, 115 which is also referred as the Triple Mode, accounting for 4.2% of the explained variance. The mature phase of AMO mode suggests that positive center is located around 50°N (Figure A.4). Prior to the mature phase, SST anomalies begin to originate in the Davis Straits and the Labrador Sea, then propagate eastward in the midlatitudes, then equatorward in the eastern basin (along the Canary Current track), and finally westward in the subtropics, which is much like the North Atlantic subtropical gyre. We cannot neglect the importance of subtropical anomalies though they are relatively weak, and the eastern basin, or the whole equatorial Atlantic, does not show noticeable signals during the evolution process. However, this kind of AMO spatial structure differs from what mentioned in previous studies (Enfield et al., 2001), which is featured by the uniform polarity of SST anomalies across the entire North Atlantic, including the tropical basin. The suggested reason explaining the difference is because the previous study did not exclude the Pacific influence, but our present analysis does, so it reflects the intrinsic component of multidecadal oscillation variability. AMO structure with large amplitude in high latitude were also identified in some previous studies (e.g., Kawamura, 1994; Mestas-Nunez and Enfield, 1999; Parker et al., 2007). The key difference of our current analysis still lies in its independence from the Pacific and Indian Oceans. 116 Figure A.4: Spatiotemporal evolution of recurrent SST variability in the Atlantic: Seasonal lead-lag regressions of the leading principal components (PC 1-4). 117 The spatial modal evolution of Atlantic Nino suggests asymmetrical anomaly distributed about the equator but tilted northwestward from the African coast to the central equatorial Atlantic. During the growing phase (AtlNino-), SST anomalies are characterized by the development from east to west; and the decaying phase (AtlNino+) shows an opposite retreat. The mature phase of Atlantic Nino structure looks much similar as the Atlantic Nino mode in Ruiz-Barradas et al. (2000), of which the asymmetry about the equator becomes quite important. The reason why the amplitudes are reduced is because the Pacific influence has been ruled out. The modal spatial evolution suggested that the main lobe is located near the Gulf Stream separation region, then propagates eastward along 40°N, while the oppositely signed side lobes decays. The mature-phase spatial structure of this tripole mode bears is similar to NAO’s impact on winter SSTs (e.g. Visbeck et al., 2001; Nigam, 2003). A.3 Datasets As in the Pacific and Atlantic PC analysis, the U.K. Met. Office’s Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST 1.1; Rayner et al., 2003) is analyzed in both the Saji et al. analysis period (1958-1998) and the century-long period (1900-2007) in Section A.4. The extended-EOF analysis of filtered Indian Ocean SSTs in the period 1900-2007 is discussed in Section A.7. The linear regressions reported in Sections A.5 and A.7 are computed for the 1902-2005 period, given the end-point truncation of the principal components in extended-EOF analysis. The Simple Ocean Data Assimilation (SODA) 2.1.6 dataset is used for computing the upper-ocean heat content (Carton et al., 2000; Carton and Giese, 2008). The SODA 118 heat content (0-381m) analysis is confined to 1958-1998, the common period of the Saji et al. analysis (1958-1998) and the ECMWF 40-year Reanalysis (1958-2001; ERA40; Uppala et al., 2005). The SODA 2.1.6 ocean reanalysis is, interestingly, driven by the ERA40 winds, which are also analyzed here. The Climate Research Unit TS 3.1 (CRU TS3.1) dataset is used in the continental precipitation analysis reported in Section A.7. All datasets are analyzed at seasonal resolution. 119 Figure A.5: Evaluation of the basis for Indian Ocean Dipole (IOD) mode of SST variability from both recent (1958-1998) and century-long (1900-2007) SST records; the former period is identical to that analyzed by Saji et al. (1999): (a) Boreal fall SST regressions of the western-IOD box (W-IOD, marked) based SST index; (b) likewise for the eastern-IOD box (E-IOD, marked) index; note, the dipole mode index (DMI, Saji et al.) is obtained by subtracting the E-IOD index from the W-IOD one. Fall regressions of the W-IOD and E-IOD indices on ENSO-filtered SSTs are shown in (c)-(d); ENSO-filtered SSTs are generated, as noted in the text, by removing the spatiotemporally-evolving ENSO influence. Corresponding regressions for the century-long period are shown in the right column. In all cases, regressions are computed on detrended SST (following footnote 1). Contour interval and shading threshold is 0.05K, and positive values are shown in warm colors; the zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black. A.4 Variability in ENSO-filtered Indian Ocean SSTs: No Evidence for a Zonal-dipole Figure A.6: Two additional constructs of ENSO variability are filtered from Indian Ocean SSTs prior to 716 computation of the regressions of the W-IOD and E-IOD based SST indices: ENSO variability is constructed 717 from the canonical and non-canonical ENSO components only (i.e., without contribution from Pacific biennial 718 variability) in the left panels, and from the widely used index of ENSO’s mature-phase variability (Nino3.4 SST 719 index) in the right panels. Boreal fall regressions are for the century-long period (1900-2007); rest as in Figure A.5. Regressions significant at the 0.05 level are contoured in black. 120 The claim for a dipole structure of interannual variability in tropical Indian Ocean SSTs rests on the zonal structure of the second-leading EOF in Saji et al. (1999) – a dipole structure – that led to the DMI index: The index is from the difference of area-averaged SST anomalies in the western (50E-70E, 10S-10N) and eastern (90E-108E, 10S-EQ) tropical Indian Ocean; the regions are marked in black in Figure A.5. Although these regions are identified from a dipole-type EOF structure, there is no assurance that they are anti-correlated as EOFs often pick up structure to maximize explained variance (e.g., Dommenget and Latif, 2002). The connectedness of the western and eastern Indian Ocean (IO) regions – the two poles of the claimed dipole – is investigated in boreal fall in Figure A.5. The Saji et al. period (1958-1998) is analyzed first using the detrended SST anomaly record.2 Regressions of the western IOD box (50E-70E, 10S-10N; the W-IOD box) (Figure A.6a) support the dipole nature of variability but with a slightly southeastward displaced eastern center vis-à-vis IOD’s eastern box. Regressions of the eastern box (90E-108E, 10S-EQ; the E-IOD box), shown in Figure A.5b, are however focused in the south-central equatorial IO rather than the western box, diminishing reciprocity between the marked dipole centers. Support for the western box as a center of action of the dipole is even less in the century-long analysis (Figure A.5f) which exhibits vanishing regressions over this region! 2 The century-long (1900-2007) anomaly record was detrended by subtracting the projections of the SST Secular Trend mode (GN2008, Figure 13 and related discussion) from the seasonal anomalies. A linearly detrended SST anomaly record could just as well be used. A sub-period (1958-1998) of the detrended record is analyzed in the left panels while the full record (1900-2007) is analyzed in the right ones of Figure A.6. 121 The connectedness of the IOD centers crumbles upon removal of the ENSO signal from the IO SSTs, irrespective of the analysis period. The spatiotemporally varying ENSO influence was estimated by multiplying the time-dependent Pacific SST PCs of canonical ENSO variability (captured as two modes, the growth (ENSO) and decay (ENSO) modes in GN2008; Figure 3), non-canonical ENSO variability (ENSONC in GN2008; Figure 5), and biennial variability (GN2008, Figure 10) with their respective regressions on contemporaneous IO SSTs in the full record (1900- 2007).3 The Indian Ocean SSTs, filtered for this ENSO influence, are referred as filtered SSTs. The analysis of filtered SSTs in the Saji et al. period (1958-1998) shows the western dipole center to be non-viable: Regressions of the western IOD box have no footprint whatsoever over the eastern box (Figure A.5c), and likewise for regressions of the eastern IOD box (Figure A.5d) over the other region. The eastern box is again found more connected to the south-central equatorial IO SSTs,4 albeit more weakly than before. The century-long (1900-2007) analysis of filtered SSTs (Figures A.5g-h) corroborates the shorter-period findings, as does an inspection of the related correlation maps (not shown). Given the premium placed on characterization of ENSO’s influence on IO SSTs, it is, perhaps, essential to evaluate the sensitivity of our findings to different 3 ENSO is taken to consist of canonical, non-canonical, and biennial variability, consistent with prevailing views: For example, Rasmusson et al. (1990) viewed ENSO as superposition of biennial and lower-frequency variability; see Figure 8 in GN2008 for the impact of biennial variability on ENSO duration. That biennial variability is an integral component of ENSO is also indicated by correlations of the observed Nino3.4 SST index and its synthetic versions based on various SST-reconstructions. For example, if only canonical ENSO modes are used in SST reconstruction, the correlation is 0.84; correlation increases to 0.92 when the non-canonical ENSO mode is additionally included in the reconstruction, and even further to 0.95 when the biennial mode is included as well. The SST principal components are available online at http://dsrs.atmos.umd.edu/DATA/NIGAM/Diab.Heating/SST-PCs/ 4 The south-central equatorial IO region cannot be a viable western dipole-center given the vanishing regressions of the eastern center in this region in the longer period analysis (Figure A.6h). 122 estimations of this influence: Two estimates are considered in Figure A.6 where the fall-season connectivity of IOD centers is investigated in the century-long period (1900-2007). The first is based on canonical and non-canonical ENSO components alone, i.e., without additional contribution from Pacific biennial variability (as in Figure A.5); this is referred as ENSO-filtered (partial). A comparison of Figures A.5g and A.6a, and 1h and A.6b shows the continued non-viability of either of the IOD centers. The second estimate is based on the widely used index of ENSO’s mature- phase variability, the Nino3.4 SST index. Although ENSO’s influence is commonly estimated in this manner, the estimation is sub-optimal for reasons stated earlier (and in Compo and Sardeshmukh, 2010). Regardless, its filtering doesn’t lead to the emergence of any support for opposite-signed connectivity of the IOD centers (Figures A.6c-d). We find that ENSO-filtered Indian Ocean SSTs do not exhibit a zonal-dipole variability structure with centers of action as defined in Saji et al. (1999). The analysis of surface temperature does not support the DMI concept, which was shown to result from the inadvertent inclusion of ENSO’s spatiotemporal influence on Indian Ocean SSTs in the Saji et al. analysis. The search for a zonal-dipole variability structure in the tropical Indian Ocean is extended into the subsurface realm in section A.6, from analysis of the upper-ocean heat content variations. The search is rooted in regressions of the eastern pole of the IOD (E-IOD), which is a center of robust dynamical variability on account of regional ocean-atmosphere interaction including upwelling off Java-Sumatra. Meyers et al. 123 (2007) and Luo et al. (2010) have used the E-IOD index to characterize the internal mode of IO variability. A.5 The 1994 and 1997 Fall SST Anomalies: IOD Events? The Dipole Mode Index is strongly positive in the fall of 1994 and 1997 (Saji et al.1999, Figure 1); 1994 was a weak El Nino year but the 1997 winter saw one of the strongest warm events on record (NOAA Climate Diagnostic Bulletin). As ENSO’s influence is aliased in the DMI definition (cf. Figures A.5-A.6 and related discussion), it is instructive to examine the structure of both raw and ENSO-filtered Indian Ocean SST and upper-ocean (0-381m) heat content anomalies in these two years; the detrended SSTs and ocean heat content anomalies are the starting point in both cases. SST anomalies are based on the century-long (1900-2007) seasonal climatology, and reconstructed from PCs obtained from an extended-EOF analysis of SSTs in the same period. The ocean heat content anomalies are, of necessity, based on the shorter 41-year period (1958-1998) – the SODA data period, which is also the Saji et al.’s analysis period. The SST anomalies are displayed in Figures A.7a-j, while the heat content ones are shown in Figures A.7k-n. The fall 1994 SST anomaly (top panel) is strongly negative off Sumatra, i.e., over the E-IOD box, but weaker and of mixed sign over the W-IOD box. The DMI is, interestingly, positive (0.60 = (−0.03) − (−0.63)) with similar-signed anomalies in the two IOD boxes! The next 3 panels display various estimates of ENSO’s contribution to the fall 1994 SST anomaly: The one obtained from accounting of the spatiotemporal influence of canonical and non-canonical ENSO variability, and 124 biennial variability (i.e., the full ENSO influence) on IO SSTs is shown in Figure A.7b, while the version without the biennial component is in panel c. The IO SST anomaly constructed from regressions of the Nino3.4 SST index is shown in Figure A.7d. All 3 estimates of ENSO’s influence on IO SSTs exhibit a zonal-dipole structure, one that would, undoubtedly, project on the DMI index. The ENSO- unrelated IO SST anomaly (i.e., observed minus ENSO’s full influence) is shown in the bottom panel. The SST anomalies off Sumatra are not as strong as before (Figure A.7a) while the ones along the Somali Coast and Arabian Sea are stronger (and of the same sign as the Sumatra ones) in the ENSO-filtered version (Figure A.7e); yielding a smaller DMI (0.25) than the one calculated from unfiltered anomalies (Figure A.7a; 0.60). The fall 1997 SST anomalies are analyzed in the right column of Figure A.7. All three estimates of the ENSO contribution to IO SSTs are now stronger; not surprising, given the unusually strong 1997 El Nino. The large-scale anomaly structure (zonal-dipole) is very similar (but for the amplitude) to that seen in the 1994 contributions (Figures A.7b-d), and also among themselves (Figures A.7g-i). The ENSO-filtered IO SST anomaly in 1997 (Figure A.7j), interestingly, does not exhibit an impressive zonal-dipole structure in the IO. DMI index from the original (Figure A.7f) and ENSO-filtered SST anomalies (Figure A.7j) is 1.05 and 0.09, respectively. A full accounting of ENSO’s contribution to the Indian Ocean surface (SST) anomalies in fall 1994 and 1997 suggests a weak IOD event in 1994, and none at all in 1997. The accounting leads to a rather different assessment in 1997, which was 125 marked as a strongly +ve IOD year in previous analyses (e.g., Saji et al., 1999), based on DMI computation off raw SST anomalies. An IOD event was also reported in fall 2006 (Horii et al., 2008), with a DMI index (0.57) similar to that in fall 1994; both weak El Nino years. Reconstruction of the IO SST anomalies indicates a strong ENSO contribution whose filtering leads to a DMI index of 0.14, or a weak IOD event, much as in fall 1994. 126 Figure A.7a-j: Origin of the Fall 1994 and 1997 Indian Ocean SST anomalies: Observed anomalies are in the top panels (a,f); ones reconstructed from ENSO’s full influence (canonical, non-canonical, biennial modes) in the next two panels (b,g); those reconstructed from ENSO’s partial influence (canonical, non-canonical modes) in the middle panels (c,h); those reconstructed using ENSO’s mature-phase index (Nino3.4 SST index) in the next from the bottom panels (d,i), while the fully ENSO-filtered anomalies are in the bottom panels (e,j). Detrended SSTs are the starting point in all cases. The western and eastern IOD box is displayed in all panels to facilitate IOD recognition. Contour interval and shading threshold is 0.15K, and positive values are shown in warm colors; the zero-contour is omitted in all panels. 127 Subsurface Anomaly Structure In view of the above-noted disagreement in characterization of the nature/structure of the IO SST anomalies in fall of 1994, and especially 1997, the related subsurface anomaly structure is examined in Figures 7k-n though plots of the upper-ocean heat content; as before, both observed and ENSO-filtered versions are shown. Interestingly, the ENSO-filtered anomaly is nearly indistinguishable from the original one in 1994, but very distinct in 1997 (a strong El Nino year). A zonal-dipole structure is evident in 1994 but with the western pole located to the right of the W- IOD box, but no dipole structure is manifest in the filtered 1997 heat-content anomaly (Figure 7n). Strong southeasterlies are present on the west coast of Sumatra in both cases (Figures 7k,m), leading to coastal upwelling and colder SSTs. The ENSO- filtered winds (Figures 7l,n) remain southeasterly but the filtered SSTs off Sumatra are cold only in fall 1994, indicating the significance of other influences on SST in fall 1997. The similarity of the observed and ENSO-filtered winds attests to the rather modest influence of ENSO on Indian Ocean winds (~1-2 m/s), consistent with earlier diagnoses of this signal (e.g., Nigam and Shen, 1993, see their Figure 2). The fall 2006 anomalies (not shown) are similar to the 1994 ones in that their ENSO-filtered versions exhibit weak IOD variability in SST but robust dipole variability in subsurface temperatures. This limited analysis (two cases) suggests that the tropical IO basin can exhibit internally-generated zonal-dipole structure, but at the subsurface. Externally-driven IO variability (e.g., from ENSO’s influence), on the other hand, can generate zonal- dipole type structure both at the surface and subsurface. If this analysis is 128 corroborated from additional observational and modeling studies, it would caution against identification of IOD events from surface analysis. Figure A.7k-n: Origin of the Fall 1994 and 1997 Indian Ocean upper-ocean (0-381m) heat content anomalies: Those from SODA 2.1.6 reanalysis (Carton et al. 2000) are in the top panels (k,m), while the fully ENSO-filtered ocean heat content anomalies are in the bottom panels (l,n). The anomalies are displayed after 4 applications of the 9-point smoother (smth9 in GrADS) with a shading threshold/interval of 2x108 J/m2; see color bar. Detrended ocean heat-content is the starting point in all cases. The western and eastern IOD box is displayed in all panels to facilitate IOD recognition. A.6 Seasonally Evolving Subsurface Variability in the Indian Ocean: ENSO’s Influence The upper-ocean (0-381m) heat content is an integrated measure of ocean temperature, reflecting input from surface fluxes, advective transports, and upwelling; it is thus a more steady measure of the upper-ocean state than SST (Merle, 1980). 129 Positive/negative ocean heat content (OHC) anomalies are generally indicative of positive/negative subsurface temperature anomalies, a deepening/shoaling thermocline, and at times, anomalously high/low SSTs. The seasonal evolution of OHC and 850 hPa wind anomalies associated with fall SST variability at the E-IOD box – the more viable dipole center – is displayed in Figure A.8; the 850 hPa winds, and not surface winds, are displayed to preclude use of extrapolated winds over adjoining landmasses. Regressions of 3 versions of the E- IOD SST index are shown: Lead/lag regressions of the ‘raw’ (i.e., unfiltered) index are in the top row; of an index constructed from just ENSO related IO variability in the middle; while those of the ENSO-filtered E-IOD index are in the bottom row. Seasonal anomalies are analyzed: The summer (JJA) anomalies (first column) are the 1-season lead regressions while the winter (DJF) ones (last column) are the 1-season lag regressions; contemporaneous (fall, SON) ones are in the middle column. Regressions of the E-IOD will naturally bring out the –ve phase IOD structure, when temperatures and OHC off Sumatra are above normal, and strong westerly anomalies blow across the western coast of Sumatra (Figure A.8b); the alongshore winds being from the northwest lead to coastal downwelling. Equatorial westerlies, on the other hand, generate downwelling over the open ocean, leading to thermocline depression over the eastern equatorial basin. Downwelling equatorial and coastally trapped Kelvin waves also contribute to the positive OHC anomalies in this region (Wyritki, 1973; Clarke and Liu, 1993). The southwesterly anomalies across the southwestern equatorial IO basin generate upwelling Rossby waves, leading to shallower thermocline and negative OHC anomalies here (Figure A.8b; Meehl et al. 130 2003); related regional feedbacks have also been proposed (e.g., Saji et al., 1999). The off-shore westerly anomalies (and related moisture transports), evident especially off Somalia in Figure A.8b, have, of course, been implicated in regional drought over eastern Africa (Saji et al., 1999; Goddard and Graham, 1999; Latif et al., 1999). The top row shows the E-IOD linked OHC anomalies in the eastern IO to peak in fall but anomalies in the southwestern equatorial IO basin attain maximum amplitude in boreal winter. The precursor (JJA) near-surface wind anomalies include easterlies over peninsular India, which weaken the climatological monsoonal flow (southwesterly), leading to diminished summer monsoon rainfall. Regardless of the season, OHC anomalies in the western basin are seldom focused in the W-IOD box (marked). The strong westerly anomalies over the tropical IO basin and stronger easterlies over the north-equatorial Pacific are part of a larger east-west atmospheric circulation anomaly centered over the Maritime Continent, as apparent from more expansive plots (not shown). The larger-scale anomaly structure is indicative of ENSO’s influence, which is analyzed by computing both the E-IOD index and its regressions from just the ENSO signal over the IO basin (Figures A.8d-f); extraction of the full ENSO signal was discussed earlier. That the ENSO signal dominates the regressions of the ‘raw’ E-IOD index is evident from the striking similarity of the top two rows. ENSO’s impressive influence on the IO is, of course, transmitted through its impact on surface winds over the IO basin (e.g., Nigam and Shen, 1993; Xie et al., 2002). Regardless, the case for filtering the externally-generated ENSO influence from IO 131 variability prior to an objective search for recurrent variability structures across this basin is reinforced by the above analysis. Finally, the structure of internally-generated variability in the IO basin is revealed in Figures A.8g-i (bottom row) via regressions of the E-IOD index; the index and regressions are both obtained from ENSO-filtered IO variability. Not unexpectedly, OHC anomalies are now present primarily in the IO basin. But somewhat surprisingly, a coherent dipole structure is manifest, with the western basin anomalies somewhat better positioned vis-à-vis the W-IOD box. The analysis supports the finding of the previous section, viz., that internally-generated variability in the tropical IO basin can exhibit a zonal-dipole structure at the subsurface. Figure A.8: Fall-centered seasonal lead-lag regressions of 3 versions of the E-IOD SST index on upper-ocean (0-381m) heat content (OHC, from SODA 2.1.6 ocean reanalysis) and 850 hPa winds (from ERA-40 atmospheric 132 reanalysis): Regressions of the ‘raw’ (i.e., unfiltered) index are in the top row; of an index constructed from just ENSO related IO variability in the middle; while those of the ENSO-filtered E-IOD index are in the bottom row. Note, the index and regressions are obtained from the same data set in all cases. The summer (JJA) anomalies (first column) are the 1-season lead regressions while the winter (DJF) ones (last column) are the 1-season lag regressions; contemporaneous (fall, SON) ones are in the middle column. The OHC regressions are shaded/contoured in warm (positive) and cool (negative) colors with an interval/ threshold of 108 J/m2: see color bar. Wind vector scales are shown below each panel. All data sets were detrended (following footnote 1) and the zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black. A.7 Recurrent Modes of Variability in Filtered Indian Ocean SSTs An objective search for the recurrent modes of variability of Indian Ocean SSTs is undertaken – in contrast to the analysis of variability patterns in the preceding sections. The distinction between mode and pattern is important as the former refers to a unique spatiotemporal variability structure while the latter to just a spatial pattern of variability that is linked, potentially, with more than one timescale. Another difference is that previous analyses factored for just ENSO’s influence on the IO basin, while the present one will also factor for the influence of Pacific decadal SST variability. Finally, the previous sections were too focused on the E-IOD because it was the more viable of the IOD poles, but it is unclear if the E-IOD region would emerge as a key variability center in an objective analysis of filtered IO SSTs. A.7.1 SST Filtering Recurrent spatiotemporal structure of SST variability in the Indian Ocean is analyzed after filtering the influence of both interannual and decadal Pacific SST variability. The influence of ENSO variability was estimated in section 3; the impact of decadal variability is, likewise, estimated from regressions of the Pan Pacific and 133 North Pacific PCs (GN2008). The first mode, with horse-shoe structure in the Pacific, exhibits connections to the tropical-subtropical Atlantic resembling the AMO. The second, capturing the 1976/77 climate-shift, is similar to Pacific Decadal Oscillation (Mantua et al., 1997) in structure but with interesting links to the IO SSTs (cf. Figure 12 in GN2008). The SST record was already detrended earlier using regressions of the nonstationary SST Secular Trend (cf. footnote 1); this mode captures the wide- spread but non-uniform warming of all basins along with a sliver of cooling in the central equatorial Pacific.5 A.7.2 Analysis Technique Table A.2: Sensitivity Analysis: T0 is the primary analysis, which was perturbed as listed below to assess robustness of the extracted variability modes. In all cases, SST was detrended using the nonstationary SST Secular Trend mode (see footnote 1), and is referred as ‘DeTrend’ here. ‘0’ in the ‘Rotated’ column indicates that loading vectors were not rotated. Name Domain Period Rotated Sampling Window SST T0 Indian Ocean 1902- 2007 2 7- season DeTrend – Pacific basin influence T1 Indian Ocean 1902- 2007 0 7- season DeTrend – Pacific basin influence T2 Indian Ocean 1902- 2007 3 7- season DeTrend – Pacific basin influence T3 Indian Ocean 1902- 2007 0 5- season DeTrend – Pacific basin influence T4 Indian Ocean 1902- 2007 0 7- season DeTrend – Pacific and Atlantic basin influences T5 Indian Ocean 1902- 2007 0 7- season DeTrend 5 The physicality of the decadal modes was evaluated using analog counts and fish recruitment records in GN2008. 134 The filtered, seasonal SST anomalies during 1902-2007 were analyzed in the Indian Ocean basin (40S-30N, 20E-120E) using the extended-EOF technique (Weare and Nasstrom, 1982); 7-season long anomaly sequences were targeted in the primary analysis (T0) where two leading loading vectors are rotated. Robustness of the variability modes was ascertained by perturbing the primary analysis: No rotation of loading vectors (T1); rotation of 3 loading vectors (T2); shorter (5-season long) sampling window (T3); Indian Ocean SSTs additionally filtered for Atlantic’s influence (T4); Indian Ocean SSTs are only detrended but not filtered for any external influences (T5); all listed in Table A.2. Table A.3 describes the sensitivity analysis results, including temporal correlation of the PCs of the primary and perturbed analyses. An observational realization of the mode is the ultimate proof of its physicality. But it is seldom that observed anomalies are composed of just one mode of variability (i.e., with all other modes suppressed at that time); of course, should this happen, an observational “analog” of that mode is encountered. The number of observational analogs of an extracted set of modes in the anomaly record is one objective measure of the “physicality” of that extraction, and the primary analysis is chosen in this manner. An observed anomaly will be deemed a modal analog if any one PC is larger than all others in that analysis by at least one unit of magnitude; note, PCs are orthonormal with or without rotation. The identification is objective and easily implemented, and Table A.4 lists the number of analogs in the six analyses (T0-T5): Considering just the first two PCs, 128 seasonal anomalies (out of the 424 analyzed) are found to be analogs of the first or second mode in the T0 analysis. 135 Table A.3: Sensitivity Results: Column 1 shows the leading modes identified in the primary analysis (T0); numbers following the name indicate the percentage variance explained by the mode and its rank, respectively. Columns 2–6 list attributes of the leading modes in the five sensitivity tests (T1–T5), with the three-slash numbers indicating correlation between the test case and primary analysis PC, the percentage variance explained by that mode, and its rank in the test analysis. Asterisks indicate ‘no match’. All correlation coefficients are significant at the 0.001 level. T0 T1 T2 T3 T4 T5 IND1/10.9/1 0.96/11.4/1 0.98/10.3/1 0.94/12.6/1 0.92/10.9/1 0.67/6.7/3 IND2/5.5/2 0.96/5.0/2 0.82/5.1/3 0.88/6.3/2 0.71/5.2/2 * Table A.4: Analog Counts: Number of “analogs” in six extended-EOF analyses of Indian Ocean SST variability during 1902-2007. An analog is identified when the absolute value of the PC of any one mode is larger than that of all others by at least one unit; note PCs are orthonormal. Mathematically, if | |PCi| - |PCj| | > 1.00 for all j not equal than i, an analog is counted. For consistent evaluation, only analogs of the two leading modes are counted. T0 T1 T2 T3 T4 T5 128 96 87 98 97 95 A.7.3 Analysis Results The two leading PCs from the primary analysis (T0) are shown in Figure A.9a. The leading PC accounts for ~11% of the variance, and represents sub-decadal variability with a notable warm-phase in the 1970s-1980s. As seen later (Figure A.10), the represented SST variability has a meridional dipole structure with a pronounced southern pole in the subtropical Southern Hemisphere (SH) basin (focused southward of Mauritius); this region was anomalously warm during the 1970s-80s. The second PC (explaining 5-6% of variance) represents, principally, interannual variability; the related loading vector (Figure A.10) is more tropically 136 focused and resembles the IOD structure in the western-central basin; the similarity with IOD however ends here, as shown and discussed later. Unlike the leading mode, this one exhibits considerable structural evolution – evolving into a meridional dipole (with a stronger northern pole) in the SH tropical-subtropical basin before dissipating (Figure A.10). Potential links with the IOD and the Subtropical IOD (SIOD; Behera and Yamagata, 2001) are investigated in Figures A.9b and A.10. The autocorrelation structure of the PCs is examined in Figure A.9b to estimate the modal time scale; a conservative estimate follows from the temporal distance of points where autocorrelation is e1 (≈ 0.37); yielding ~6 years for PC-1 and 2-3 years for PC-2. The autocorrelation structure of the IOD and SIOD indices is also displayed in Figure A.9b, for context; both, evidently, represent much shorter timescale (~1 year) variability. The correlation of PC-1 and the IOD (SIOD) index is 0.14 (0.47), significant at the 0.001 level. The corresponding PC-2 correlations are −0.11 and 0.22, at the 0.05 significance level. The IOD is thus unrelated to the extracted modes which exhibit some pattern (but not timescale) similarity to SIOD variability. Spatiotemporal evolution of the SST variability modes in the T0 analysis is displayed in the first two columns of Figure A.10. Seasonal lead-lag regressions of the fall-season (SON) PCs are shown in order to compare with the IOD pattern (3rd column) which is robust in fall; the fourth column depicts evolution of the Subtropical IOD using its winter (DJF) based index. The 5-season evolution (time running downward) clearly conveys the variability time scales: significantly longer than 5 seasons for PC-1, just about this period for PC-2, and shorter than 5-seasons for both IOD and SIOD; consistent with the autocorrelation analysis (Figure A.10b). 137 The evolution of the leading mode is not seasonally sensitive in view of its ~6 year time scale but others are, as manifest from the regressions centered on other seasons (not shown). As noted earlier, the leading mode is focused in the SH with large amplitudes in the western subtropics. The second mode, in contrast, has footprints in the northern basin as well, and shows significant spatiotemporal development: evolving from a single-signed anomaly pattern in the Tropics in spring- fall into a meridional dipole in the SH subtropics in boreal winter (DJF); note the similarity of the winter patterns in the SH6 . The second mode dissipates by the following spring (MAM) whereas the first mode shows no such sign. A comparison of column 2 and 3 indicates both similarities and differences between the second mode of this analysis and the IOD. Key differences include the lack of modal amplitude in the western basin (including Arabian Sea) in fall (SON), and absence of zonal-dipole structure in the Tropics (e.g., no modal amplitude off Sumatra). The similarity accrues from the amplitude focus in the central basin in boreal fall-winter (SON-DJF), especially structure in the SH Tropics (e.g., northwest- to-southeast amplitude tilt). The regression patterns are most similar in boreal winter (DJF) when both exhibit some resemblance with the Subtropical IOD pattern in the SH. The analysis indicates that the SH subtropical Indian Ocean exhibits SIOD-type variability patterns, but on both decadal and interannual timescales. The interannual one undergoes significant metamorphosis, beginning with single-signed anomalies in the Tropics and ending in a SIOD-type structure in boreal winter (DJF). 6 An extended-EOF analysis is well positioned to distinguish variability structures having similar spatial footprint but different time scales and evolution. 138 Figure A.9: Principal components (PCs) of the two leading modes of recurrent spatiotemporal SST variability in the Indian Ocean are displayed in the upper panel. The PCs are obtained from rotated, extended-EOF analysis of Pacific-uninfluenced, non-seasonal Indian Ocean SST variability during 1902- 2007; PC-1 and PC-2 explains 10.9% and 5.5% of the variance, respectively. Tick marks on the vertical axis are drawn every three standard deviations. The autocorrelation structure of the twos PCs and the IOD and Subtropical-IOD (SIOD) indices is shown in the lower panel. 139 Figure A.10: Spatiotemporal evolution of recurrent SST variability in the Indian Ocean: Seasonal lead-lag regressions of the leading principal components 140 (PC-1 and PC-2) are shown in the left two columns, with time increasing downward. Fall-centered regressions (i.e., of the fall PC values) are shown to facilitate comparison with the IOD structure which is robust in boreal fall (SON). Indian Ocean SSTs are filtered for the Pacific’s influence (as in T0 analysis) prior to computation of regressions. Fall-centered regressions of the IOD, and winter-centered (DJF) regressions of the Subtropical IOD on unfiltered Indian Ocean SSTs are shown in the far right columns. The western and eastern IOD box is displayed in all panels to facilitate IOD recognition. Regressions are for the 1902-2007 period, and detrended SSTs are the starting point in all cases. Contour interval and shading threshold is 0.05K, and positive values are shown in warm colors; zero-contour is omitted in all panels. Regressions significant at the 0.05 level are contoured in black. A.7.4 Link with Continental Precipitation Precipitation links of the IO SST variability modes are shown in Figure A.11, beginning with boreal summer (JJA, top panels); contemporaneous correlations are shown to gauge significance of the links. PC-1, representing sub-decadal SST variability in the subtropical SH basin, is linked with diminished summer rainfall over sub-Saharan Africa (including Sahel). The links are weak elsewhere except Southeast Asia where correlations are significantly positive. During boreal fall (SON), PC-1 continues to be linked with rainfall deficits over western-central Sahel, but positive links emerge over the Congo basin and to its southeast; rainfall variability in these regions has been linked to higher SST in the southwestern Indian Ocean (Behera and Yamagata, 2001; Xie and Arkin, 1996; Reason, 2001). In austral summer (DJF), PC-1 is positively linked with increased rainfall in northern–central Australia; the correlation distribution suggests a southward excursion of the monsoon rainfall belt. Precipitation links of the interannual SST mode (PC-2) are generally weaker (Figure A.11, right column). The links have some structure in fall (SON) when Sahel 141 is wetter and eastern equatorial Africa marginally drier; southeastern Australia is wetter during this time (austral spring). Figure A.11: Seasonal regressions/correlations of the Indian Ocean SST principal components on the Pacific-uninfluenced Indian Ocean SST and adjoining continental precipitation: JJA (top), SON (middle), and DJF (bottom). SST regressions are shaded and contoured in warm (positive) and cool (negative) colors with an interval/threshold of 0.05K. Correlations with CRU TS3.1’s half- degree continental precipitation are shown with a shading threshold of 0.10 and an interval of 0.05; see color bar. All data sets were detrended (following footnote 1) and the zero-contour is omitted in all panels. Regressions/correlations significant at the 0.05 level 798 are contoured in black. 142 A.8 Conclusion The viability of the Indian Ocean Dipole (IOD) mode of SST variability is investigated almost a dozen years after it was proposed (Saji et al., 1999), and just about a decade after intense discussion of its physicality (Allan et al., 2001; Dommenget and Latif, 2002; Hastenrath, 2002; Behera, et al. 2003; Yamagata et al., 2003; and others). At times, the discussion became lively, as with Philander’s (2003) metaphoric use of a cow for IOD’s complex structure. Given the close scrutiny, especially of IOD’s relationship with ENSO in these critiques and follow-on studies, there is some obligation to state the reasons for revisiting this issue. Our primary interest, in fact, was assessment of the impact of Indian Ocean (IO) SSTs on Chinese summer rainfall, as part of a broader effort to understand the reasons for the drying of northern China in recent decades. Towards this end, we investigated the structure of recurrent spatiotemporal variability of IO SSTs in the 20th century, factoring for Pacific SSTs’ influence (including ENSO). To our surprise, the IOD was not among the extracted modes in our surface (SST) analysis, piquing our interest in Saji et al.’s study. Interestingly, it is neither new observational data nor deployment of refined analysis techniques or improved climate models that lead to the finding that IOD is, in part, an artifact of ENSO’s influence on the IO basin. The zonal-dipole SST structure – the defining IOD attribute – crumbles when ENSO’s influence on IO SSTs is removed prior to any analysis of IO SST variability. Saji et al. (1999) took their leading, seasonally-invariant SST EOF, with monopolar structure in the IO basin, to 143 be ENSO’s entire influence,7 and proceeded to define the Dipole Mode Index (DMI) on the basis of the second EOF’s structure. As the EOFs are temporally orthogonal, Saji et al. viewed the IOD (based on DMI) to be ENSO-independent. ENSO’s full influence on IO SSTs was, of course, not captured by Saji et al.’s first mode. Although capturing ENSO’s full influence is non-trivial, as discussed in Guan and Nigam (2008) and Compo and Sardeshmukh (2010), and earlier in section 3 and 6a of this study, we show that a zonal-dipole variability structure is untenable even when a rudimentary estimate of ENSO’s influence on IO SSTs (from Nino3.4 SST index regressions) is filtered in advance of IO SST variability analysis. Not surprisingly, Dommenget (2011) reached a similar conclusion, viz., variations of the IO zonal SST gradient (or DMI) are not independent of ENSO, using the DEOF analysis technique (Dommenget, 2007) where a first-order auto-regressive model is used to generate the null-hypothesis for the EOF patterns. This analysis, like ours, revealed no coherent connection between the two DMI centers at the surface. Meyers et al. (2007) also sought to remove the ENSO signal from the IO SSTs prior to the IO analysis, using a lagged-EOF technique that seeks to capture a spatiotemporally evolving signal (e.g., ENSO) as a single mode. The efficacy of this technique in capturing ENSO’s complex variability was however not assessed by the authors; for example, from correlations of the observed and reconstructed Nino3.4 SST indices. The method involves restoring phase lags – a daunting task when the variability mode exhibits multiple timescales, not all of which are a priori known. Some differences in assessment of the IOD-ENSO relationship among analyses are 7 This too would not be different from the physicist representing the cow as a sphere in Philander’s brief perspective on the IOD-ENSO debate titled Of Dipoles and Spherical Cows (2003). 144 thus expected. Note, our extended-EOF analysis (Weare and Nasstrom, 1982) is equivalent to the widely used multichannel singular spectrum analysis (von Storch and Zwiers, 1999). Unable to find a zonal-dipole mode of internal variability at the surface (SST), we extended the search to the subsurface realm, specifically, upper-ocean (0-381m) heat content (OHC) variations. The search was first conducted in context of the two recent fall-period IO SST anomalies (1994 and 1997), both of which are classified as +ve IOD events (e.g., Saji et al., 1999), with DMI index of +0.60 and +1.05, respectively. Deconstruction of the fall SST anomalies revealed a DMI index of +0.25 and −0.09 in the ENSO-filtered SST anomalies, i.e., a weak IOD event in 1994 and none at all in 1997. Deconstruction of the OHC anomalies did indicate a zonal-dipole structure in 1994 but with the western pole positioned to the right of the western IOD box; no dipole was however manifest in 1997 in the ENSO-filtered OHC. This limited analysis suggests that internally-generated variability in the tropical IO basin can exhibit zonal-dipole structure but, perhaps, only at the subsurface. This finding, based on 2 cases, was corroborated from regressions of the E-IOD index (based on SST anomalies in the eastern IOD box) on OHC over a 40-year period. Finally, we report on an objective search for the recurrent modes of variability of IO SSTs from which the influence of both ENSO and Pacific decadal SST variability has been removed. The reported search for the internally-generated modes of IO variability is different from prior analyses which typically focus on recurrent spatial patterns, in contrast with spatiotemporal patterns in the present search. Interestingly, none of the extracted modes was found to exhibit a zonal-dipole 145 structure in SST. The second leading mode does have a dipole-like structure but in the meridional direction and in boreal fall/winter, when it resembles the Subtropical IOD (SIOD) pattern. The resemblance however stops there as this mode’s evolution timescale is considerably longer than the SIOD’s. The seasonally evolving continental precipitation links of the two SST variability modes contain interesting features over sub-Saharan Africa, the Congo basin, and northern-central Australia. Our study calls into question the edifice built on the DMI index, including constructs such as EQUINOO (IOD-related zonal surface winds at the equator, Gadgil et al., 2004). 146 Bibliography Ablowitz, M. J., J. Hammack, D. Henderson, and C. M. Schober, 2001: Long-time dynamics of the modulational instability of deep water waves. Physica D, 152, 416-433. Allan, R. J., D. Chambers, W. Drosdowsky, H. Hendon, M. Latif, N. Nicholls, I. Smith, R. Stone, C. Roger and Y. Tourre, 2001: Is there an Indian Ocean dipole and is it independent of the El Niño-Southern Oscillation?, CLIVAR Exchanges, 6 (3 (no.21)), 18-22. Alpert, J. C., and V. K. Kumar, 2007: Radial wind super-obs from the WSR-88D radars in the NCEP operational assimilation system. Mon Weather Rev, 135, 1090-1109. Anderson, J., T. Hoar, K. Raeder, H. Liu, N. Collins, R. Torn, and A. Avellano, 2009: The Data Assimilation Research Testbed a Community Facility. B Am Meteorol Soc, 90, 1283-1296. Annan, J. D., 2004: On the orthogonality of bred vectors. Mon Weather Rev, 132, 843-849. Ansmann, G., R. Karnatak, K. Lehnertz, and U. Feudel, 2013: Extreme events in excitable systems and mechanisms of their generation. Phys Rev E, 88. Arecchi, F. T., U. Bortolozzo, A. Montina, and S. Residori, 2011: Granularity and Inhomogeneity Are the Joint Generators of Optical Rogue Waves. Phys Rev Lett, 106. Ashok, K., Z. Y. Guan, and T. Yamagata, 2003: A look at the relationship between the ENSO and the Indian Ocean Dipole. J Meteorol Soc Jpn, 81, 41-56. Ashok, K., S. K. Behera, S. A. Rao, H. Y. Weng, and T. Yamagata, 2007: El Niño Modoki and its possible teleconnection. J Geophys Res-Oceans (1978–2012), 112. 147 Balci, N., A. L. Mazzucato, J. M. Restrepo, and G. R. Sell, 2012: Ensemble Dynamics and Bred Vectors. Mon Weather Rev, 140, 2308-2334. Banfield, D., B. J. Conrath, J. C. Pearl, M. D. Smith, and P. R. Christensen, 2000： Thermal tides and stationary waves on Mars as revealed by Mars Global Surveyor Thermal Emission Spectrometer, J. Geophys. Res., 105, 9521–9537. Banfield, D., B. J. Conrath, M. D. Smith, P. R. Christensen, and R. J. Wilson, 2003: Forced waves in the Martian atmosphere from MGS TES nadir data, Icarus, 161, 319–345. Baquero-Bernal, A., M. Latif, and S. Legutke, 2002: On dipolelike variability of sea surface temperature in the tropical Indian Ocean. J. Climate, 15, 1358-1368. Barlow, M., S. Nigam, and E. H. Berbery, 2001: ENSO, Pacific decadal variability, and US summertime precipitation, drought, and stream flow. J Climate, 14, 2105-2128. Baschek, B., and J. Imai, 2011: Rogue Wave Observations Off the US West Coast. Oceanography, 24, 158-165. Bascom, W., 1964: Waves and Beaches, Doubleday, p58-59. Behera, S. K., R. Krishnan, and T. Yamagata, 1999: Unusual ocean-atmosphere conditions in the tropical Indian Ocean during 1994. Geophys. Res. Lett., 26, 3001-3004. Behera, S. K., and T. Yamagata, 2001: Subtropical SST dipole events in the southern Indian ocean. Geophys. Res. Lett., 28, 327-330. Behera, S. K., S. A. Rao, H. N. Saji, and T. Yamagata, 2003: Comments on “A Cautionary Note on the Interpretation of EOFs”. J. Climate, 16, 1087-1093. Behera, S. K., J. J. Luo, S. Masson, S. A. Rao, H. Sakuma, and T. Yamagata, 2006: A CGCM study on the interaction between IOD and ENSO. J. Climate, 19, 1688-1705. Boussinesq, J., 1877: Essai sur la théorie des eaux courantes, Imprimerie nationale. 148 Bowler, N. E., 2006: Comparison of error breeding, singular vectors, random perturbations and ensemble Kalman filter perturbation strategies on a simple model. Tellus A, 58, 538-548. Cai, M., E. Kalnay, and Z. Toth, 2003: Bred vectors of the Zebiak-Cane model and their potential application to ENSO predictions. J Climate, 16, 40-56. Calini, A., and C. M. Schober, 2002: Homoclinic chaos increases the likelihood of rogue wave formation. Phys Lett A, 298, 335-349. Calini, A., N. M. Ercolani, D. W. McLaughlin, and C. M. Schober, 1996: Mel'nikov analysis of numerically induced chaos in the nonlinear Schrodinger equation. Physica D, 89, 227-260. Cane, M. A., and Coauthors, 1997: Twentieth-century sea surface temperature trends. Science, 275, 957-960. Chiang, J. C. H., and D. J. Vimont, 2004: Analogous Pacific and Atlantic meridional modes of tropical atmosphere-ocean variability. J Climate, 17, 4143-4158. Chung, C., and S. Nigam, 1999: Weighting of geophysical data in Principal Component Analysis. J Geophys Res-Atmos, 104, 16925-16928. Carrassi, A., A. Trevisan, and F. Uboldi, 2007: Adaptive observations and assimilation in the unstable subspace by breeding on the data-assimilation system. Tellus A, 59, 101-113. Carter, J. D., and C. C. Contreras, 2008: Stability of plane-wave solutions of a dissipative generalization of the nonlinear Schrodinger equation. Physica D, 237, 3292-3296. Carton, J. A., G. Chepurin, and X. Cao, 2000: A Simple Ocean Data Assimilation analysis of the global upper ocean 1950-1995 Part 2: results. J. Phys. Oceanogr., 30, 311-326. Carton, J. A., and B.S. Giese, 2008: A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA). Mon Weather Rev, 136, 2999-3017. 149 Chabchoub, A., N. P. Hoffmann, and N. Akhmediev, 2011: Rogue Wave Observation in a Water Wave Tank. Phys Rev Lett, 106. Chambers, D., B. Tapley, and R. Stewart, 1999: Anomalous warming in the Indian Ocean coincident with El Nino. J Geophys Res-Oceans (1978–2012), 104, 3035-3047. Chapman, S. and R. S. Lindzen, 1970: Atmospheric Tides, D. Reidel Press, Dordrecht, Holland, 200 pp. Cho, A., 2011: Ship in Bottle, Meet Rogue Wave in Tub. Science, 332, 774-774. Christensen, P. R., and Coauthors, 1992: Thermal Emission Spectrometer Experiment - Mars-Observer Mission. J Geophys Res-Planet, 97, 7719-7734. Christensen, P. R., and Coauthors, 2001: Mars Global Surveyor Thermal Emission Spectrometer experiment: Investigation description and surface science results. J Geophys Res-Planet, 106, 23823-23871. Clarke, A. J., and X. Liu, 1993: Observations and dynamics of semiannual and annual sea levels near the eastern equatorial Indian Ocean boundary. J. Phys. Oceanogr., 23, 386-399. Compo, G. P., and P. D. Sardeshmukh, 2010: Removing ENSO-Related Variations from the Climate Record. J. Climate, 23, 1957-1978. Conrath, B. J., J. C. Pearl, M. D. Smith, W. C. Maguire, P. R. Christensen, S. Dason, and M. S. Kaelberer, 2000: Mars Global Surveyor Thermal Emission Spectrometer (TES) observations: Atmospheric temperatures during aerobraking and science phasing. J Geophys Res-Planet, 105, 9509-9519. Corazza, M., and Coauthors, 2003: Use of the breeding technique to estimate the structure of the analysis "errors of the day". Nonlinear Proc Geoph, 10, 233- 243. Coy, L., S. D. Eckermann, K. W. Hoppel, and F. Sassi, 2011: Mesospheric Precursors to the Major Stratospheric Sudden Warming of 2009: Validation and 150 Dynamical Attribution Using a Ground-to-Edge-of-Space Data Assimilation System. Journal of Advances in Modeling Earth Systems, 3, M10002. Danforth, C. M., E. Kalnay, and T. Miyoshi, 2007: Estimating and correcting global weather model error. Mon Weather Rev, 135, 281-299. Deser, C., and M. L. Blackmon, 1995: On the Relationship between Tropical and North Pacific Sea-Surface Temperature-Variations. J Climate, 8, 1677-1680. Dommenget, D., and M. Latif, 2002: A cautionary note on the interpretation of EOFs. J. Climate, 15, 216-225. Dommenget, D., V. Semenov, and M. Latif, 2006: Impacts of the tropical Indian and Atlantic Oceans on ENSO. Geophys. Res. Lett., 33(11). Dommenget, D., 2007: Evaluating EOF modes against a stochastic null hypothesis. Climate Dyn., 28, 517-531. Dommenget, D., 2011: An objective analysis of the observed spatial structure of the tropical Indian Ocean SST variability. Climate Dyn., 36, 2129-2145. Dudley, J. M., G. Genty, and B. J. Eggleton, 2008: Harnessing and control of optical rogue waves in supercontinuum generation. Opt Express, 16, 3644-3651. Eluszkiewicz, J., J. L. Moncet, M. W. Shephard, K. Cady-Pereira, T. Connor, and G. Uymin, 2008: Atmospheric and surface retrievals in the Mars polar regions from the Thermal Emission Spectrometer measurements. J Geophys Res- Planet, 113. Enfield, D. B., A. M. Mestas-Nunez, and P. J. Trimble, 2001: The Atlantic multidecadal oscillation and its relation to rainfall and river flows in the continental US. Geophys Res Lett, 28, 2077-2080. Evans, E., and Coauthors, 2004: Rise undergraduates find that regime changes in Lorenz's model are predictable. B Am Meteorol Soc, 85, 520-524. 151 Evensen, G., 1994: Sequential Data Assimilation with a Nonlinear Quasi-Geostrophic Model Using Monte-Carlo Methods to Forecast Error Statistics. J Geophys Res-Oceans, 99, 10143-10162. Forbes, J. M., 1995: Tidal and planetary waves. The Upper Mesosphere and Lower Thermosphere: A Review of Experiment and Theory, AGU, 67-87. Forbes, J. M., and M. E. Hagan, 2000: Diurnal Kelvin wave in the atmosphere of mars: Towards an understanding of 'stationary' density structures observed by the MGS accelerometer. Geophys Res Lett, 27, 3563-3566. Forbes, J. M., A. F. C. Bridger, S. W. Bougher, M. E. Hagan, J. L. Hollingsworth, G. M. Keating, and J. Murphy, 2002: Nonmigrating tides in the thermosphere of Mars. J Geophys Res-Planet, 107. Forbes, J. M., and S. Miyahara, 2006: Solar semidiurnal tide in the dusty atmosphere of Mars. J Atmos Sci, 63, 1798-1817. Forget, F., Y. Wanherdrick, and S. R. Lewis, 2001: Validation of the Mars General Circulation Model and Climate Database with new spacecraft observations, in Work Package 7, Tech. Note 11369/95/NL/JG, Eur. Space Agency, Paris. Forristall, G., 2005: Understanding rogue waves: Are new physics really necessary, in: Rogue Waves: Proc. 14th ‘Aha Huliko’a Hawaiian Winter Workshop, pp. 29-35. Frauen, C., and Dommenget, D., 2012: Influences of the tropical Indian and Atlantic Oceans on the predictability of ENSO. Geophys. Res. Lett., 39(2). Gadgil, S., P. Vinayachandran, P. Francis, and S. Gadgil, 2004: Extremes of the Indian summer monsoon rainfall, ENSO and equatorial Indian Ocean oscillation. Geophys. Res. Lett., 31. L12213, doi:10.1029/2004GL019733. Goddard, L., and N. E. Graham, 1999: Importance of the Indian Ocean for simulating rainfall anomalies over eastern and southern Africa. J Geophys Res-Atmos (1984–2012), 104, 19099-19116. 152 Gramstad, O., and K. Trulsen, 2007: Influence of crest and group length on the occurrence of freak waves. J Fluid Mech, 582, 463-472. Greybush, S. J., E. Kalnay, T. Miyoshi, K. Ide, and B. R. Hunt, 2011: Balance and Ensemble Kalman Filter Localization Techniques. Mon Weather Rev, 139, 511-522. Greybush, S. J., and Coauthors, 2012: Ensemble Kalman filter data assimilation of T hermal Emission Spectrometer temperature retrievals into a Mars GCM. J Geophys Res-Planet, 117. Greybush, S. J., E. Kalnay, M. J. Hoffman, and R. J. Wilson, 2013: Identifying Martian atmospheric instabilities and their physical origins using bred vectors. Q J Roy Meteor Soc, 139, 639-653. Gu, D. F., and S. G. H. Philander, 1995: Secular Changes of Annual and Interannual Variability in the Tropics during the Past Century. J Climate, 8, 864-876. Guan, B., and S. Nigam, 2008: Pacific sea surface temperatures in the twentieth century: An evolution-centric analysis of variability and trend. J. Climate, 21, 2790-2809. Guan, B., and S. Nigam, 2009: Analysis of Atlantic SST Variability Factoring Interbasin Links and the Secular Trend: Clarified Structure of the Atlantic Multidecadal Oscillation. J. Climate, 22, 4228-4240. Hagan, M. E., and R. G. Roble, 2001: Modeling diurnal tidal variability with the National Center for Atmospheric Research thermosphere-ionosphere- mesosphere-electrodynamics general circulation model. J Geophys Res- Space Physics, 106, 24869-24882. Hallerberg, S., D. Pazó, J. M. López, and M. A. Rodríguez, 2010: Logarithmic bred vectors in spatiotemporal chaos: Structure and growth. Phys Rev E, 81. 153 Hamill, T. M., J. S. Whitaker, and C. Snyder, 2001: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter. Mon Weather Rev, 129, 2776-2790. Hamilton, K., and R. R. Garcia, 1986: Theory and observations of the short-period normal mode oscillations of the atmosphere. J Geophys Res-Atmos, 91, 11867-11875. Hammack, J. L., and H. Segur, 1974: Korteweg-Devries Equation and Water Waves .2. Comparison with Experiments. J Fluid Mech, 65, 289-314. Hammack, J. L., and H. Segur, 1978: Korteweg-Devries Equation and Water-Waves .3. Oscillatory Waves. J Fluid Mech, 84, 337-358. Harlim, J., and A. J. Majda, 2008: Mathematical strategies for filtering complex systems: Regularly spaced sparse observations. J Comput Phys, 227, 5304- 5341. Hastenrath, S., 2002: Dipoles, temperature gradient, and tropical climate anomalies. Bull. Amer. Meteor. Soc., 83, 735-738. Henderson, K. L., D. H. Peregrine, and J. W. Dold, 1999: Unsteady water wave modulations: fully nonlinear solutions and comparison with the nonlinear Schrodinger equation. Wave Motion, 29, 341-361. Hinson, D. P., and R. J. Wilson, 2002: Transient eddies in the southern hemisphere of Mars. Geophys Res Lett, 29. Hinson, D. P., M. Pätzold, R. J. Wilson, B. Pätzold, S. Tellmann, and G. L. Tyler, 2008: Radio occultation measurements and MGCM simulations of Kelvin waves on Mars. Iracus, 193(1). Hoffman, M. J., E. Kalnay, J. A. Carton, and S. C. Yang, 2009: Use of breeding to detect and explain instabilities in the global ocean. Geophys Res Lett, 36. 154 Hoffman, M. J., J. Eluszkiewicz, D. Weisenstein, G. Uymin, and J. L. Moncet, 2012: Assessment of Mars atmospheric temperature retrievals from the Thermal Emission Spectrometer radiances. Icarus, 220, 1031-1039. Hoffman, M. J., and Coauthors, 2010: An ensemble Kalman filter data assimilation system for the martian atmosphere: Implementation and simulation experiments. Icarus, 209, 470-481. Hohmann, R., U. Kuhl, H. J. Stockmann, L. Kaplan, and E. J. Heller, 2010: Freak Waves in the Linear Regime: A Microwave Study. Phys Rev Lett, 104. Hoppel, K. W., S. D. Eckermann, L. Coy, G. E. Nedoluha, D. R. Allen, S. D. Swadley, and N. L. Baker, 2013: Evaluation of SSMIS Upper Atmosphere Sounding Channels for High-Altitude Data Assimilation. Mon Weather Rev, 141, 3314-3330. Horii, T., H. Hase, I. Ueki, and Y. Masumoto, 2008: Oceanic precondition and evolution of the 2006 Indian Ocean dipole. Geophys. Res. Lett., 35, L03607, doi:10.1029/2007GL032464. Houben, H., 1999: Assimilation of Mars Global Surveyor meteorological data. Moon and Mars, 23, 1899-1902. Houtekamer, P. L., and H. L. Mitchell, 2001: A sequential ensemble Kalman filter for atmospheric data assimilation. Mon Weather Rev, 129, 123-137. Huang, B., and J. Shukla, 2007a: Mechanisms for the Interannual Variability in the Tropical Indian Ocean. Part I: The Role of Remote Forcing from the Tropical Pacific. J. Climate, 20, 2917-2936. Huang, B., and J. Shukla, 2007b: Mechanisms for the Interannual Variability in the Tropical Indian Ocean. Part II: Regional Processes. J. Climate, 20, 2937- 2960. 155 Hunt, B. R., E. J. Kostelich, and I. Szunyogh, 2007: Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. Physica D, 230, 112-126. Islas, A., and C. M. Schober, 2011: Rogue waves, dissipation, and downshifting. Physica D, 240, 1041-1054. Izumo, T., and Coauthors, 2010: Influence of the state of the Indian Ocean Dipole on the following year's El Nino. Nat Geosci, 3, 168-172. Jansen, M. F., D. Dommenget, and N. Keenlyside, 2009: Tropical atmosphere–ocean interactions in a conceptual framework. J. Climate, 22, 550-567. Kalnay, E., 2003: Atmospheric modeling, data assimilation, and predictability, Cambridge university press. Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-year reanalysis project. B Am Meteorol Soc, 77, 437-471. Kalnay, E., and S. C. Yang, 2010: Notes and Correspondence Accelerating the spin- up of Ensemble Kalman Filtering. Q J Roy Meteor Soc, 136, 1644-1651. Kalnay, E., H. Li, T. Miyoshi, S. C. Yang, and J. Ballabrera-Poy, 2007: 4-D-Var or ensemble Kalman filter? Tellus A, 59, 758-773. Kaplan, A., Y. Kushnir, and M. A. Cane, 2000: Reduced space optimal interpolation of historical marine sea level pressure: 1854-1992. J Climate, 13, 2987-3002. Kassam, A. K., and L. N. Trefethen, 2005: Fourth-order time-stepping for stiff PDEs. Siam J Sci Comput, 26, 1214-1233. Kawamura, R., 1994: A Rotated Eof Analysis of Global Sea-Surface Temperature Variability with Interannual and Interdecadal Scales. J Phys Oceanogr, 24, 707-715. Keller, J. D., A. Hense, L. Kornblueh, and A. Rhodin, 2010: On the Orthogonalization of Bred Vectors. Weather Forecast, 25, 1219-1234. 156 Kharif, C., and E. Pelinovsky, 2003: Physical mechanisms of the rogue wave phenomenon. Eur J Mech B-Fluid, 22, 603-634. Kharif, C., E. Pelinovsky, T. Talipova, and A. Slunyaev, 2001: Focusing of nonlinear wave groups in deep water. Jetp Lett+, 73, 170-175. Kibler, B., and Coauthors, 2010: The Peregrine soliton in nonlinear fibre optics. Nat Phys, 6, 790-795. Kirkpatrick, K., B. Schlein, and G. Staffilani, 2011: Derivation of the Two- Dimensional Nonlinear Schrodinger Equation from Many Body Quantum Dynamics. Am J Math, 133, 91-130. Kleinbohl, A., R. J. Wilson, D. Kass, J. T. Schofield, and D. J. McCleese, 2013: The semidiurnal tide in the middle atmosphere of Mars. Geophys Res Lett, 40, 1952-1959. Korteweg, D.J., G. De Vries, 1895: XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 39, 422-443. Kug, J.-S., and I.-S. Kang, 2006: Interactive feedback between ENSO and the Indian Ocean. J. Climate, 19, 1784-1801. Latif, M., D. Dommenget, M. Dima, and A. Grötzner, 1999: The role of Indian Ocean sea surface temperature in forcing east African rainfall anomalies during December-January 1997/98. J. Climate, 12, 3497-3504. Lee, C. C. W., J. Savarino, and M. H. Thiemens, 2001: Mass independent oxygen isotopic composition of atmospheric sulfate: Origin and implications for the present and past atmosphere of Earth and mars. Geophys Res Lett, 28, 1783- 1786. Lee, C., W. G. Lawson, M. I. Richardson, J. L. Anderson, N. Collins, T. Hoar, and M. Mischna, 2011: Demonstration of ensemble data assimilation for Mars using 157 DART, MarsWRF, and radiance observations from MGS TES. J Geophys Res-Planet, 116. Levitus, S., J. I. Antonov, T. P. Boyer, H. E. Garcia, and R. A. Locarnini, 2005: EOF analysis of upper ocean heat content, 1956-2003. Geophys Res Lett, 32. Lewis, S. R., and P. R. Barker, 2005: Atmospheric tides in a Mars general circulation model with data assimilation. Planetary Atmospheres, Ionospheres, and Magnetospheres, 36, 2162-2168. Lewis, S. R., P. L. Read, B. J. Conrath, J. C. Pearl, and M. D. Smith, 2007: Assimilation of thermal emission spectrometer atmospheric data during the Mars Global Surveyor aerobraking period. Icarus, 192, 327-347. Li, H., E. Kalnay, T. Miyoshi, and C. M. Danforth, 2009: Accounting for Model Errors in Ensemble Data Assimilation. Mon Weather Rev, 137, 3407-3419. Lin, S. J., 2004: A "vertically Lagrangian'' finite-volume dynamical core for global models. Mon Weather Rev, 132, 2293-2307. Lindzen, R. S., 1978: Effect of daily variations of cumulonimbus activity on the atmospheric semidiurnal tide. Mon Weather Rev, 106, 526-533. Lindzen, R. S., 1990: Dynamics in Atmospheric Physics, Cambridge University Press, New York, 310pp. Lindzen, R. S., and S. Chapman, 1969: Atmospheric tides. Sp. Sci. Revs., 10, 3-188. Longuet-Higgins, M. S., 1952: On the Statistical Distribution of the Heights of Sea Waves. J Mar Res, 11, 245-266. Lorenc, A. C., R. S. Bell, and B. Macpherson, 1991: The Meteorological-Office Analysis Correction Data Assimilation Scheme. Q J Roy Meteor Soc, 117, 59- 89. Lorenz, E. N., 1963: Deterministic Nonperiodic Flow. J Atmos Sci, 20, 130-141. Lorenz, E. N., 2006: Computational periodicity as observed in a simple system. Tellus A, 58, 549-557. 158 Luo, J. J., S. Behera, Y. Masumoto, H. Sakuma, and T. Yamagata, 2008: Successful prediction of the consecutive IOD in 2006 and 2007. Geophys. Res. Lett., 35, L14S02. Luo, J. J., R. Zhang, S. K. Behera, Y. Masumoto, F.-F. Jin, R. Lukas, and T. Yamagata, 2010: Interaction between El Nino and Extreme Indian Ocean Dipole. J. Climate, 23, 726-742. Madeleine, J. B., F. Forget, E. Millour, T. Navarro, and A. Spiga, 2012: Interaction between El Nino and Extreme Indian Ocean Dipole. Geophys Res Lett, 39, L23202. Magnusson, L., M. Leutbecher, and E. Kallen, 2008: Comparison between Singular Vectors and Breeding Vectors as Initial Perturbations for the ECMWF Ensemble Prediction System. Mon Weather Rev, 136, 4092-4104. Malomed, B. A., D. Mihalache, F. Wise, and L. Torner, 2005: Spatiotemporal optical solitons. J Opt B-Quantum S O, 7, R53-R72. Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R. C. Francis, 1997: A Pacific interdecadal climate oscillation with impacts on salmon production. Bull. Amer. Meteor. Soc., 78, 1069-1079. Mayo, N., 1997: Ocean Waves-Their Energy and Power, Physics Teacher, 35, 352. Mclaughlin, D. W., and C. M. Schober, 1992: Chaotic and Homoclinic Behavior for Numerical Discretizations of the Nonlinear Schrodinger-Equation. Physica D, 57, 447-465. McCleese, D. J., and Coauthors, 2007: Mars Climate Sounder: An investigation of thermal and water vapor structure, dust and condensate distributions in the atmosphere, and energy balance of the polar regions. J Geophys Res-Planet, 112. 159 Meehl, G. A., J. M. Arblaster, and J. Loschnigg, 2003: Coupled ocean-atmosphere dynamical processes in the tropical Indian and Pacific Oceans and the TBO. J. Climate, 16, 2138-2158. Merle, J., 1980: Seasonal heat budget in the equatorial Atlantic Ocean. J. Phys. Oceanogr., 10, 464-469. Mestas-Nunez, A. M., and D. B. Enfield, 1999: Rotated global modes of non-ENSO sea surface temperature variability. J Climate, 12, 2734-2746. Meyers, G., P. McIntosh, L. Pigot, and M. Pook, 2007: The years of El Nino, La Nina, and interactions with the tropical Indian Ocean. J. Climate, 20, 2872- 2880. Minobe, S., 1997: A 50-70 year climatic oscillation over the North Pacific and North America. Geophys Res Lett, 24, 683-686. Monahan, A. H., 2001: Nonlinear principal component analysis: Tropical Indo- Pacific sea surface temperature and sea level pressure. J Climate, 14, 219-233. Mitchell, D. M., L. Montabone, S. Thomson, and P. L. Read, 2014: Polar vortices on Earth and Mars: A comparative study of the climatology and variability from reanalyses. Q J Roy Meteor Soc, n/a-n/a. Miyoshi, T., 2011: The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter. Mon Weather Rev, 139, 1519-1535. Miyoshi, T., and S. Yamane, 2007: Local ensemble transform Kalman filtering with an AGCM at a T159/L48 resolution. Mon Weather Rev, 135, 3841-3861. Miyoshi, T., Y. Sato, and T. Kadowaki, 2010: Ensemble Kalman Filter and 4D-Var Intercomparison with the Japanese Operational Global Analysis and Prediction System. Mon Weather Rev, 138, 2846-2866. Montabone, L., S. R. Lewis, and P. L. Read, 2005: Interannual variability of Martian dust storms in assimilation of several years of Mars global surveyor 160 observations. Planetary Atmospheres, Ionospheres, and Magnetospheres, 36, 2146-2155. Montabone, L., S. R. Lewis, and P. L. Read, 2011: Mars Analysis Correction Data Assimilation (MACDA): MGS/TES v1.0. NCAS British Atmospheric Data Centre, doi:10.5285/78114093-E2BD-4601-8AE5-3551E62AEF2B. Montabone, L., S. R. Lewis, L. Steele, P. L. Read, T. Ruan, M. D. Smith, D. Kass, A. Kleinböhl, J. T. Schofield, J. H. Shirley, and D. J. McCleese, 2012: Mars analysis correction data assimilation: a multi-annual reanalysis of atmospheric observations for the red planet. In: 4th World Climate Research Programme International Conference on Reanalyses, Silver Spring, Maryland, USA. Montabone, L. et al., 2014: The Mars Analysis Correction Data Assimilation (MACDA) Dataset V1.0. Geoscience Data Journal. DOI: 10.1002/gdj3.13. Montmessin, F., F. Forget, P. Rannou, M. Cabane, and R. M. Haberle, 2004: Origin and role of water ice clouds in the Martian water cycle as inferred from a general circulation model. J Geophys Res-Planet, 109. Murtugudde, R., J. P. McCreary, and A. J. Busalacchi, 2000: Oceanic processes associated with anomalous events in the Indian Ocean with relevance to 1997- 1998. J Geophys Res-Oceans (1978–2012), 105, 3295-3306. Nakamura, H., G. Lin, and T. Yamagata, 1997: Decadal climate variability in the North Pacific during the recent decades. B Am Meteorol Soc, 78, 2215-2225. Nayfeh, A.H., and B. Balachandran, 1995: Chaos, in: Applied Nonlinear Dynamics, Wiley-VCH Verlag GmbH, pp. 277-421. Newman, C. E., P. L. Read, and S. R. Lewis, 2004: Investigating atmospheric predictability on Mars using breeding vectors in a general-circulation model. Q J Roy Meteor Soc, 130, 2971-2989. Nicholls, N., and W. Drosdowsky, 2001: Is there an equatorial Indian Ocean SST dipole, independent of the El Niño Southern Oscillation? 81st American 161 Meteorological Society Annual Meeting, Albuquerque, New Mexico, USA, 14-19 January 2001. Nigam, S., and H.-S. Shen, 1993: Structure of Oceanic and Atmospheric Low- Frequency Variability over the Tropical Pacific and Indian Oceans. Part I: COADS observations. J. Climate, 6, 657-676. Nigam, S., M. Barlow, and E. H. Berbery, 1999: Analysis links Pacific decadal variability to drought and streamflow in United States. Eos, Trans. Amer. Geophys. Union, 80, 621–625. Nigam, S., 2003: Teleconnections. Encyclopedia of Atmospheric Sciences, J. R. Holton et al., Eds., Academic Press, 2243–2269. Parker, D., C. Folland, A. Scaife, J. Knight, A. Colman, P. Baines, and B. W. Dong, 2007: Decadal to multidecadal variability and the climate change background. J Geophys Res-Atmos, 112. Nore, C., M. E. Brachet, and S. Fauve, 1993: Numerical Study of Hydrodynamics Using the Nonlinear Schrodinger-Equation. Physica D, 65, 154-162. Onorato, M., A. R. Osborne, M. Serio, and S. Bertone, 2001: Freak waves in random oceanic sea states. Phys Rev Lett, 86, 5831-5834. Osborne, A. R., and T. L. Burch, 1980: Internal Solitons in the Andaman Sea. Science, 208, 451-460. Osborne, A. R., M. Onorato, and M. Serio, 2000: The nonlinear dynamics of rogue waves and holes in deep-water gravity wave trains. Phys Lett A, 275, 386- 393. Pazó, D., M. A. Rodríguez, and J. M. López, 2010: Spatio-temporal evolution of perturbations in ensembles initialized by bred, Lyapunov and singular vectors. Tellus A, 62, 10-23. 162 Pazó, D., M. A. Rodríguez, and J. M. López, 2011: Maximizing the Statistical Diversity of an Ensemble of Bred Vectors by Using the Geometric Norm. J Atmos Sci, 68, 1507-1512. Pazó, D., J. M. López, and M. A. Rodríguez, 2013: The geometric norm improves ensemble forecasting with the breeding method. Q J Roy Meteor Soc, 139, 2021-2032. Pelinovsky, E., T. Talipova, and C. Kharif, 2000: Nonlinear-dispersive mechanism of the freak wave formation in shallow water. Physica D, 147, 83-94. Peña, M., and E. Kalnay, 2004: Separating fast and slow modes in coupled chaotic systems. Nonlinear Proc Geoph, 11, 319-327. Philander, S. G., 2003: Of Dipoles and Spherical Cows. Bull. Amer. Meteor. Soc., 84, 1424-1424. Pires, C., R. Vautard, and O. Talagrand, 1996: On extending the limits of variational assimilation in nonlinear chaotic systems. Tellus A, 48, 96-121. Polavarapu, S., T. G. Shepherd, Y. Rochon, and S. Ren, 2005: Some challenges of middle atmosphere data assimilation. Q J Roy Meteor Soc, 131, 3513-3527. Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, 1992: Numerical Recipes in C: The Art of Scientific Computing, Cambridge Univ. Press. Primo, C., M. A. Rodríguez, and J. M. Gutiérrez, 2008: Logarithmic bred vectors. A new ensemble method with adjustable spread and calibration time. J Geophys Res-Atmos, 113. Rasmusson, E. M., X. Wang, and C. F. Ropelewski, 1990: The biennial component of ENSO variability. J. Mar. Syst., 1, 71-96. Rasmusson, E. M., and T. H. Carpenter, 1982: Variations in Tropical Sea-Surface Temperature and Surface Wind Fields Associated with the Southern Oscillation El-Nino. Mon Weather Rev, 110, 354-384. 163 Rayner, N. A., and Coauthors, 2003: Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J Geophys Res-Atmos (1984–2012), 108. Reason, C., R. Allan, J. Lindesay, and T. Ansell, 2000: ENSO and climatic signals across the Indian Ocean Basin in the global context: Part I, Int. J. Climatol., 20, 1285-1327. Reason, C., 2001: Subtropical Indian Ocean SST dipole events and southern African rainfall. Geophys. Res. Lett., 28, 2225-2227. Richardson, M. I., A. D. Toigo, and C. E. Newman, 2007: PlanetWRF: A general purpose, local to global numerical model for planetary atmospheric and climate dynamics. J Geophys Res-Planet, 112. Rogberg, P., P. L. Read, S. R. Lewis, and L. Montabone, 2010: Assessing atmospheric predictability on Mars using numerical weather prediction and data assimilation. Q J Roy Meteor Soc, 136, 1614-1635. Rosenthal, W., and S. Lehner, 2005: Results of the MAXWAVE project, in: Proc. 14th Aha Huliko Winter Workshop, Honolulu, Hawaii. Ruiz-Barradas, A., J. A. Carton, and S. Nigam, 2000: Structure of interannual-to- decadal climate variability in the tropical Atlantic sector. J Climate, 13, 3285- 3297. Russell, J. S.,1844: Report on waves, in: 14th meeting of the British Association for the Advancement of Science, pp. 390. Saji, N. H., B. N. Goswami, P. N. Vinayachandran, and T. Yamagata, 1999: A dipole mode in the tropical Indian Ocean. Nature, 401, 360-363. Sakazaki, T., M. Fujiwara, X. Zhang, M. E. Hagan, and J. M. Forbes, 2012: Diurnal tides from the troposphere to the lower mesosphere as deduced from TIMED/SABER satellite data and six global reanalysis data sets. J Geophys Res-Atmos, 117, D13108. 164 Sankey, D., S. Ren, S. Polavarapu, Y. J. Rochon, Y. Nezlin, and S. Beagley, 2007: Impact of data assimilation filtering methods on the mesosphere. J Geophys Res-Atmos, 112, D24104. Segur, H., and J. L. Hammack, 1982: Soliton Models of Long Internal Waves. J Fluid Mech, 118, 285-304. Singleton, T. 2011: Data assimilation experiments with a simple coupled ocean- atmosphere model. Ph. D. thesis, University of Maryland, College Park, 116 pp. Skourup, J., N. E. O. Hansen, and K. K. Andreasen, 1997: Non-Gaussian extreme waves in the Central North Sea. J Offshore Mech Arct, 119, 146-150. Smith, M. D., J. C. Pearl, B. J. Conrath, and P. R. Christensen, 2001: One Martian year of atmospheric observations by the Thermal Emission Spectrometer. Geophys Res Lett, 28, 4263-4266. Solli, D. R., C. Ropers, P. Koonath, and B. Jalali, 2007: Optical rogue waves. Nature, 450, 1054-U1057. Som, B. K., M. R. Gupta, and B. Dasgupta, 1979: Coupled Non-Linear Schrodinger Equation for Langmuir and Dispersive Ion-Acoustic Waves. Phys Lett A, 72, 111-114. Sparrow, M. K., 1982: Digital Coding of Single Fingerprints - a New Approach for the Computer-Age. J Police Sci Admin, 10, 206-217. Stansell, P., 2005: Distributions of extreme wave, crest and trough heights measured in the North Sea. Ocean Eng, 32, 1015-1036. Strang, G., 1968: On Construction and Comparison of Difference Schemes. Siam J Numer Anal, 5, 506-&. Swinbank, R., R. L. Orris, and D. L. Wu, 1999: Stratospheric tides and data assimilation. J Geophys Res-Atmos, 104, 16929-16941. 165 Szunyogh, I., and Coauthors, 2008: A local ensemble transform Kalman filter data assimilation system for the NCEP global model. Tellus A, 60, 113-130.Tang, Y. M., and Z. W. Deng, 2010: Low-dimensional nonlinearity of ENSO and its impact on predictability. Physica D, 239, 258-268. Tang, Y. M., and Z. W. Deng, 2011: Bred Vector and ENSO Predictability in a Hybrid Coupled Model during the Period 1881-2000. J Climate, 24, 298-314. Toth, Z., and E. Kalnay, 1993: Ensemble Forecasting at Nmc - the Generation of Perturbations. B Am Meteorol Soc, 74, 2317-2330. Toth, Z., and E. Kalnay, 1997: Ensemble forecasting at NCEP and the breeding method. Mon Weather Rev, 125, 3297-3319. Tsonis, A. A., and J. B. Elsner, 1992: Nonlinear Prediction as a Way of Distinguishing Chaos from Random Fractal Sequences. Nature, 358, 217-220. Uppala, S. M., and Coauthors, 2005: The ERA‐40 re‐analysis. Quarterly Q. J. R. Meteorol. Soc., 131(612), 2961-3012. Visbeck, M. H., J. W. Hurrell, L. Polvani, and H. M. Cullen, 2001: The North Atlantic Oscillation: Past, present, and future. P Natl Acad Sci USA, 98, 12876-12877. Vitanov, N. K., A. Chabchoub, and N. Hoffmann, 2013: Deep-Water Waves: on the Nonlinear Schrodinger Equation and its Solutions, Journal of Theoretical and Applied Mechanics, 43, 43-54. arXiv:1301.0990, doi:10.2478/jtam-2013- 0013. Weare, B. C., and J. S. Nasstrom, 1982: Examples of Extended Empirical Orthogonal Function Analyses. Mon Weather Rev, 110, 481-485. Webster, P. J., A. M. Moore, J. P. Loschnigg, and R. R. Leben, 1999: Coupled ocean- atmosphere dynamics in the Indian Ocean during 1997-98. Nature, 401, 356- 360. 166 Whitaker, J. S., T. M. Hamill, X. Wei, Y. C. Song, and Z. Toth, 2008: Ensemble data assimilation with the NCEP Global Forecast System. Mon Weather Rev, 136, 463-482. White, B. S., and B. Fornberg, 1998: On the chance of freak waves at sea. J Fluid Mech, 355, 113-138. Wilson, R. J., 2000: Evidence for diurnal period Kelvin waves in the Martian atmosphere from Mars Global Surveyor TES data. Geophys Res Lett, 27, 3889-3892. Wilson, R. J., 2002: Evidence for nonmigrating thermal tides in the Mars upper atmosphere from the Mars Global Surveyor Accelerometer Experiment. Geophys Res Lett, 29. Wilson, R. J. 2011: Dust cycle modeling with the GFDL Mars general circulation model, 4th International workshop on the Mars atmosphere: modeling and observations, Paris. Wilson, R. J., and K. Hamilton, 1996: Comprehensive model simulation of thermal tides in the Martian atmosphere. J Atmos Sci, 53, 1290-1326. Wilson, R. J., S. R. Lewis, L. Montabone, and M. D. Smith, 2008: Influence of water ice clouds on Martian tropical atmospheric temperatures. Geophys Res Lett, 35. Wilson, R. J., and M. I. Richardson, 2000: The Martian atmosphere during the viking mission, I - Infrared measurements of atmospheric temperatures revisited. Icarus, 145, 555-579. Wilson, R J., and S. D. Guzewich, in press: Influence of Water Ice Clouds on Nighttime Tropical Temperature Structure as Seen by the Mars Climate Sounder. Geophys Res Lett. doi:10.1002/2014GL060086. 5/14. 167 Wu, L., Z. Liu, R. Gallimore, R. Jacob, D. Lee, and Y. Zhong, 2003: Pacific decadal variability: The tropical Pacific mode and the North Pacific mode. J Climate, 16, 1101-1120. Wu, R., and B. P. Kirtman, 2004: Understanding the impacts of the Indian Ocean on ENSO variability in a coupled GCM. J. Climate, 17, 4019-4031. Wyrtki, K., 1973: Physical oceanography of the Indian Ocean. The biology of the Indian Ocean, Springer, 18-36. Xiao, W. 2013: Study of Directional Ocean Wavefield Evolution and Rogue Wave Occurrence Using Large-Scale Phase-Resolved Nonlinear Simulations. Ph. D. thesis, Massachusetts Institute of Technology, 134 pp. Xie, P., and P. A. Arkin, 1996: Analyses of global monthly precipitation using gauge observations, satellite estimates, and numerical model predictions. J. Climate, 9, 840-858. Xie, S.-P., H. Annamalai, F. A. Schott, and J. P. McCreary, 2002: Structure and mechanisms of South Indian Ocean climate variability. J. Climate, 15, 864- 878. Yamagata. T, S. K. Behera, S. A. Rao, Z. Guan, K. Ashok, and H. N. Saji, 2003: Comments on “Dipoles, temperature gradients, and tropical climate anomalies”. Bull. Amer. Meteor. Soc., 84, 1418-1422. Yamagata. T, S. K. Behera, J. J. Luo, S. Masson, M. R. Jury, and S. A. Rao, 2004: Coupled Ocean-Atmosphere Variability in the Tropical Indian Ocean. AGU Book Ocean-Atmosphere Interaction and Climate Variability, C. Wang, S.-P. Xie and J.A. Carton (eds.), Geophys. Monogr., 147, AGU, Washington D.C., 189-212. Yang, S. C., E. Kalnay, M. Cai, and M. M. Rienecker, 2008: Bred vectors and tropical Pacific forecast errors in the NASA coupled general circulation model. Mon Weather Rev, 136, 1305-1326. 168 Yang, S. C., C. Keppenne, M. Rienecker, and E. Kalnay, 2009: Application of Coupled Bred Vectors to Seasonal-to-Interannual Forecasting and Ocean Data Assimilation. J Climate, 22, 2850-2870. Yang, S. C., M. Cai, E. Kalnay, M. Rienecker, G. Yuan, and Z. Toth, 2006: ENSO bred vectors in coupled ocean-atmosphere general circulation models. J Climate, 19, 1422-1436. Yu, J. Y., and K. M. Lau, 2005: Contrasting Indian Ocean SST variability with and without ENSO influence: A coupled atmosphere-ocean GCM study. Meteorol. Atmos. Phys., 90, 179-191. Zabusky, N. J., 1984: Computational Synergetics. Phys Today, 37, 36-45. Zabusky, N. J., and M. D. Kruskal, 1965: Interaction of Solitons in a Collisionless Plasma and Recurrence of Initial States. Phys Rev Lett, 15, 240-&. Zhang, K. Q., A. P. Ingersoll, D. M. Kass, J. C. Pearl, M. D. Smith, B. J. Conrath, and R. M. Haberle, 2001: Assimilation of Mars Global Surveyor atmospheric temperature data into a general circulation model. J Geophys Res-Planet, 106, 32863-32877. Zhang, R. H., and A. J. Busalacchi, 2005: Interdecadal change in properties of El Nino-Southern Oscillation in an intermediate coupled model. J Climate, 18, 1369-1380. Zurek, R. W., 1976: Diurnal Tide in Martian Atmosphere. J Atmos Sci, 33, 321-337. Zurek, R. W., 1988: Free and Forced Modes in the Martian Atmosphere. J Geophys Res-Atmos, 93, 9452-9462. Zurek, R. W., J. R. Barnes, R. M. Haberle, J. B. Pollack, J. E. Tillman, and C. B. Leovy. 1992: Dynamics of the atmosphere of Mars. Mars, ed. H. Kieffer et al., U. Ariz. Press, Tucson, 835-933.