UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Enso simulation and prediction using hybrid coupled models with data assimilation Tang, Youmin 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2001-611825.pdf [ 8.38MB ]
Metadata
JSON: 831-1.0053123.json
JSON-LD: 831-1.0053123-ld.json
RDF/XML (Pretty): 831-1.0053123-rdf.xml
RDF/JSON: 831-1.0053123-rdf.json
Turtle: 831-1.0053123-turtle.txt
N-Triples: 831-1.0053123-rdf-ntriples.txt
Original Record: 831-1.0053123-source.json
Full Text
831-1.0053123-fulltext.txt
Citation
831-1.0053123.ris

Full Text

E N S O S I M U L A T I O N A N D P R E D I C T I O N USING H Y B R I D C O U P L E D M O D E L S W I T H D A T A ASSIMILATION by Y o u m i n T a n g B . Sc. (Meteoro logy) N a n j i n g Ins t i tu te of Meteoro logy , 1984 M . Sc. (Meteoro logy) N a n j i n g Ins t i tu te of Me teo ro logy , 1987 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y i n T H E F A C U L T Y O F G R A D U A T E S T U D I E S Department of E A R T H A N D O C E A N S C I E N C E S W e accept this thesis as conforming to the requi red s t andard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A A p r i l 2001 © Y o u m i n T a n g , 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1Z1 Date: Abstrac t The possibility of using a nonlinear empirical atmospheric model for hybrid coupled atmosphere-ocean modelling has been examined by using a neural network (NN) model for predicting the contemporaneous wind stress field from the upper ocean state in the tropical Pacific. Upper ocean heat content (HC) from a 6-layer dynamic ocean model was a better predictor of the wind stress than the (observed or modelled) sea surface temperature (SST). The results showed that the N N model generally had slightly better skills in predicting the contemporaneous wind stress than the linear regression (LR) model in the off-equatorial tropical Pacific and in the eastern equatorial Pacific. Next, the N N and LR atmospheric models were separately coupled to the dynamic ocean model, yielding respectively a hybrid coupled model with a nonlinear atmosphere (NHCM) and one with a linear atmosphere (LHCM). The POP (Principal Oscillation Pattern) analysis shows that the N H C M can produce more realistic ENSO oscillatory behaviour than the L H C M . The phase-locking between the peak of the warm events and the seasonal cycle is more realistically distributed in the N H C M than the very narrow phase-locking in the L H C M . Sensitivity experiments show that in the absence of external forcing, neither N H C M nor L H C M displays the irregular behaviour of ENSO oscillations, suggesting that stochastic forcing (instead of chaos) is likely to cause the irregularity of ENSO. ENSO prediction skills in the two hybrid coupled models have also been investigated under two initialization schemes. The stress-from-ocean initialization scheme, which considers the ocean-feedback in the initial conditions, generally has better predictive skills than the stress-only scheme, where the ocean is simply forced by the observed ii wind stress. The stress-from-ocean scheme also manifests less decadal variability in the forecast skills than the stress-only scheme. From 1964-1998, with the stress-from-ocean initialization scheme, the forecast correlation skill at 12-month lead time for the NIN03 SST anomalies (SSTA) is about 0.55 for the N H C M and 0.50-0.55 for the L H C M . The main difference in forecast skills between the N H C M and the L H C M occurs in the 1990s, where the N H C M has better skills. A nonlinear canonical correlation analysis of the wind stress and the SSTA shows that the relation between the two was indeed more nonlinear in the 1990s than in the 1980s, which would give the N H C M an advantage over the L H C M in the 1990s. The impact of assimilating different types of data on ENSO simulations and predic-tions was investigated by separately assimilating the SST, sea surface height anomalies (SSHA), upper ocean heat content anomalies (HCA), and wind stress, with the 3-D Var (variational assimilation) technique. For ENSO prediction, H C A assimilation is the best in improving the correlation skills of the prediction for lead times greater than 5 months. The improvement from SST assimilation is the best for lead times of 4 months or shorter; at longer lead times, SST assimilation degrades the predictions. Wind-stress assimilation is generally the second best for lead times 6 months or longer, while SSHA assimilation generally produces prediction skills only slightly better than without data assimilation. In general, data assimilation yielded less significant improvements to ENSO prediction in the 1990s than in 1980s, which was explained upon examining the H C A in the western and in the eastern equatorial Pacific. In summary, this thesis has initiated the fusion of neural network techniques into dynamical models. Using an N N model for the atmosphere, it has produced the first H C M with a nonlinear empirical atmospheric component, and showed that the nonlinear atmosphere could have advantages over a linear atmosphere in modelling ENSO vari-ability and in ENSO prediction. This study has introduced N N for the assimilation of in non-prognostic variables (e.g. HCA) by using N N to relate the non-prognostic variable to prognostic variables, which are then assimilated into the model. While the full 4-D Var H C M is beyond the scope of this thesis, a neural-dynamical hybrid approach under 4-D Var has been developed to study the simple Lorenz system. iv Table of Contents Abstract ! ii List of Tables ix List of Figures xi List of Acronyms xxv Preface xxvii Acknowledgements xxix 1 Introduction 1 1.1 ENSO theory 1 1.1.1 ENSO and its characteristics 1 1.1.2 Bjerknes hypothesis 2 1.1.3 Delayed-Oscillator mechanism 4 1.1.4 ENSO chaos and stochastic forcing . 5 1.2 ENSO models and ENSO prediction 6 .1.3 Data assimilation in ENSO prediction 8 1.4 Nonlinear physical processes in the tropical atmosphere 10 1.4.1 Physical processes associated with SST in the tropical atmosphere 10 1.4.2 Nonlinearities and linearization 12 1.5 Objectives of the thesis 13 2 The construction of two hybrid models 17 2.1 Introduction 17 2.2 Ocean Model 18 2.2.1 The modelling strategy 18 2.2.2 The modelling skills 20 2.3 Atmospheric models 22 2.3.1 EOFs for predictors and predictands 22 2.3.2 Neural Network Model 27 2.3.3 Linear model 30 2.4 Results from atmospheric models 30 2.4.1 The predictors 30 2.4.2 Prediction skills of the N N and L R models 33 2.5 The ocean model driven by the empirical wind 39 2.5.1 SST skills 39 2.5.2 Heat content Skills 45 3 Hybrid coupled models of the tropical Pacific — Interannual variability 50 3.1 Introduction 50 3.2 The Coupling strategy 51 3.3 Simulations from the standard coupling experiments 53 3.3.1 N H C M simulations 53 3.3.2 Some differences between N H C M and L H C M 61 3.4 Sensitivity experiments 69 4 Hybrid coupled models of the tropical Pacific — E N S O prediction 74 4.1 Introduction 74 4.2 The strategy of prediction experiments 74 v i 4.3 Prediction experiments with 'stress-only' initializations 76 4.4 Prediction experiments with 'stress-from-ocean' initialization 91 4.4.1 Case study 92 5 Impact of data assimilation on E N S O simulations and predictions 103 5.1 Introduction 103 5.2 The assimilation system 104 5.3 Methodology for assimilating different types of data 107 5.3.1 Data 107 5.3.2 SSHA observations 107 5.3.3 H C A observations .; •'. 108 5.3.4 Assimilation strategies 114 5.4 Impact of data assimilation on ENSO simulation 114 5.5 The impact of data assimilation on ENSO prediction 118 6 Cons t ruc t ing hybr id coupled models v ia 4-D variational data assimila-t ion 130 6.1 Introduction 130 6.2 Hybrid neural-dynamical Lorenz model 131 6.3 Simple hybrid model experiments 132 6.4 The hybrid model using variational data assimilation 136 6.5 Determining N N parameters via variational assimilation 138 6.5.1 The initial guesses 138 6.5.2 Impact of assimilation window in the highly nonlinear case . . . 142 6.5.3 SC assimilation 145 6.6 Determining dynamical parameters and initial conditions by the hybrid model 156 vii 7 Summary and Discussion 162 8 References 171 Appendices 185 A Ocean model formulation 185 B Discrete form of Lorenz hybrid model for coding 188 C Summary of Experiments 193 v m List of Tables 2.1 Correlation between the observed SST and that from 10 models. The results of the first 9 models are taken from Palmer and Anderson (1994). The final model is the one used in this study. Model A , described in Wu described in Davey et al(1994), is also a 2|-layer model. While model D and E are versions of the G F D L Modular Ocean Model, with resolution of l | ° x l | ° and | ° latitude x l | ° longitude respectively 21 2.2 The contributions of the first 10 principal components in % 25 2.3 Cross-validated correlation between the predicted wind stress E O F time series and the observed wind stress E O F time series for the first 3 modes, using the observed SST, the model SST, the model H C , and the model HC+SST as predictors. Results are shown for both the N N model and the L R model, and for the zonal (a;) and meridional (y) components of the wind stress 31 5.1 Information on the datasets used 107 7.1 Correlation between the observed SST NIN03 index and that from Lamont model, described in Chen et al (1995, 1998) and the N H C M with HC assimilation 169 A . I Values of the parameters used in the 6-layer ocean model 187 ix Summary of main experiments and results List of Figures 1.1 E O F time series for SST and zonal wind stresses. Solid line is for wind stresses and dashed line for SST. Upper: E 0 F 1 ; middle: E 0 F 2 ; bottom: E 0 F 3 . A l l time series are smoothed by a 5-month running mean filter. . 14 2.1 Location of the various oceanic regions used in the analysis 20 2.2 Statistics for the model SST anomalies relative to the COADS SST data:(a) correlation (contour interval=c.i.=0.1) and (b) RMS error (c.i.=0.1°C). . 23 2.3 Time evolution of SST anomalies from the ocean model forced with F S U observed stress (solid line) and observations (dash line) in regions EQ2, NIN03 and NIN04 from top to bottom. The Y-axis units are ° C 24 2.4 The time series of the first E O F mode, the second E O F mode and third E O F mode from top to bottom. The left panel is for zonal wind and right panel for meridional wind 26 2.5 A n example of a neural network model, where there are four neurons in the input layer, three in the hidden layer, and one in the output layer. The parameters and Wj are the weights, and bj and b are the biases. The parameters bj and b can also be regarded as the weights for constant inputs of value 1 28 2.6 The first E O F mode for (a) the heat content (HC) from the ocean model driven by the F S U wind stress, (b) the F S U zonal wind stress and (c) the meridional wind stress. Negative contours are shown as dashed curves and zero contours as dotted curves. Contour interval = 0.01 32 xi 2.7 (a) The first mode of the zonal wind stress plotted as a function of the third mode of the SST, and (b) the second mode of the zonal wind stress versus the first mode of the SST. The observations are shown as dots; the simulated zonal wind stress mode from the N N model is indicated by the small circles and from the L R model by the straight lines 34 2.8 Correlations of the FSU wind stress and the idealized wind stress (recon-structed from the first 3 E O F modes of the F S U wind stress) for (a) the zonal stress and (b) the meridional stress. Contour interval = 0.1 35 2.9 Cross-validated anomaly correlation between the predicted wind stress by N N (with model HC as predictors) and the idealized wind stress (i.e. the F S U wind stress with only the first 3 E O F modes): (a) zonal stress, and (b) meridional stress. Contour interval=0.1. The seasonal cycle has been removed prior to the analysis, and in all following maps 37 2.10 The sum of squared error (SSE) of the predicted wind stress by N N verified against the idealized wind stress: (a) zonal stress, and (b) meridional stress. Contour interval=0.2 N 2 m - 4 38 2.11 The correlation between the wind stress predicted by the N N model and the idealized stress minus the correlation between the wind stress predicted by the L R model and the idealized wind stress, for: (a) zonal stress, and (b) meridional stress. Contour interval = 0.02. The zero contours are the dotted curves, while the negative contours are the dashed curves. Positive values means the N N predicted wind stress is outperforming the L R one. 40 2.12 The SSE of the N N wind stress minus the SSE of the L R wind stress, (both verified against the idealized wind stress), for: (a) zonal stress, and (b) meridional stress. Contour interval=0.02 N 2 m - 4 . Negative values means the N N predicted wind stress is outperforming the L R one 41 xii 2.13 A comparison of the ocean model driven by the F S U wind stress and that driven by the idealized stress (i.e. with only the first 3 EOFs), (a) SST correlation between the two models (c.i.=0.1) and (b) SST RMS error (c.i.=0.05°C); (c) HC correlation (c.i.=0.1) and (d) HC RMS error (c.i.=0.005°C) 43 2.14 A cross-validated comparison of the ocean model SST between that driven by the idealized wind stress and that driven by the empirical wind stress from the N N model (with HC as predictors), by (a) correlation (contour interval=0.05), and (b) RMS error (contour interval=0.1°C) 44 2.15 The skill differences in the model SST between the model driven by the N N model wind stress and that driven by the L R model stress, with both model SSTs verified against the standard SST, i.e. the SST from the model driven by the idealized wind stress: (a) correlation skill difference (Contour interval = 0.02), and (b) RMS error difference (contour interval=0.02°C). Positive regions in (a) indicate N N ahead of LR, while negative region in (b) indicate N N ahead of LR 46 2.16 Cross-validated skill differences in the model HC between the ocean model driven by N N model wind stress and that driven by the L R model stress, both verified against the model HC driven by the idealized wind stress: (a) correlation difference (c.i.=0.02), (b) RMS error difference (c.i.=0.05°C). Positive regions in (a) indicate N N ahead of LR, while negative regions in (b) indicate N N ahead of L R 48 2.17 Time evolution of observed SST anomalies in NIN03 (solid curve) and model HC averaged over 124°E-70°W, 5°N-5°S (dot-dashed curve) forced by the F S U observed wind stress. Both curves were normalised and smoothed by a 3-point running mean 49 xiii 3.1 A schematic diagram of coupling strategy for the coupled models 52 3.2 (a) The SST climatology of N H C M from the run of 97 years (c i . = 1°C) and (b) The observed SST climatology during 1961-1990 (c i . = 1°C) and (c) the difference, i.e. (a) minus (b). (C.i.= 0.5°C) 55 3.3 Time-longitude distribution of the mean seasonal cycle of (a) SST and (b) zonal wind stress, repeated once with the annual mean removed, for the N H C M . Negative contours are shown as dashed curves and zero contours as dotted curves. C.i.=0.2°C for (a) and 0.02 N m " 2 for (b) 56 3.4 (a) EOF1 of SST from the N H C M , (b) EOF2 of SST from the N H C M , (c) EOF1 of the observed SST and (d) EOF2 of the observed SST. Negative contours are shown as dashed curves and zero contours as dotted curves. C i . = 0.01°C for (a),(b) and (d), and c i . = 0.005°C for (c) 57 3.5 (a) EOF1 of the zonal wind stress from the N H C M , (b) EOF2 of the zonal wind stress from the N H C M , (c) EOF1 of the F S U zonal wind stress and (d) EOF2 of the FSU zonal wind stress. C i . = 0.01 N m " 2 for (a), (c) and (d), and c i . = 0.005 N m~ 2 for (b) 58 3.6 Spatial patterns of the ENSO-related POP mode for N H C M . Patterns P (real) are shown on the left, and patterns — Q (imaginary) appear on the right. SST anomalies are in panels (a) and (b), HC anomalies are in panels (c) and (d), and zonal wind stress are in panels (e) and (f). C.i.= 0.005°C for (a),(b),(c) and (d); C.i.= 0.005 N n T 2 for (e) and (f) 60 3.7 (a) The SST climatology of L H C M from the run of 97 years (c i . = 1°C), (b) The SST climatology of L H C M minus the observed SST climatology during 1961-1990 (c i . = 0.5°C) and (c) The SST climatology of N H C M minus the SST climatology of L H C M (ci . = 0.1°C) 63 xiv 3.8 The model error of subsurface ocean temperature at the depth of around 100m verified against the idealised subsurface temperature for (a) L H C M and (b) N H C M . Negative contours shown as dashed curves denote that the subsurface model temperatures are lower than the 'idealized data', c i . = 0.02°C 64 3.9 Time-longitude diagrams of ocean heat content anomalies along the equa-tor from (a) N H C M , (b) control run forced with observed F S U and (c) L H C M . C.i.= 0.1°C. Fig. 3.9b was smoothed by a 5-point running mean in time . . 66 3.10 Time-longitude diagrams of SST anomalies along the equator from (a) N H C M , (b) control run forced with observed F S U and (c) L H C M . C.i.= 1°C. Fig. 3.10b was smoothed by a 5-point running mean in time 67 3.11 Spatial patterns of the ENSO-related POP mode for L H C M . Patterns P (real) are shown on the left, and patterns — Q (imaginary) appear on the right. SST anomalies are in panels (a) and (b), HC anomalies in panels (c) and (d), and zonal wind stress in panels (e) and f. C.i.= 0.005 in all panels 68 3.12 The NINO 3 index found in N H C M and L H C M in the standard case, (a) The time series of NINO 3 index, with solid curve for N H C M and dot-dash curve for L H C M . A histogram of the number of warm events peaking in Jan. (month 1) to Dec. (month 12) for (b) N H C M , and (c) L H C M . . . . 70 xv 3.13 Analyses of the NIN03 index for N H C M and L H C M at the critical cou-pling parameter 8. N H C M are shown on the left (8 = 1.68), and L H C M appear on the right (8 = 1.5). (a) and (b) are the histograms of the number of warm events peaking in month 1 to 12 (i.e. January to Decem-ber); (c) and (d) are the power spectra; (e) and (f) are the phase planes reconstructed from the time series 73 4.1 Correlation skills for the predictions of SST anomalies in the NIN03 area as a function of the forecast lead time for (a) the hindcast period from 1964-1990, and (b) the forecast period from 1990-1998, by N H C M (solid), L H C M (dashed) and persistence (dot-dashed). Both HCMs use the 'stress-only' initialization scheme 77 4.2 Predicted NIN03 SST anomalies (°C) by the N H C M from 1964-1990 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for observations . . . 79 4.3 Predicted NIN03 SST anomalies (°C) by the N H C M from 1990-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for observations 80 4.4 Predicted NIN03 SST anomalies (°C) by the L H C M from 1990-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for observations 81 4.5 Predictive (correlation) skills of the NIN03 SSTA from (a) L H C M and (b) N H C M , for the 1990s (solid), 1980s (circled), 1970s (dashed), and the 1960s (dot-dashed), as a function of the lead time. Both coupled models use the 'stress-only' initialization scheme 83 xvi 4.6 Nonlinear C C A (NLCCA) mode 1 for the 1980s shown in the P C (prin-cipal component) space of the zonal wind stress anomalies. The N L C C A analyzed the relation between the first three PCs of the zonal wind stress anomalies and the first three PCs of the SSTA. The data are shown as dots, and their projection onto the N L C C A mode 1 are the circles, and the classical C C A solution is shown as a straight line for comparison. The solutions are shown in (a) the PC1-PC2 space, (b) the PC1-PC3 space, (c) the PC2-PC3 space, and (d) a 3-D view. Following the principle of par-simony, the architecture of the N L C C A was set for minimal nonlinearity, i.e. only 2 hidden neurons used, and the weight penalty parameter was 0.1 in all cost functions (Hsieh 2001) 85 4.7 N L C C A mode 1 for the 1980s in the P C space of the SSTA 86 4.8 N L C C A mode 1 for the 1990s in the P C space of the zonal wind stress. . 87 4.9 N L C C A mode 1 for the 1990s in the P C space of the SSTA 88 4.10 Seasonal dependence of the predictive skills for the NIN03 SSTA dur-ing 1964-1998 from (a) L H C M and (b) N H C M . Both coupled models use the stress-only initialization scheme initialized on January 1 (dot-dashed), April 1 (solid), July 1 (dotted) and October 1 (dashed) 89 4.11 Seasonal variation of variances of NINO3 index for (a) observed SST, (b) observed HC , (c) observed zonal (solid) and meridional (dashed) wind stress and (d) model SST taken from control run forcing with F S U winds. 90 4.12 Contrast ofthe predictive skills for the NIN03 SSTA during 1964-1990 be-tween the stress-only initialization scheme (solid line) and the stress-from-ocean initialization scheme (dashed line) for (a) L H C M and (b) N H C M . 93 xvii 4.13 Contrast of the predictive skills for the NIN03 SSTA during 1990-1998 between the stress-only initialization scheme (solid line) and the stress-from-ocean initiahzation scheme (dashed line) for (a) L H C M and (b) N H C M . 94 4.14 Predicted NIN03 SST anomalies (°C) by N H C M using the stress-from-ocean initiahzation scheme during 1964-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Sohd line is for the predicted values and dash line for the observations 95 4.15 Predicted NIN03 SST anomalies (°C) by L H C M using the stress-from-ocean initiahzation scheme during 1964-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Sohd line is for the predicted values and dash line for the observations 96 4.16 Correlation skills of the predicted NIN03 SSTA for (a) the hindcasts period from 1964-1990, and (b) the forecast period from 1990-1998, by N H C M (sohd), L H C M (dashed). Both H C M models used the stress-from-ocean initialization scheme 97 4.17 Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1982. The observed SSTA from January 1982 to March 1983 is shown in the middle panel. The ordinate ranging from 1 to 15 is the calendar month from January 1982 to March 1983 99 4.18 Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1997. The evolution of the observed SSTA from January 1997 to March 1998 is shown in the middle panel 100 x v i n 4.19 Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1998. The evolution of observed SSTA from January 1998 to March 1999 is shown in the middle panel 101 5.1 The time series of (a) the first principal component (PC) (i.e. E O F time coefficient) and (b) the second P C , for the observed SSHA (solid line), the SSHA (dotted line) from the control run (with the ocean model forced by the F S U wind stress), and the first interface HI A (dash fine) from the control run 109 5.2 The spatial patterns of the first two E O F modes from the control model run for the SSHA (a and b) and for the HI A (c and d). Solid fines indicate positive contours, dashed fines, negative contours, and dotted lines, the zero contour 110 5.3 Scatter plot of the first P C of the H C A versus (a) the second P C of the SSTA, and (b) the first P C of the H2A, the thickness anomalies of the second layer. The observations are shown as dots, the N N relation is indicated by the small circles and the linear regression by the straight line. 112 5.4 Time-longitude plots of the ocean heat content anomalies along the equa-tor during 1980-1998 from (a) the observations and (b) the control run. (c) The H C A calculated from the prognostic variables estimated by the cross-validated ensemble N N model using (b) as predictors. Values prior to 1980 are not shown, as the quality of the observed HC is poor 113 xix 5.5 Correlations between the observed tropical Pacific SSTA and the SSTA from the ocean model forced by the F S U wind during 1980-1989, (a) with-out data assimilation, (b) with SSTA assimilation, (c) with SSHA assimi-lation, and (d) with H C A assimilation 116 5.6 Correlations between the observed tropical Pacific SSTA and the SSTA from the ocean model forced by the F S U wind during 1990-1998, (a) with-out data assimilation, (b) with SST assimilation, (c) with SSHA assimila-tion, and (d) with H C A assimilation 117 5.7 Time evolution of observed and modeled Niiio3 SSTA during January 1980-December 1997. Solid line denotes observed SSTA in all figures and dash line is from the ocean model run (a) without assimilation, (b) with SST assimilation, (c) with SSHA assimilation, and (d) with H C A assimi-lation 119 5.8 Correlation between observed and predicted SST anomalies in the Niho3 region, as a function of lead time for (a) the 1980s and (b) the 1990s. Solid line is from the prediction without data assimilation. The correlation from the predictions with data assimilations are shown by the dash line for the wind stress assimilation, asterisks for SST assimilation, dot-dash for SSHA assimilation and circles for H C A assimilation 121 5.9 (a) Correlation between the monthly mean heat content anomalies in the Niho3 and Nirio4 regions as a function of time lag between the two time series, with H C A from (a) the model with H C A assimilation and (b) ob-servations. Dashed line is for the 1980s and solid line for the 1990s. A positive time lag means that the signal in Nifio4 is leading the signal in Niiio3 123 xx 5.10 Time-longitude diagrams of predicted SSTA along the equator for several assimilation experiments for the 1982/83 warm event. The predictions are all initialized on January 1, 1982. For comparison, the observed SSTA and the SSTA from the model without data assimilation are shown. The ordinate with values 1 to 15 is the calendar month from January 1982 to March 1983 125 5.11 Same as in Fig. 5.10 but for the 1988/89 cool event. The prediction was initialized on January 1, 1988 127 5.12 Same as in Fig. 5.10 but for the 1997 warm event. The prediction was initialized on January 1, 1997 128 6.1 The Lorenz system integrated over 15000 time steps for the weakly non-linear case (left column) and the highly nonlinear case (right column). . . 133 6.2 The ensemble mean of 25 NNs (dot-dash curve) for approximating the third Lorenz equation (solid curve) for the weakly nonlinear case (left column) and the highly nonlinear case (right column). Top panels (a) and (b) are for the training period of 3000 time steps, and bottom panels (c) and (d) for test period of 1000 time steps 135 6.3 The hybrid model (dot-dash curve) and the true Lorenz model (solid curve) integrated from identical initial conditions for the weakly nonlinear case (left column) and the highly nonlinear case (right column) 136 6.4 Schematic diagram for the (a) SC (Strong Continuity constraint), (b) N C (No Continuity constraint) and (c) PSC (Partial Strong Continuity con-straint) assimilation schemes. The solid curves are the model trajectory, dots, the observed values and dash lines, the model errors to be minimized. Here a model trajectory starts from an observed value 140 xxi 6.5 The NC assimilation results, with the model training results shown on the left panels, and predictions over the test period on the right panels; where the solid curves show the true data and dot-dash curves the outputs of the hybrid model. The dot-dash curves overlap the solid curves in the left panels 141 6.6 The PSC assimilation results with AS= 200 time steps. The model train-ing results are shown on the left, and model predictions on the right, with solid curves for the true data and dot-dash curves for the outputs of the hybrid model 143 6.7 Same as Fig. 6.6 but for AS= 1000 time steps 144 6.8 Correlation skills (averaged over 10 experiments) for the variables X, Y and Z at various assimilation windows (T = 100, 300, 500 and 800 time steps) during (a) the training period, and (b) the test period 146 6.9 (a) Correlation skills and (b) R E E for the variables X, Y and Z during the training and test periods (averaged over 100 experiments), for the weakly nonlinear case, where the SC assimilation retrieves the N N parameters. The error bars show the standard deviation. Results from the simple hybrid model without data assimilation (Section 3) are also shown for comparison 149 6.10 Bar charts from 100 assimilation experiments for the highly nonlinear case showing the number of experiments with (a) correlation above a particular correlation level and (b) R E E below a particular R E E level, during the test period 151 xxi i 6.11 A typical experiment with class 1 results for the highly nonlinear case, showing (a) the fitting during the training period (left panels) and (b) the forecasting during the testing period (right panels). Sohd curves denote the true value, while dot- dashed curves denote model results 152 6.12 Same as Fig. 6.11 but for a typical experiment with class 2 (i.e. moderately successful forecast) results 153 6.13 Same as Fig. 6.11, but for a typical experiment with class 3 (i.e. unsuc-cessful forecast) results 154 6.14 The system trajectories in the X — Y — Z space during the training pe-riod (thin curves) and the testing period (thick curves) for (a) the class 1 experiment and (b) the class 3 experiment 155 6.15 The results from retrieving the dynamical parameters a, b and the initial conditions of the hybrid model, where the initial guesses were the true values scaled down by 20% 157 6.16 The correlation skill averaged over 100 experiments where the N N param-eters, the dynamical parameters a and 6, and the initial conditions were retrieved for the weakly nonlinear case. Error bars indicate the standard deviation 158 6.17 Bar charts from 100 assimilation experiments for the highly nonlinear case showing the number of experiments with (a) correlation above a particular correlation level and (b) R E E below a particular R E E level, during the test period. The N N parameters, the dynamical parameters a and b, and the initial conditions were retrieved 160 xxin .1 V a r i a t i o n s of the n o r m a l i z e d cost func t ion ( J / J 0 ) (sol id l ine) a n d n o r m a l -i z e d gradient ( | |g | | | | g o | | _ 1 ) (dash l ine) w i t h the number o f i te ra t ions for (a) the w e a k l y nonl inear case of F i g . 6.7 and (b) the s t rongly nonl inear case of F i g . 6.11. T h e bet ter i n i t i a l guesses i n (a) meant tha t the cost f unc t i on a n d its gradients had re la t ive ly less dis tance to drop to convergence t h a n i n (b) XXIV List of Acronyms AS Assimilation Segments EOF Empirical Orthogonal Function ENSO E l Nino/Southern Oscillation 4D-Var Four-Dimensional Variational Data Assimilation G C M General Circulation Model HC Heat Contents ITCZ Inter-Tropical Convergence Zone LR Linear Regression L H C M Linear Hybrid Coupled Model N H C M Nonlinear Hybrid Coupled Model NN Neural Network NC No Continuity Constraint Assimilation RMSE Root Mean Square Error PCA Principal Component Analysis PSC Partially Strong Continuity Constraint Assimilation SC Strong Continuity Constraint Assimilation xxv SSE Sum of Square Error SSH Sea Surface Height SSHA Sea Surface Height anomalies SST Sea Surface Temperature SSTA Sea Surface Temperature anomalies SVD Singular Value Decomposition. TAO Tropical Atmosphere Ocean 3D-Var Three-Dimensional Variational Data Assimilation XBT Expendable bathythermograph xxvi Preface This thesis presents my efforts to improve ENSO (El Nino-Southern Oscillation) sim-ulation and prediction using hybrid coupled models (HCMs). The goals of this thesis are to investigate two essential problems in hybrid coupled modelling— (i) atmospheric model construction and (ii) oceanic model initialization. Towards the first goal, a non 1 linear technique, neural networks, is used to construct a nonlinear empirical atmosphere instead of the traditional linear atmosphere in the H C M . For the second goal, the 3D-Var data assimilation technique is applied to the H C M . While the full 4-D var H C M is beyond the scope of this thesis, a neural-dynamical hybrid approach under 4-D var has been developed to study the simple Lorenz system. Because of the weakly nonlinear nature of the tropical ocean-atmosphere system at ENSO time scales, the difference between the nonlinear atmosphere H C M and the lin-ear atmosphere H C M is not dramatic. But this study provides an extension and an alternative viewpoint to the hybrid modelling strategy— supplying information on the linearity and nonlinearity of the tropical Pacific ocean-atmosphere interaction, which has implications beyond the focus of this study. The goals of this thesis were arrived at gradually during the four years of research work. Only towards the end, and after five journal papers have been written (Tang et al., 2001; Tang and Hsieh, 2001a; Tang, 2001; Tang and Hsieh, 2001b and Tang and Hsieh, 2001c), did the final structure of the thesis crystallize to a clear shape. Tang et al (2001) is presented in Chapter 2, introducing the ocean model and developing the two atmosphere models. Tang (2001) and Tang and Hsieh (2001b) cover Chapter 3 and Chapter 4, respectively, comparing the two HCMs in terms of ENSO simulation (Chapter 3) and xxvii ENSO prediction (Chapter 4). The 3-D var and 4-D var data assimilation, presented Chapter 5 and Chapter 6, respectively, were from Tang and Hsieh (2001c, 2001a). xxvm Acknowledgements I would like to thank my supervisor Dr. William Hsieh for his help and guidance during the course of this work. He has provided such a creative academic atmosphere that I have freely developed my interest and built up my knowledge in my studies. I am sincerely thankful for the constructive suggestions and comments from the mem-bers of my committee— Dr. Susan AUen and Dr. Lionel Pandolfo. As well, I would like to acknowledge Dr. Richard Pawlowicz and Dr. Phil Austin for their help in the under-standing of ocean dynamics and tropical atmospheric process. None of this work could have been accomplished without the computer support pro-vided by Denis Laplante. I am grateful to Dr. Keith Haines at Edinburgh University, Dr. Magdalena Alonso Balmaseda and Dr. David Anderson at E C M W F for kindly providing us the ocean model. In particular, Dr. Haines and Dr. Balmaseda gave the author insightful recommendations and advice about the ocean model. Special thanks go to my friends, Dr. Adam Monahan at Humboldt-Universitaet zu Berlin and Dr. X u Li at Maryland University, for many fruitful and simulating discussions which greatly benefit the work. Dr. Adam Monahan also spent a lot time to read the earliest version of this thesis and gave many detailed comments. Also my colleagues, past and present, improved this work with useful comments and suggestions. Thank you, Jingxi Lu, Benyang Tang, Yuval, Aiming Wu. Most of all, I would like to express my love and gratitude to my family. Without their understanding and support, this would be impossible for me to finish the thesis in 4 years. xxix Chapter 1 Introduction 1.1 E N S O theory 1.1.1 ENSO and its characteristics ENSO is an acronym for the E l Nino-Southern Oscillations, itself a combination of two terms. Originally, E l Nino referred to an annual weak warm current that flows south-ward along the coast of southern Ecuador and northern Peru about Christmas time, but subsequently became associated with the unusually large warming that occur at the equatorial eastern Pacific and off the coast of Peru every few years, changing the local and regional ecology. Southern Oscillation (SO) referred to the interannual seesaw pres-sure fluctuations between the Indian Ocean and the eastern tropical Pacific— "When pressure is high in the Pacific Ocean, it tends to be low in the Indian Ocean from Africa to Australia" (Philander 1990). It became clear by the beginning of the 1980s that the E l Nino and SO are intimately related, with the large-scale warming (El Nino) coordinated with a negative SO phase, i.e., lower pressure in the eastern-central tropical Pacific and higher pressure in the western tropical Pacific-Indian Ocean compared with normal conditions. The seesaw structure over the tropic Pacific is associated with the relaxation of the trade winds over the equatorial Pacific and a southward displacement ofthe ITCZ (Intertropical Convergence Zone). Hence scientists coined the term ENSO to describe this large-scale, coupled atmosphere-ocean interaction. The robust tropical features of the observed warm ENSO 1 Chapter 1. Introduction 2 event are as below (Battisti and Sarachik 1996): (i) quasi-stationary warm SST anomalies appear in the eastern and central equatorial Pacific; (ii) at an event onset, a relaxation of the trade winds is associated with positive SST anomalies; (iii) the thermocline deepens in the east and shoals in the west along the equator; (iv) the oscillation of ENSO is irregular, with a quasi-period of 3-7 years; (v) the peak anomaly usually occurs in winter (i.e. phase-locking to the annual cycle). In the next subsections, we will briefly summarize several of the most influential hypotheses about the ENSO mechanism. 1.1.2 Bjerknes hypothesis Although Walker documented the interannual variations in the circulation of the trop-ical and global atmosphere and proposed the famous concept of Southern Oscillation at the turn of this century, the study of ENSO in view of the relation between inter-annual oceanographic and meteorological variations in the tropical Pacific did not start until Bjerknes proposed the well-known "trade-wind" theory. In his famous paper "At-mospheric teleconnections from the equatorial Pacific" (Bjerknes 1969), he organized the intertwined elements of the Southern Oscillation, E l Nino, and large-scale air-sea interactions into a new conceptual framework supported by plausible dynamic and ther-modynamic reasoning. Under the normal action of the easterly trade winds, a system of westward surface currents accumulate in the upper waters of the western equatorial Pacific, with a deep thermocline (150-200m), in contrast to the shallow thermocline in the eastern boundary (30-50m). According to the Ekman theory, the easterly winds cause a steady movement of the water away from the equator (Ekman transport) into both hemispheres, creating a horizontal divergence of surface water that in turn induces an equatorial upwelling of water from below. Upwelling produces an elongated east-west tongue of cool surface Chapter 1. Introduction 3 water in the central and eastern equatorial ocean, where the winds are strong and cooler water is brought up from beneath the shallow thermocline. The cooling does not occur in the western Pacific where the winds are weaker and the thermocline is deep (150-200m). When the easterly trade winds decrease, the forces that maintain the zonally asym-metric thermocline and the upwelled cold tongue disappear, so the thermocline in the east becomes depressed and equatorial upwelling weakens. Meanwhile, the upper layer of warm water moves equatorward and eastward, resulting in a readjustment, or "relax-ation", of the ocean system away from the normal state. A l l of these local and remote responses to the wind contribute to the rises in sea level and sea surface temperature (SST) at the eastern equatorial Pacific that are characteristic of E l Nino. Because of the deep thermocline in the western Pacific, the readjustment of ocean system there is not as pronounced as in the east. A similar argument but with opposite conditions and consequences can be made to explain the cold phase of ENSO, i.e., the La Nina (also called the E l Viejo). In this paradigm, ENSO arises as a self-sustained cycle in which anomalies of SST in the Pacific cause the trade winds to strengthen or slacken and that in turn drives the ocean circulation changes that produce anomalous SST. For example, when the easterly trade wind weakens or withdraws, the warming in the east reduces the zonal SST gradient and hence weakens the Walker circulation, resulting in weaker trade winds at the equator. Such a positive feedback between SST and trade winds leads to unstable interactions between the ocean and the atmosphere in the equatorial Pacific, so anomalous episodes of ENSO develop. However, Bjerknes did not include the mechanisms for moving the coupled system from a phase with warm SST anomalies to one with cold SST anomalies, nor the mechanism for stopping the warming (the positive feedback) in the eastern Pacific. There was also uncertainty on where the initial changes in SST (or in the trade wind) come from in order to start the unstable growth. These questions had puzzled scientists Chapter 1. Introduction 4 until the delayed-oscillator theory was proposed at the late 1980s. 1.1.3 Delayed-Oscil lator mechanism The equator is an efficient waveguide for atmospheric and oceanic waves trapped to propagate along the equator due to the vanishing Coriolis force at the equator. These waves play an important role in the adjustment of the equatorial ocean and tropical atmosphere in response to changes in the surface forcing. Two types of waves, the equatorial Kelvin and Rossby waves, are of particular importance for the adjustment of the ocean on seasonal to interannual time scales. The Kelvin wave is a special gravity wave, which propagates eastward with a speed of approximately 2-3 m s _ 1 and has its maximum amplitude centred at the equator. The Rossby waves, whose propagation is prompted by the latitudinal variation of the Coriolis parameter (the /3-effect), propagate westward at a rate of 0.6-0.8 m s - 1 near the equator. The oceanic Kelvin and Rossby waves propagate energy and momentum received from the wind stress and they provide the oceanic memory that is so important to interannual variability and ENSO. Similar waves exist in the atmosphere, but their propagation rates are far greater than the oceanic counterparts. Therefore, the adjustment time scale of the tropical atmosphere to changes in SST is much shorter (10 days or less) than the adjustment time scale of the equatorial ocean (approximately six months) to changes in the wind stress. The short adjustment time of the atmosphere allows one to make the assumption (good to zeroth order) that the atmosphere is in a statistical equilibrium with the SST on time scales longer than a few months. Thus, the memory of the state of the climate system primarily resides in the ocean. The delayed-oscillator mechanism is based on the theory of equatorial waves and was Chapter 1. Introduction 5 proposed by Schopf and Suarez (1988) and Battisti and Hirst (1989). In the delayed-oscillator theory, the ENSO cycle ascribes a central role to equatorial oceanic wave pro-cesses which affect sea surface temperature (SST) and subsequent ocean-atmosphere in-teractions through redistribution of upper ocean heat content. According to the delayed action oscillator scenario, the easterly wind anomalies over the western Pacific prevailing during the cold phase of the ENSO cycle force an upwelhng Kelvin wave packet, which propagates eastward along the equator and causes coohng at the sea surface in the east-ern Pacific, where the thermocline is shallow. The ocean response to the easterly winds in the west, however, consists also of a downwelling Rossby wave packet, which propa-gates westward. This Rossby waves response has its strongest signals off the equator. It does not influence the SST in these regions because the thermocline is deep in the western Pacific. The Rossby waves then reflect at the western boundary into a packet of downwelling Kelvin waves, which have maximum amplitudes at the equator and propa-gate eastward. Once it has propagated far enough into the eastern Pacific, it is able to affect the SST, and a positive SST anomaly develops, which can grow by unstable air-sea interactions into the warm phase of ENSO. Thereafter the sequence of events repeats itself but with reversed signs. In this view, ENSO is a low-frequency basin-wide mode of oscillation which is perfectly predictable. 1.1.4 E N S O chaos and stochastic forcing While the delayed-oscillator answers the questions on how the ENSO cycle moves from a warm phase to a cold phase and on what stops the warming, it still fails to explain the irregularity and the phase-locking in the ENSO cycle. In the early 1990s, some studies with theoretical models and intermediate coupled models showed that the chaotic behaviour inherent to nonhnear interactions between annual cycle and coupled ENSO mode are responsible for the irregularity and phase locking. For example, Jin et al. Chapter 1. Introduction 6 (1994) and Tziperman et al. (1995) found independently that nonlinear interaction of the annual cycle and the coupled ENSO mode leads to ENSO chaos. Chang et al. (1995, 1996) investigated the relationship between such a frequency-locking behaviour and the transition to chaos in several ENSO coupled models. They found that a further increase of the coupling strength could lead to chaotic behaviour of the coupled ENSO mode. However, there has been considerable debate as to whether chaos is a main cause of ENSO irregularity and phase-locking because some coupled models, e.g., the coupled general circulation model (Philander et al. 1992), lack low-order chaotic behaviour. The role of noise forcing on ENSO cycle has been addressed at various times during T O G A (Tropical Ocean-Global Atmosphere), especially since the late 1990s. Recent studies show that the effects of realistic noise apphed to a hybrid coupled model or an intermediate coupled model in a regime that would otherwise be periodic are sufficient to produce irregularity generally consistent with observed ENSO signals. (Blanke et al. 1997; Eckert and Latif 1997). The spectral signature in the presence of noise typically differs from models with irregularity induced purely by chaos, which tend to retain fairly sharp dominant spectral peaks. 1.2 E N S O models and E N S O prediction Over the last decade, a considerable number of models have been developed to investigate the ENSO mechanism and to improve ENSO prediction (Neelin et al. 1998; Latif et al. 1998). The operational ENSO predictions several seasons in advance are performed at present by about a dozen of forecast models. The skills attained by the models are generally comparable, with around 0.6 correlation at a lead time of 6 months. Models for ENSO prediction can be categorized into purely statistical models, hybrid coupled models and coupled ocean-atmosphere models. Advanced statistical techniques Chapter 1. Introduction 7 have been widely used, ranging from linear models (Barnett and Preisendorfer 1987) to nonlinear models (Tangang et al. 1998a, b; Hsieh and Tang 1998). Statistical models have several drawbacks as discussed in Latif et al. (1998): (1) forecast skills degrade relatively quickly with increasing lead time; (2) artificial skill remains a problem in all statistical forecast schemes; and (3) a statistical relationship does not always imply a causal relationship. In coupled ocean-atmosphere models, the forecast skills drop more slowly with in-creasing lead time than the statistical models. The problems associated with the coupled models are: (1) climate drift exists in almost all CGCMs (coupled general circulation models); (2) some physical processes still have to be parameterized; (3) proper initial-ization of the coupled model is difficult; and (4) the cost is much greater than statistical models. The hybrid-coupled model (HCM) combines a dynamical ocean model with a statis-tical atmospheric model for coupled interaction studies (Barnett et al. 1993; Syu and Neelin 2000; Balmaseda et al. 1994). H C M shares the advantages of both dynamical models and statistical models, although some drawbacks associated with the dynamical model and statistical model are also retained. The main advantages of a hybrid model are: (1) climate drift problem avoided, (2) easier understanding of the coupled mecha-nisms and lower computing cost than a full C G C M (Blank et al 1997), (3) comparable or even better simulation and prediction skills relative to a C G C M (Palmer and Anderson, 1994; Kang and Kug 2000). Currently, all empirical atmospheric models used in HCMs have been linear statistical models. Although the dynamic response of the atmosphere and ocean in the tropical Pacific is quasi-linear (Tang et al. 2000; Tang et al. 2001), the coupled interactions between the two components seem more complicated than the relations extracted from data and from uncoupled model responses. In this thesis, two HCMs of the tropical Chapter 1. Introduction 8 Pacific were developed and studied for their interannual variability. The first, the 'linear' hybrid coupled model (LHCM) has a 6-layer nonlinear dynamical ocean model coupled to a linear regression atmospheric model. The second, the 'nonlinear' H C M (NHCM) has the same ocean model but a nonlinear neural network atmospheric model. The differences between the N H C M and L H C M in ENSO prediction are also studied. 1.3 D a t a assimilation in E N S O prediction An important issue affecting the ENSO prediction skills is the initialisation scheme of the prediction models. When Zebiak and Cane (1987) first made ENSO predictions with dy-namical coupled model, they used the observed surface wind stress in the tropical Pacific to force the ocean model for several years prior to the prediction data. Predictions were then made by running the coupled model forward in time from the established ocean initial conditions. Due to the fast response of the atmosphere to the ocean, the initial conditions of the atmosphere model was a 'slave' to the ocean state as expressed through sea surface temperature (SST). This initialization scheme has been widely adopted in subsequent intermediate coupled models and coupled general circulation models (Bal-maseda et al. 1994; Latif et al. 1993; Kleeman 1993; Kirtman and Zebiak 1997) and can produce a reasonable prediction skill of around 0.6 correlation at a lead time of 6 months. The above initialization scheme only uses the surface wind stress information and ignores the model feedback from the ocean to the atmosphere. Although it has shown a great deal of success in establishing the initial conditions for prediction models, it has also been of interest to initialise the prediction models considering the feedback from the ocean to the atmosphere and using other observations, in particular the direct observa-tions of ocean states. The hybrid coupled model of Barnett et al. (1993) is initialised with observed SSTs and obtained similar hindcast skills. Chen et al. (1995) developed Chapter 1. Introduction 9 an ocean initialization scheme for the Zebiak and Cane (1987) coupled model with the consideration of feedback from ocean to atmosphere through a nudging approach. Fur-ther, they assimilated the sea level data using a similar nudge scheme for the coupled model (Chen et al. 1998) improving the forecast skills of the Lamont model. Kleeman et al. (1995) assimilated the 400-m depth-averaged equatorial temperature anomaly into an intermediate complexity coupled model by using the adjoint variational technique. In their work, the subsurface information was used as proxy data and an empirical relation between it and the model prognostic variables were used to initialise the model states. The results indicated that such an initialization approach can lead to a marked increase of hindcast skills, in particular for long lead times of over 8 months. A similar approach using the proxy data to assimilate subsurface structures were also performed by Fischer et al. (1997), in which the sea level observations are inserted into the 3-dimensional tem-perature field with a statistically based relation. J i et al. (1998) and Rosati et al. (1997), on the other hand, used the observed three-dimensional temperatures to directly initialize the coupled general circulation models through a 3D Var assimilation system, and found that initialization of the ocean by the assimilation of observed subsurface temperature data generally leads to improved prediction skills. However, as mentioned in Latif et al. (1998), compared with the field of numerical weather prediction, techniques of data assimilation and initialization in ENSO prediction are still at a relatively early stage of development. We are at present unclear on a basic issue, i.e. which observations are most useful in initialising ENSO predictions and how do the prediction skills differ for assimilating different types of data? It seems that we need a systematic comparison of assimilating different types of data on ENSO predictability based on an identical assimilation system and an identical forecast model. Since ENSO is a fundamental oscillatory mode of the coupled ocean-atmosphere system and that the memory of the coupled system resides in the ocean thermocline state (Neelin et al. 1998), Chapter 1. Introduction 10 one might expect that such a systematic comparison can provide insights for most ENSO models which share the same basic ENSO dynamics. However, at the present stage, the most useful data to be assimilated might still be model dependent—this is not totally satisfactory but we are still at a relatively early stage in understanding ENSO and its predictability (Latif et al. 1998). 1.4 Nonlinear physical processes in the tropical atmosphere Nonlinear physical processes that will be discussed here are only confined to the atmo-spheric response (coupled) processes to the ocean (SST), as at low frequencies the atmo-sphere over the tropical Pacific can be regarded as a strongly forced quasi-equilibrium system governed by the state of the tropical Pacific SST (Latif et al, 1998). Some other physical processes, such as the feedback of the extra-tropical atmosphere to the tropical atmosphere, the internal variability in the tropical atmosphere etc., are neglected because they make small contributions to the ENSO oscillations. 1.4.1 Physical processes associated with SST in the tropical atmosphere Observations indicate latent heat release to be the primary energy source in the tropics. The latent heat release occurs in association with convective cloud systems and verti-cal motion from the external forcing specified by sea surface temperatures (SST). The scale analysis of large-scale tropical motions show that the diabatic heating is mainly nonlinearly balanced by vertical motion (Holton 1992; Hoskins and Pearce 1983). A comparison of the SST pattern with the convective patterns reveal that, over the ocean, moist air rises where SSTs are highest while dry air subsides where surface waters are cool (Philander, 1990). The warmest waters are in the western tropical Pacific and in a band of latitudes just north of the equator where the ITCZ (Intertropical Convergence Chapter 1. Introduction 11 Zone) is located. The large-scale convection over the warmest surface waters is associated not with increased local evaporation from the ocean but with the convergence of moist air onto those regions. High SST is a necessary condition for active deep convection, the intensity of which increases with SST in the range from 26°C to 30°C. However, convec-tion can have a large spatial variability which SST does not have and its occurrence is intermittent, indicating that relation between SST and deep convection is very complex (Delecluse et al., 1998). Changes in the strength and position of the convection over the warm pool can affect the sea level pressure. Variations in the pressure gradients in turn are correlated with variations in the trade winds over the Pacific ocean by the Hadley circulation and Walker circulation variations. During the warm phase of ENSO (El Nino), the ITCZ moves equatorward, the South Pacific Convergence Zone moves northward, and the convergence zone over the western Pacific moves eastward. The zonal Walker circulation is weak; the trade wind weakens and retreats eastward during these periods. The opposite happens during the cold phase of ENSO ( La Nina). Therefore, the tropical SST changes influencing the atmosphere is through the surface fluxes of moisture, heat, and momentum and a readjustment of the tropical circulation in a thermally direct sense. These fluxes depend directly on the SST, but they also depend nonlinearly on the wind speed, the near-surface air temperature and the near-surface specific humidity. The response of the tropical atmosphere to SST is associated with important processes such as the variations in the ITCZ, the Walker cell and the Hadley cell, in position and extent. These physical processes should yield highly nonlinear behaviour. In other words, the processes involved can be complex and may involve feedback mechanisms which make the processes anything but straightforward (Zebiak, 1985). Hence it is important to be able to model the nonlinear response of atmospheric winds to SST. Chapter 1. Introduction 12 1.4.2 Nonlinearities and linearization For a simple view on how the tropical SST anomalies affect the atmosphere, the linear theory for small perturbations of a resting atmosphere, i.e. with the heating rate assumed to be small enough for linear theory to apply, was first developed by Gil l (1980). Using the linear shallow water equations, Gill investigated the response of the tropical atmo-sphere to a given distribution of heating and obtained some realistic tropical atmospheric characteristics. In the field of ENSO studies, most intermediate atmospheric models fol-low the essence of the Gill model (Zebiak 1986, Zebiak and Cane 1987; Battisti 1988). A typical Gil l model can be written as (Neelin et al, 1998): £aua - fva + dxp = 0 (1.1) eava + fua + dyp = 0 (1.2) eap + c2a(dxua + dyva) = -Q (1.3) where ea is an inverse damping time, ua and va are low-level velocities in the atmosphere, and p is a pressure perturbation scaled by the mean density, which is sometimes taken proportional to a layer mean temperature perturbation Ta. The forcing, Q, and the effective phase speed, c a , have different interpretations in different approaches. The success of these simple atmospheric models is partly because, to a first approx-imation, what is required of them in ENSO theory is modest, a rough representation of the anomalous wind stress response to the SST anomaly. Because the ocean response to the wind stress tends to be non-local, it turns out that the coupled response can, in some circumstances, be fairly forgiving of atmosphere errors in position or extent of wind stress anomaly (Neelin et al, 1998). Nonetheless, it was a piece of good fortune for ENSO studies that the complexity of tropical atmospheric dynamics could be even roughly captured in simple models. For mid-latitudes, the current lack of such models is Chapter 1. Introduction 13 a primary hurdle to advancements in mid-latitude coupled theory. However, even for the Gill models, the response of tropical atmosphere to the heat forcing is still nonlinear by nonlinearly parameterising the forcing Q (Gill, 1980; Zebiak, 1986; Battisti 1988). When the forcing Q is specified from SST, the problem of atmo-spheric response becomes harder because the known relation between SST and convection is far from perfect. The fact that the coupled models with a linear atmosphere can realistically model ENSO behaviour (e.g., Barnett et al. 1993; Balmaseda et al. 1994) indicates that the nonlinearity in the tropical atmosphere at very low frequencies is not strong. Fig. 1.1 shows the leading three E O F time series for the zonal wind stress and the SST, indicating that the first E O F time series of the zonal wind stress has a good linear relation to the first E O F time series of the SST, but the linear relations dramatically weaken for the higher E O F modes. The same happens for the relations between the meridional component of the wind stress and SST (not shown). Our results (see Chapter 2) show that the skill improvements are very small when we predict the first E O F mode of a wind stress component from the upper ocean state using a nonlinear neural network atmosphere instead of a linear regression atmosphere. But there is improvement using a neural network model when predicting the second and third E O F modes of a wind stress component. 1.5 Objectives of the thesis In the thesis, the main objective is to develop a better hybrid coupled model of the tropical Pacific, and to improve ENSO prediction. One important aspect affecting H C M performances is the construction of the atmo-spheric model, i.e. the method for estimating the surface wind stress for a given upper Chapter 1. Introduction 14 Figure 1.1: E O F time series for SST and zonal wind stresses. Solid line is for wind stresses and dashed line for SST. Upper: E 0 F 1 ; middle: E0F2 ; bottom: E 0 F 3 . A l l time series are smoothed by a 5-month running mean filter. Chapter 1. Introduction 15 ocean state and the choice of representative ocean states as predictors. Previous empiri-cal atmospheric models ranged from the hnear correlation model (Latif et at 1990, Latif and Villwock 1991) to the SVD model (Syu et al. 1995). However, all the atmosphere models used for H C M so far are based on hnear statistical theory. This motivates the following questions: How big are the errors from assuming a hnear wind stress response to the upper ocean changes over the seasonal time scale? How much improvement will occur if nonlinear theory is applied to construct the empirical atmospheric model in a H C M ? Another important aspect affecting model predictability is the initialization scheme of the numerical model. How to optimally choose the initial conditions, i.e., how to assimilate data, especially subsurface heat content data into an intermediate complexity layered ocean model with low vertical resolution is the another objective of the thesis. The effects of assimilating different types of data on model predictability is also investigated. The last objective is to try solving a potential problem inherent to the traditional hybrid model. In the traditional hybrid model, the atmospheric model is constructed from the observed data or the output of the ocean model. Such a constructed statisti-cal equation cannot be guaranteed to be consistent with the oceanic model dynamics, i.e. potential inconsistency exists between the constructed statistical equation and the other dynamical equations. This might lead to dramatic errors after a period of model integration. Therefore, in this thesis, a new approach is proposed which combines a neural network and four dimensional variational data assimilation to simulate a missing dynamical equation. This thesis is structured as follows: Chapter 2 describes the construction of two atmospheric models and discuss the response of the ocean to the two atmospheric models. Chapter 3 displays the ENSO characteristics in the two hybrid coupled models. The performances of the two coupled models in terms of ENSO prediction are examined in Chapter 1. Introduction 16 Chapter 4. The issues as to how to assimilate different types of data into our coupled model and how the assimilations affect ENSO simulation and prediction are investigated in Chapter 5. Chapter 6 formulates a general procedure for optimally constructing a missing dynamical equation using N N and variational data assimilation, and presents some numerical experiments. The summary and discussions are given in Chapter 7. Chapter 2 The construction of two hybrid models 2.1 Introduction The term hybrid coupled model (HCM) was introduced by Neelin (1990) to describe a coupled model consisting of a simple atmospheric model and a complex oceanic model. Latif and Villwock were motivated by this concept to couple an empirical (statistical) atmosphere model to an oceanic G C M (Latif and Villwock 1990). Since then, the term hybrid coupled models usually refers to models which couple a statistical atmospheric model to a dynamical ocean model (e.g. Latif and Villwock 1990; Syu et al. 1995; Barnett et al. 1993; Balmaseda et al. 1994). This design uses the fact that the ocean has the long-term memory in the coupled atmosphere-ocean system. The fast adjustment of the atmosphere to the ocean variables such as the sea surface temperature (SST) and upper ocean heat content (HC) motivates the use of a steady-state statistical model for the atmosphere. A l l HCMs are based on the assumption that for monthly or longer time scales, contemporaneous correlation between wind stress and oceanic variables is associated with the atmosphere's rapid non-local adjustment to the oceanic anomaly patterns throughout the basin (Syu et al. 1995). In this chapter, we will develop two coupled models, i.e. one intermediate complexity dynamical ocean model with two different statistical atmospheric models. One atmo-spheric model is linear from linear regression and the other is nonlinear from neural networks. Our focus will be on the two most important issues affecting the performance 17 Chapter 2. The construction of two hybrid models 18 of the atmospheric models, i.e. the methodology in constructing the atmospheric model, and the choice of predictors in the atmospheric model. The response of the dynamical ocean to the two atmospheric models will also be investigated. This chapter is structured as follows: Section 2.1 briefly describes the dynamical ocean model and its performance in terms of ENSO simulations. Section 2.2 constructs the empirical atmospheric models, including linear and non-linear models. Section 2.3 compares the estimated wind stress using the linear and nonlinear models. Section 2.4 presents the ocean model forced with the wind stress from the two empirical atmospheric models. 2.2 Ocean M o d e l 2.2.1 The model l ing strategy The ocean model used in this research, one of intermediate complexity, originated from Anderson and McCreary (1985) and Balmaseda et al. (1994, 1995), but is extended to six active layers in this study. It consists of the depth averaged primitive equations in six active layers, overlying a deep inert layer. The model allows for an exchange of mass, momentum and heat at each layer interface by a parameterisation of entrainment. The model uses an Arakawa C grid layout, with a resolution of 1.5° x 1.5°, covering an extension from 30°N - 30°S in latitude and from 123°E - 69°W in longitude. The time step for integration is two hours. The boundaries are closed, with no-slip conditions. It has been shown to be a useful tool for studying the ENSO problem and can compare favourably with much more elaborate GCMs (Palmer and Anderson 1994). The model equations and some parameters are shown in Appendix A. The ability of the ocean model to reproduce the interannual variability of SST has been tested by forcing the ocean model with the observed F S U (Florida State University) Chapter 2. The construction of two hybrid models 19 observed wind stress (Goldenberg and O'Brien 1981) during the period 1961-1990. The wind stress used to force the model was computed from the pseudo-stress as r = /v_C0|U|U (2.1) where U is the wind speed (m/s), p a = 1.2 kg m 3 is the density of air, Ca is the drag coefficient, taken to have the value C_=0.0015 (Balmaseda et al. 1994), and fx is a parameter controlling the strength of the forcing. The model was first spun up for 100 years with monthly climatological wind stress and heat flux Qg as forcing fields. The climatological wind stress used consisted of the seasonal mean of the F S U stress, while the heat flux was represented by the Oberhuber (1988) heat flux Qo plus a relaxation term to T0, the observed SST, i.e. Qa = Q0 + X(T-T0), (2.2) where T is the model SST, Q0 and T 0 are the monthly climatological heat flux and SST, respectively, and A (which is negative) controls the rate of relaxation to the observed SST and varies with individual grids. The values of A are taken from the Oberhuber climatology (1988). After the 100-year spinup by the seasonal forcing, the model seasonal climatology was obtained. We then made a 30-year model control run, with forcing by the F S U wind stress from January 1961 to December 1990. The initial state was taken to be December of the model seasonal cycle. For this integration, the surface heat formulation was modified to Qs = £ o + A ( T m - T o ) + 0 . 2 A ( r - T m ) , (2.3) where Tm represents the model SST climatology. The factor 0.2 allows a weaker re-laxation to the model seasonal cycle, so that the model SST anomalies have similar magnitudes to those observed (Balmaseda et al. 1994). Chapter 2. The construction of two hybrid models 20 3 0 N 15N 30S 5 0 W 120W 9 0 W 6 0 W Figure 2.1: Location of the various oceanic regions used in the analysis. 2.2.2 The model l ing skills For easier comparisons and discussions hereafter, the ocean domain was divided into several standard regions, Nihol-4 and Eql-3 (Fig. 2.1). The very good performance of the model is seen from the statistics given in Table 2.1. In order to allow more direct comparison with the other models, we smoothed our time series by a 5-month running mean filter. The model has generally relatively good correlation skills in the whole equatorial Pacific except in the area of Ninol+2. In order to further quantify the performance of the model over the whole model domain, a point by point comparison between model and observations has been carried out by calculating the correlation and the RMS ( Root Mean Square) error, as in Miller et al (1993). Good correlation skills appear in the equatorial Pacific band of 15°S and 15°N, where the maximum correlation (0.6-0.8) is achieved over a relatively broad area in the central Pacific, decreasing smoothly westward and eastward (Fig. 2.2a). In the RMS error map (Fig. 2.2b), the largest error is located in the eastern Pacific, along the coast of South America. In the remainder of the basin, the RMS error is much smaller, with Chapter 2. The construction of two hybrid models 21 Table 2.1: Correlation between the observed SST and that from 10 models. The results of the first 9 models are taken from Palmer and Anderson (1994). The final model is the one used in this study. Model A , described in Wu et al (1994), is a l | layer model with specified mean climatology. Model B , is a 2|-layer model, described in Balmaseda et al (1994). Model C, described in Davey et al(1994), is also a 2|-layer model. While model D and E are versions of the G F D L Modular Ocean Model, with resolution of l | ° x l | ° and | ° latitude x l | ° longitude respectively. Model Region EQ3 Niho4 EQ2 Niho3 EQ1 Nihol+2 Cane-Zebiak - 0.46 - 0.60 - 0.68 Max-Planck Institute - 0.76 - 0.74 - 0.59 O P Y C - 0.72 - 0.63 - 0.46 G F D L - 0.81 - 0.69 - 0.57 A 0.43 0.64 0.77 0.73 0.69 0.54 B 0.27 0.59 0.77 0.75 0.69 0.40 C 0.34 0.55 0.67 0.62 0.51 0.26 D 0.60 0.76 0.79 0.55 0.55 0.54 E 0.59 0.76 0.79 0.65 0.60 0.58 Model used here 0.55 0.73 0.82 0.80 0.75 0.38 Chapter 2. The construction of two hybrid models 22 typical values of 0.4-0.6 °C. Compared with GCMs (Miller et at 1993), the ocean model presented here shares with the GCMs the high RMS error near the South American coast. In the western Pacific the RMS errors are smaller than those in the other models. Fig. 2.3 is the time evolution of the SST anomalies in the area of Nifio3, Niho4, and EQ2 from the model and from the observations. It shows that the model simulation is in very good agreement with the observations in these areas. 2.3 Atmospheric models 2.3.1 EOFs for predictors and predictands Two empirical atmospheric models were constructed: One was the traditional hnear regression (LR) model widely used in HCMs (Barnett et al 1993); the other was a non-linear regression model, a neural network. As potential predictors for both atmospheric models, several oceanic variables were chosen, namely the observed SST, the model SST and HC from the ocean model forced with the observed wind stress. The time period taken for the model construction was from 1964-1990, since in the first 3 years, the output of the ocean model was greatly affected by the ocean initial conditions and had poor agreement with observations. The observed SSTs were from the Comprehensive Ocean Atmosphere Data Set (COADS, Slutz et al. 1985; Smith et al. 1996). HC is defined here as the integral of the temperatures over the first two layers, calculated from H C i = T M HC = J2HCh (2.4) where Tr is the temperature of the bottom layer, and Hinit(i) is the initial thickness of layer i. The actual model depths in the two layers, climatologically monthly averaged over 30 year simulations, range from 89m to 144m in the first layer and from 134m to 188m in the second layer, with variations in season and in space. Chapter 2. The construction of two hybrid models 23 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Figure 2.2: Statistics for the model SST anomalies relative to the COADS SST data:(a) correlation (contour interval=c.i.=0.1) and (b) RMS error (c.i.=0.1°C). Chapter 2. The construction of two hybrid models 24 Figure 2.3: Time evolution of SST anomalies from the ocean model forced with F S U observed stress (sohd line) and observations (dash line) in regions EQ2, NIN03 and NIN04 from top to bottom. The Y-axis units are ° C. Chapter 2. The construction of two hybrid models 25 Table 2.2: The contributions of the first 10 principal components in % The EOFs mode Obs. SST Mod. HC Mod. SST •y 22.91597 8.233033 5.775219 4.890564 4.015489 3.757824 3.070947 2.894300 2.803019 2.457046 1 2 3 4 5 6 7 8 9 10 50.13675 60.17405 44.09614 11.31180 14.96291 14.42467 5.907232 10.23942 8.322747 4.483858 3.862239 3.926422 3.820368 1.903810 3.047551 3.321202 1.400397 2.473582 2.390177 1.121487 2.052638 2.323771 0.721858 1.838482 1.712472 0.616893 1.607294 1.529594 0.462901 1.179851 14.87040 11.42141 8.899261 7.412765 4.464581 3.539740 3.350200 3.117069 2.697852 2.283948 As in other studies (Barnett et al. 1993; Balmaseda et al. 1994), an E O F (Empirical Orthogonal Function) analysis was first applied to each dataset to extract the predictors and predictands. The oceanic predictor field T(x, t), and the predictand field r(x,<), the zonal or meridional component of the wind stress, were expressed by E O F analysis as where n is the mode number and the seasonal cycle had been removed for both fields prior to the E O F analysis. Table 2.2 is variance contribution explained by the first 10 E O F modes for the variables used in statistical atmospheric model. As shown in the table, for the model SST and HC, the first 3 E O F modes accounted for over 65% of total variance, whereas for the observed SST about 67% of total vari-ance was explained by the first 3 E O F modes. The variance contribution by individual modes became rather small after the first 3 modes. In contrast, the first 3 wind stress EOFs explained only about 35% of the total variance, due to presence of high frequency oscillations and noise in the wind stress field. The first three wind E O F modes captured the main low frequency signals, e.g. ENSO (Fig. 2.4), and are highly correlated with the T ( X , i ) = E «n(*K(x) T ( X , i ) = E /M<)/n(x) (2.5) n n Chapter 2. The construction of two hybrid models 26 64 66 68 70 72 74 76 78 80 82 84 86 88 90 64 66 68 70 72 74 76 78 80 82 84 86 88 90 0.4 r 64 66 68 70 72 74 76 78 80 82 84 86 88 90 64 66 68 70 72 74 76 78 80 82 84 86 88 90 0.4 r 64 66 68 70 72 74 76 78 80 82 84 86 88 90 64 66 68 70 72 74 76 78 80 82 84 86 88 90 Figure 2.4: The time series of the first E O F mode, the second E O F mode and third E O F mode from top to bottom. The left panel is for zonal wind and right panel for meridional wind. observed SST anomaly averaged over the NIN03 area. Hence, following the suggestions of Latif et al (1990) and Goswami and Shukla (1991), we used the first 3 E O F modes of oceanic variables T as predictors, and the first 3 E O F modes of zonal or meridional wind stress as predictands, in constructing both the linear and the non-linear models. Chapter 2. The construction of two hybrid models 27 2.3.2 Neural Network Model A feed-forward neural network ( N N ) is a non-parametric statistical model for extracting nonlinear relations in data. A common N N model configuration is to place between the input and output variables (also called 'neurons'), a layer of'hidden neurons' (Fig. 2.5). The hyperbolic tangent function is used as the transfer function (Bishop 1995, pl27). The value of the jth hidden neuron is yj = t a n h ( ^ wi]xi + bj), (2.6) i where X{ is the i-th input, W{j the weight parameters and bj the bias parameters. The output neuron is given by z = Yl™Jyj + l). (2.7) j A cost function J = ({z- zohsf) (2.8) measures the mean square error between the model output z and the observed values z0\,s. The parameters Wij, Wj, bj and b are adjusted as the cost function is minimized. The procedure, known as network training, yields the optimal parameters for the network. As in Tangang et al (1998a,b), steepest descent with momentum and adaptive learning rates was used during the optimization, without any constraints imposed. The random number generator in M A T L A B which generates uniformly distributed random numbers on the interval (0.0, 1.0) was used to initialize these parameters. For an N N with m i inputs and m2 hidden neurons, the number of model parameters is m i x m2+2 x m 2 + l . In this chapter, we set the number of hidden neurons to 3, so the N N has 16 parameters. The 3 input neurons were the first 3 E O F time series ctn(t) (either for SST or for HC), and the single output neuron was one of the (zonal or meridional) wind stress E O F time series 3n{t), i.e. a separate network was used to Chapter 2. The construction of two hybrid models 28 Input Layer Figure 2.5: A n example of a neural network model, where there are four neurons in the input layer, three in the hidden layer, and one in the output layer. The parameters Wij and Wj are the weights, and bj and b are the biases. The parameters bj and b can also be regarded as the weights for constant inputs of value 1. predict each of the wind stress modes. There was no time lag between the predictors and the predictand. One critical element in N N design is the size of the architecture, i.e., the number of hidden neurons. A smaller than optimal network fails to learn all the relevant information from the data, i.e. underfits, whereas a larger than optimal network often ends up learning the noise in the data, i.e. overfits. The strategy we employed to determine the number of hidden neurons was based on cross-validating the correlation skills of the network. Chapter 2. The construction of two hybrid models 29 A cross-validation procedure involves dividing the data record into several segments, selecting one segment for the test data and the rest for training data. The N N model is built using the training data only, and model predictions are tested on the test data. Next, another segment is selected as the test data, and a new model version is built. This is repeated until the entire data record has been used for testing. The cross-validation scheme ensures that no training data are used for testing the prediction skills, hence the artificial skill associated with the overfitting problem (von Storch and Zwiers, 1999) can be effectively eliminated by the cross-validation scheme. In this study, the record of 1964-1990 was divided into 3 segments, the first 7 years and the remaining two 10-years periods. We tested different N N models by increasing the number of hidden neurons. Cross-validation skills (correlation and RMS error) were monitored for each network. As the number of hidden neurons increased starting from 1 to 3, the cross-validation skills were enhanced. These skills were little changed as the number of hidden neurons was increased to 4-6. But the skills were greatly degraded with further increase in the number of hidden neurons. Hence we used 3 hidden neurons in the atmospheric model. In the case where both SST and HC are used as predictors, the network had 6 input neurons but still only 3 hidden neurons. Cross validation was also used to choose the the error goal (i.e. the level of J to terminate the optimization) and the maximum number of iterations (epochs) in the training process. Finally, an ensemble of 25 NNs with random initial parameters were used and the final output of the N N model was actually the ensemble average of the 25 individual model outputs. Ensembles can greatly alleviate the instability problems associated with N N modelling (Hsieh and Tang, 1998). Chapter 2. The construction of two hybrid models 30 2.3.3 Linear model The linear regression model is similar to that of Barnett et al (1993). Assume that the wind stress anomalies r are a linear response to T, we have T(x, t) = Y, « n ( * ) e n ( x ) r(x, *) = £ /3„(*)/„(x), (2.9) Pn = J2Cmai (2-10) i=l where T, CK, /3, e, m and n have the same meaning as in (2.5); and i is a mode number. The matrix of regression coefficients (C) relating the two datasets is defined by < Cti[3n > (2.11) where < ... > denotes a time average. The statistical estimate of the surface wind- stress field (r) is obtained from ai(t) = Y/T{x,t)ei{x) (2.12) M) = TlCniai{t) (2-13) i and finally r(*,0 = E&(0/nOO (2-14) n 2.4 Results from atmospheric models 2.4.1 The predictors Most HCMs used either simulated SST from ocean models or observed SST to estimate the wind stress (Barnett et al. 1993; Syu et al. 1995). Whether SST is the best predictor for the wind stress is debatable. Observations show that the SSTs do not reflect the changing subsurface temperatures in the tropical western Pacific, where subsurface temperature anomalies and thermocline displacements have an important role in the Chapter 2. The construction of two hybrid models 31 Table 2.3: Cross-validated correlation between the predicted wind stress E O F time series and the observed wind stress E O F time series for the first 3 modes, using the observed SST, the model SST, the model HC, and the model HC+SST as predictors. Results are shown for both the N N model and the L R model, and for the zonal (x) and meridional (y) components of the wind stress. Predictors Obs SST Mod.SST HC HC+SST Predictands N N L R N N L R N N L R N N L R rx ft 0.74 0.72 0.81 0.81 0.89 0.89 0.88 0.88 rx ft 0.57 0.52 0.54 0.52 0.76 0.73 0.75 0.70 rx ft 0.38 0.24 0.11 0.06 0.40 0.35 0.29 0.26 Ty ft 0.86 0.83 0.83 0.81 0.86 0.86 0.84 0.86 Tv ft 0.53 0.43 0.47 0.46 0.66 0.66 0.57 0.62 Ty ft 0.47 0.45 0.21 0.19 0.24 0.21 0.28 0.24 ocean-atmosphere couphng processes (e.g. White and Pazan 1987, Latif and Graham 1992). Furthermore, SST is noisier than the HC, probably resulting in the wind stress being less well correlated with SST than with HC. In addition, as the ocean memory resides in the thermocline depth, HC is a more stable variable than SST, leading to HC having a higher correlation with wind stress. Therefore, using the upper ocean heat content (HC, defined in Eq. 2.4) as predictors might produce higher skills than SST. The first E O F modes for the HC, the zonal and meridional wind stress (Fig. 2.6), are the modes associated with the ENSO oscillation-where the HC anomaly shifts east-west along the equator (Fig. 2.6a), the zonal wind anomaly develops in the western equatorial Pacific (Fig. 2.6b), and the trade winds show anomalous convergence along the equator (Fig. 2.6c). Because of their ENSO nature, the anomalies in these E O F modes are all mainly confined to within 15°N - 15°S. To explain the anomalies outside this narrow equatorial belt, the second, third, and even higher modes are needed. Table 2.3 shows the cross-validated skills attained by the N N and L R models with the observed SST, model SST, model HC, and model SST+HC serving as predictors for the Chapter 2. The construction of two hybrid models 32 3 0 N P 30S 30N 15N 15S 30S 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Figure 2.6: The first E O F mode for (a) the heat content (HC) from the ocean model driven by the F S U wind stress, (b) the F S U zonal wind stress and (c) the meridional wind stress. Negative contours are shown as dashed curves and zero contours as dotted curves. Contour interval = 0.01. Chapter 2. The construction of two hybrid models 33 wind stress E O F time series fln (n = 1,2,3). Here SST + HC does not mean combined EOFs, but that their separate E O F time series normalized by the standard deviations are the predictors. HC as predictors generally had the highest skills as expected, whereas both observed SST and model SST generally attained the lowest skills. SST+HC predic-tors attained lower skills than HC alone. This suggests that the first 3 EOFs of HC well represent the ocean status, and more SST EOFs as input only bring additional noise and overfitting. For the first zonal wind stress mode, the model SST actually did better than the observed SST, probably because the ocean model acted as a complicated space-time filter, thereby removing some noise in the SST (Latif and Graham 1992). For the other wind stress modes, the model SST did not do as well as the observed SST. In general, N N did not predict the first zonal or meridional wind stress mode much better than LR; only for the second and third modes did N N seem to have an edge over LR. Some small differences between LR and N N can be seen in Fig. 2.7, where the N N curves almost overlap the straight lines, indicating weakly nonlinear relations. 2.4.2 Prediction skills of the N N and L R models Since we only tried to predict the first 3 wind stress modes, the prediction skills were compared in later sections against an idealized wind stress -the wind stress reconstructed from only the first 3 E O F modes of the FSU wind stress. Such a comparison more objectively evaluates the skill of the atmospheric model, as it excludes the noisy higher modes which are not modelled. The correlation map between the idealized wind stress and the F S U wind stress shows that the idealized field is a reasonable representation, especially in the western and central Pacific (Fig. 2.8). As the model HC has the best prediction skills for the wind stress, henceforth we will only use the HC as predictor. The reconstructed stress field was obtained by multiplying the predicted E O F time series (either by N N or LR) to the first 3 E O F modes. The Chapter 2. The construction of two hybrid models 34 Figure 2.7: (a) The first mode of the zonal wind stress plotted as a function of the third mode of the SST, and (b) the second mode of the zonal wind stress versus the first mode of the SST. The observations are shown as dots; the simulated zonal wind stress mode from the NN model is indicated by the small circles and from the LR model by the straight lines. Chapter 2. The construction of two hybrid models 35 a 150E 180E 150W 120W 9 0 W Figure 2.8: Correlations of the F S U wind stress and the idealized wind stress (recon-structed from the first 3 E O F modes of the F S U wind stress) for (a) the zonal stress and (b) the meridional stress. Contour interval = 0.1. Chapter 2. The construction of two hybrid models 36 cross-validated correlation and the sum of squared error (SSE) of the reconstructed stress field from N N model verified against the idealized field are shown in Fig. 2.9 and Fig. 2.10. As seen in Fig. 2.9, for zonal stress, the best skills occurred in the equatorial western and eastern Pacific, whereas the worst occurred off the coast of South America, in the Australian summer monsoon region and in the subtropical North Pacific. The high skill areas are associated with the anomalous SST areas during ENSO events where the response of zonal stress to the ocean status is strong due to active couphng. The lower correlation along the coast of South America might be attributed to the fact that wind-stress is almost 'white' in this region (Goldenberg and O'Brien, 1981, Latif and Flugel 1991). For meridional stress, the highest skills were found in the Intertropical Convergence Zone (ITCZ) and the South Pacific Convergence Zone (SPCZ) areas. The map of SSE (Fig. 2.10) indicated that the estimation of the amplitude by the N N model was good, especially in eastern Pacific ocean. Large errors only occurred at the ITCZ and SPCZ areas, in contrast to the correlation map (Fig. 2.9), where these two areas have good skills (especially for the meridional stress). The active couphng in these areas induced large anomalous variations in the stress, which generated large amphtude errors even though the phase errors were small, producing good correlations but large SSE. Differences between the prediction skills of the N N model and the L R model for the period of 1964-1990 are shown in Figs. 2.11 and 2.12. The correlation skill differences (Fig. 2.11) between N N and L R were very small, though N N skills were indeed slightly ahead of L R skiUs in most areas. For the zonal stress (Fig. 2.11a), the N N model outper-formed L R in almost the whole subtropical domain of 15°N - 30°N and near the Niho3 region. That the improvements occurred in these regions can be understood from our earlier finding (Table 2.3) that the N N and the L R had the same skill in predicting the mode 1 zonal stress (Fig. 2.6), but the N N had an edge over L R for modes 2 and 3. Hence Chapter 2. The construction of two hybrid models 37 a I ' J > . < . .1. - . — 'I • - •-• ^ •'• .1 - -1 5 0 E 1 8 0 E 1 5 0 W 1 2 0 W 9 0 W Figure 2.9: Cross-validated anomaly correlation between the predicted wind stress by N N (with model H C as predictors) and the idealized wind stress (i.e. the F S U wind stress with only the first 3 E O F modes): (a) zonal stress, and (b) meridional stress. Contour interval=0.1. The seasonal cycle has been removed prior to the analysis, and in all following maps. Chapter 2. The construction of two hybrid models 38 Figure 2.10: The sum of squared error (SSE) of the predicted wind stress by NN verified against the idealized wind stress: (a) zonal stress, and (b) meridional stress. Contour interval=0.2 N 2 m~4. Chapter 2. The construction of two hybrid models 39 only in regions outside the main anomaly area of mode 1 (i.e. the western equatorial Pacific in Fig. 2.6) would the N N appear to have slightly better skills than LR, as found in Fig. 2.11a. For the meridional stress (Fig. 2.11b), N N did slightly better in the eastern Pacific away from the equator. The difference of the SSE between the two models (Fig. 2.12) indicated that the N N model slightly outperformed the L R model for much of the domain except for the western equatorial region of 150°E - 160°W centred at 5°S. The slight advantage of N N over L R is manifested more clearly in the error map (Fig. 2.12) than in the correlation map (Fig. 2.11), suggesting that model nonlinearity may be slightly more useful in estimating the amplitude than in estimating the phase of the wind stress anomalies. The very small skill differences between N N and L R follows from the fact that the equatorial dynamical system is almost linear, so a nonlinear model does not give much better results than a linear model. Tang et al. (2000) found that with sea level pressure (SLP) as predictors for the SST anomalies, N N slightly outperformed L R in the Niho3 region, but not so in the Niho4 region, suggesting that nonlinearity is quite weak in the eastern-central equatorial Pacific, but even weaker in the western-central equatorial Pacific. Here, Fig. 2.12 and to a lesser extent Fig. 2.11 are consistent with the Tang et al. (2000) finding. 2.5 The ocean model driven by the empirical wind 2.5.1 SST skills To assess the effect of the empirical wind stress from the N N and the L R models on the ocean model, we ran the ocean model twice, with forcing by the two model predicted wind stresses during the period 1964-1990. The outputs of the ocean model forced by the idealized wind stress (i.e. the first 3 EOFs of the FSU data) were later used to verify Chapter 2. The construction of two hybrid models 40 a 150E 180E 150W 120W 90W b . . ,l, ,) 1 i - .r — .Q.02T- . . . L ... 150E 180E 150W 120W 90W Figure 2.11: The correlation between the wind stress predicted by the N N model and the idealized stress minus the correlation between the wind stress predicted by the L R model and the idealized wind stress, for: (a) zonal stress, and (b) meridional stress. Contour interval = 0.02. The zero contours are the dotted curves, while the negative contours are the dashed curves. Positive values means the N N predicted wind stress is outperforming the L R one. Chapter 2. The construction of two hybrid models 41 Figure 2.12: The SSE ofthe N N wind stress minus the SSE ofthe L R wind stress, (both verified against the idealized wind stress), for: (a) zonal stress, and (b) meridional stress. Contour interval=0.02 N 2 m - 4 . Negative values means the N N predicted wind stress is outperforming the L R one. Chapter 2. The construction of two hybrid models 42 the skills from the empirical wind stress. Fig. 2.13 compares the SST and HC from the ocean model driven by the idealized wind stress with that driven by the full F S U stress, showing a generally close relation, especially in the western and central equatorial Pacific. This justifies the use of the ocean model driven by the idealized stress as the standard for comparing the empirical winds. Fig. 2.14 shows the skills of the SST from the ocean model forced with the empirical wind stress from the N N model. As seen in the correlation map (Fig. 2.14a), the skill is good, with a correlation of over 0.8 covering much of the model domain. The highest skill, up to over 0.9, occurred at the central equatorial Pacific region. A comparison of ocean models driven by the empirical wind stress from the N N model and that from the L R model was then made. Fig. 2.15a depicts the difference in the model SST correlation skills between the N N model and the L R model when they were each verified against the standard SST. For the equatorial western and central Pacific, the difference was negligible. As mentioned earlier, this is probably due to the mainly hnear dynamics in the equatorial western and central Pacific. The largest differences occurred in the off-equatorial areas and in the eastern equatorial Pacific, where the N N winds tended to outperform the L R winds. The maximum correlation difference reached 0.34 in the northwest region around 150°E and 15°N. The positive correlation skill areas in the off-equatorial regions of Fig. 2.15a roughly coincided with the positive zonal wind stress skills attained by the N N over LR in Fig. 2.11a. There was no such agreement between Fig. 2.15a and Fig. 2.11a in the eastern equatorial Pacific, where the skill differences between N N and L R were small for the zonal stress, but relatively large for the SST. This could be due to the fact that in equatorial western Pacific, the oceanic response is mainly locally forced by the wind stress, whereas in the eastern Pacific, equatorial Kelvin wave propagation, upwelling and vertical mixing are thought to predominantly control the SST variabihty (Battisti 1988). Chapter 2. The construction of two hybrid models 43 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Figure 2.13: A comparison of the ocean model driven by the F S U wind stress and that driven by the idealized stress (i.e. with only the first 3 EOFs), (a) SST correlation be-tween the two models (c.i.=0.1) and (b) SST RMS error (c.i.=0.05°C); (c) HC correlation (c.i.=0.1) and (d) HC RMS error (c.i.=0.005°C). Chapter 2. The construction of two hybrid models 44 a 150E 180E 150W 120W 9 0 W b . . .1- . i j i ' > • • 150E 180E 150W 120W 9 0 W Figure 2.14: A cross-validated comparison of the ocean model SST between that driven by the idealized wind stress and that driven by the empirical wind stress from the N N model (with H C as predictors), by (a) correlation (contour interval=0.05), and (b) R M S error (contour interval=0.1°C). Chapter 2. The construction of two hybrid models 45 The difference in the SST RMS error between the ocean models driven by the N N and L R wind stress (Fig. 2.15b) generally agreed with Fig. 2.15a, i.e. higher correlation skiU areas corresponded with lower RMS error regions, and vice versa. In summary, using a nonlinear empirical atmospheric model to drive the ocean might bring modest benefits for SST simulation in the off-equatorial tropical regions and in the eastern equatorial Pacific. Improvements in the equatorial western and central Pacific would be unlikely as the dynamics in these regions appeared very nearly hnear from other studies (e.g. Tang et al. 2000). 2.5.2 Heat content Skills The HC redistribution in western tropical Pacific can lead to the evolution of SST anoma-lies in eastern Pacific and has been known to be an important factor in the evolution of many ENSO episodes. In particular, the HC anomalies over the equatorial band 5°N to 5°S can be a very good precursor for the SST anomalies in the Nino3 region (Zebiak 1989; Latif and Graham et al. 1992; Balmaseda et al. 1994). We therefore examined the HC in the ocean model forced with the wind stress from the N N model and from the L R model. The HC from the ocean model forced with the idealized stress was taken as the verification standard. The HC skills from the ocean forced with the N N model wind stress was generally better than those with the LR stress for the western Pacific basin, in particular in the western and central Pacific over the equatorial band of 10°N to 10°S (Fig. 2.16). From the ocean model driven by the idealized stress, the HC anomalies averaged over the whole equatorial Pacific (124 oE-70°W,5 oN-5°S) led the observed Nino3 SST anomalies by few months (Fig. 2.17). Their correlations were 0.73 and 0.72 with the HC leading by 3 and 4 months respectively. The 3-month lag correlations of the observed Nifio3 SST anomahes and the HC index averaged over whole equatorial Pacific for the ocean model forced by Chapter 2. The construction of two hybrid models 46 a . ———rr. .< J .'. .r< \ . ' . ? . . .' i . . . t . . 150E 180E 150W 120W 90W Figure 2.15: The skill differences in the model SST between the model driven by the NN model wind stress and that driven by the LR model stress, with both model SSTs verified against the standard SST, i.e. the SST from the model driven by the idealized wind stress: (a) correlation skill difference (Contour interval = 0.02), and (b) RMS error difference (contour interval=0.02°C). Positive regions in (a) indicate NN ahead of LR, while negative region in (b) indicate NN ahead of LR. Chapter 2. The construction of two hybrid models 47 the N N wind and by the L R wind did indicate the N N wind to have slightly better skill (0.62 for N N and 0.58 for LR). In summary, this chapter investigated the possibility of using a nonlinear empirical atmospheric model for hybrid coupled modelling by developing a neural network (NN) model for predicting the contemporaneous wind stress field from the ocean state, and comparing the N N model with a linear regression (LR) model. The oceanic response to the two atmospheric models were also examined. The N N model generally had slightly better skills than the simple LR model in the off-equatorial tropical Pacific and in the eastern equatorial Pacific. Compared with the oceanic response to the L R atmosphere, the oceanic response to the N N atmosphere generated slightly better SST simulation skills in the off-equatorial tropical Pacific and in the eastern equatorial Pacific, and better H C simulation skills in the western and central equatorial Pacific, Chapter 2. The construction of two hybrid models 48 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Figure 2.16: Cross-validated skill differences in the model HC between the ocean model driven by N N model wind stress and that driven by the L R model stress, both verified against the model HC driven by the idealized wind stress: (a) correlation difference (c.i.=0.02), (b) RMS error difference (c.i.=0.05°C). Positive regions in (a) indicate N N ahead of LR, while negative regions in (b) indicate N N ahead of LR. Chapter 2. The construction of two hybrid models 49 Figure 2.17: Time evolution of observed SST anomalies in NIN03 (solid curve) and model HC averaged over 124°E-70°W, 5°N-5°S (dot-dashed curve) forced by the F S U observed wind stress. Both curves were normalised and smoothed by a 3-point running mean. Chapter 3 Hybrid coupled models of the tropical Pacific — Interannual variability 3.1 Introduction During the last decades, the study of a variety of coupled models, including simple models (e.g., Suarez and Schopf 1988; Hirst 1986), intermediate coupled models (e.g., Zebiak and Cane 1987; Anderson and McCreary, 1985), hybrid coupled models (e.g., Barnett et al 1993; Balmaseda et al. 1994,1995; Neelin 1990), and fuUy coupled general circulation models (GCM) (e.g.,Philander et al. 1992; Latif et al. 1993), led to a better description of the phenomenon and understanding of the underlying cychc nature of ENSO (Schopf and Suarez 1988; Battisti and Hirst 1988; Jin et al. 1994; Tziperman et al. 1994; Chang et al. 1995,1996). As emphasized in previous chapters, the hybrid coupled model has been widely applied in the understanding and prediction of ENSO (Barnett et al. 1993; Davey et al. 1994; Chang et al. 1995; Syu et al. 1995; Balmaseda et al. 1994, 1995; Blanke et al 1997; Eckert and Latif 1997). An important aspect affecting the H C M performance is the construction ofthe empirical atmospheric model, e.g. what method was used for estimating the surface wind stress field from a given ocean state. A l l empirical atmospheric models used in HCMs have so far been linear statistical models. In Chapter 2, the possibility of improving the empirical atmospheric model by a nonlinear approach using neural networks (NN) was investigated. However, in Chapter 2, the ocean model was not coupled with atmospheric model and all simulations were driven from uncoupled experiments. This chapter is an 50 Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 51 extension of Chapter 2. In this chapter, we couple the ocean and atmosphere, examine the dynamical behaviour of the H C M with a nonlinear atmosphere (henceforth the N H C M ) , and a hnear atmosphere (LHCM), and investigate the differences between N H C M and L H C M . This chapter is structured as follows: Section 3.2 briefly describes the coupling strat-egy. Section 3.3 displays the performances by the N H C M and compares N H C M with L H C M . Section 3.4 discusses the dynamical behaviour of both coupled models. 3.2 The Coupl ing strategy After an atmospheric model is constructed, a coupling between the atmosphere and the ocean can be implemented. The SST anomalies from the ocean model are projected onto the first three E O F modes to extract the predictors. An atmospheric model (LR or NN) is then performed to reconstruct the zonal wind stress r_, and a similar one for the meridional component r y . As common with statistical models, the variance of the predicted variables are lower than the variance in the observed variables, hence the wind-stress estimates are scaled up to their observed variance by an adjusting scale factor of 1.2, as in Barnett et al. (1993). The stress anomalies are further multiplied by a scalar parameter o, i.e. the relative coupling parameter, to examine the impact of the coupling strength on ENSO oscillations. The reconstructed wind stress anomalies are added to the climatological wind stress field to force the ocean model. The couphng interval is 15 days as in Balmaseda et al. (1994, 1995). Fig. 3.1 demonstrate the schematic diagram of coupling process. The ensemble of atmospheric models derived with the cross-validation scheme is used in the couphng. Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 52 Given initial conditions for ocean model Extract EOF time coefficients for SST and HC (to be used as predictors in the atmos. model) Run Atmos. Model to predict wind stress + Climatological wind stress f Drive ocean model with wind stress F i g u r e 3.1: A schemat ic d i ag ram of coup l ing s trategy for the coup led models . Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 53 3.3 Simulations from the standard coupling experiments In the section, we present results of the standard coupling experiments. Standard cou-pling is defined as having the coupling coefficient S equal to unity. Two coupled models, the dynamical ocean with either a nonlinear atmosphere or a linear atmosphere, are in-vestigated. Each coupling begins from the end of the ocean climatology run of 100 years and is integrated forward for 100 years. As the two models require in general about 3 years to reach their equilibrium, the integration of first 3 years will be removed in the discussions. The seasonal cycle has been removed prior to E O F and P O P (principal oscil-lation pattern) analyses in order to explore the model's interannual variability. We also use the subsurface simulations from the control ocean model forced with the observed F S U wind stress (called 'idealized' data) instead of actual observed subsurface data as a standard to verify the coupled model's subsurface behaviour. 3.3.1 N H C M simulations Fig. 3.2 is the N H C M SST climatology averaged over the integration of 97 years and the observed climatology during 1961-1990, showing that the N H C M can well generate realistic climatology in both spatial structures and amplitudes although its 'warm pool' in the equatorial western Pacific and 'cold tongue' in the equatorial eastern Pacific are both slightly stronger than observed. The N H C M also captured well the pronounced seasonal cycle of SST in the eastern Pacific, although the simulated amplitudes are smaller than the observations (Fig. 3.3a). The seasonal cycle of the zonal wind stress also compares well with the observations (Fig. 3.3b). Note that the atmospheric model only calculates the anomalous wind stress, which is added to the observed climatological wind stress to drive the ocean model. Hence the good climatology simulation for the ocean with N H C M comes mainly from the contributions of the uncoupled response of the ocean to Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 54 the observed climatological winds. In a coupled model with an anomaly atmospheric component, the contributions from the coupled feedbacks to the model climatology are generally much smaller than from the uncoupled response of the ocean to the observed climatology winds. With the seasonal cycle removed, Fig. 3.4 compares the two leading E O F modes of the modelled SST from N H C M with those derived from the observed SST during 1964-1990. The observed SST field is taken from the COADS. The first two EOFs of the modelled SST account for 53% and 24% of the total variance respectively, which compares with 50% and 11% of the total variance accounted for by the two leading modes from the observations. As seen in Fig. 3.4, the leading EOFs for the simulated SST (upper panel of Fig. 3.4) generally resembled the observed modes, except that the simulated modes appear to be more narrowly confined to the equator, with less variability at the eastern boundary. These are common defects in all HCMs, and even in coupled G C M models (Barnett et al 1993; Chang et al 1996)— the weak eastern boundary variability is often blamed on the poor parameterization of vertical mixing in the eastern equatorial Pacific Ocean. The first two E O F modes of the simulated zonal wind stress and those from the observed F S U zonal wind stress (Fig. 3.5) reveal that the structures of simulated modes are in good agreement with those from observations. The variance accounted for by the first two modes of the N H C M modelled zonal wind stress is 81% and 10% respectively. The major dynamical processes associated with the ENSO cycle in the tropical Pacific ocean are the westward off-equatorial propagation of information and the eastward prop-agation of information along the equator (Barnett et al 1993, Philander 1990). These information propagations mainly result from the propagation of Rossby waves and equa-torial Kelvin waves in the tropic Pacific. In order to explore the main features of propa-gating modes in the coupled model, a POP analysis has been made using the 3 leading Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 55 Figure 3.2: (a) The SST chmatology of N H C M from the run of 97 years (c i . = 1°C) and (b) The observed SST chmatology during 1961-1990 (c i . = 1°C) and (c) the difference, i.e. (a) minus (b). (C.i.= 0.5°C). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 56 150 200 250 150 200 250 Longitude Longitude Figure 3.3: Time-longitude distribution of the mean seasonal cycle of (a) SST and (b) zonal wind stress, repeated once with the annual mean removed, for the N H C M . Negative contours are shown as dashed curves and zero contours as dotted curves. C.i.=0.2°C for (a) and 0.02 N m " 2 for (b). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 57 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Figure 3.4: (a) E 0 F 1 of SST from the N H C M , (b) EOF2 of SST from the N H C M , (c) EOF1 of the observed SST and (d) EOF2 of the observed SST. Negative contours are shown as dashed curves and zero contours as dotted curves. C i . = 0.01°C for (a),(b) and (d), and c i . = 0.005°C for (c). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 58 30N 15N 15S 3 0 S 30N 15N 150E 180E 150W 120W 90W 15S 30S 150E 180E 150W 120W 90W 30N 15N 15S 30S 30N 150E 180E 150W 120W 90W 15S 30S 150E 180E 150W 120W 9 0 W Figure 3.5: (a) E 0 F 1 of the zonal wind stress from the N H C M , (b) E 0 F 2 of the zonal wind stress from the N H C M , (c) E 0 F 1 of the FSU zonal wind stress and (d) E 0 F 2 of the F S U zonal wind stress. C i . = 0.01 N m " 2 for (a), (c) and (d), and c i . = 0.005 N rn"2 for (b). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 59 combined E O F modes of three fields from the coupled model— SST, upper ocean heat content (HC) and zonal wind stress. Several measures, e.g. the oscillation period T, the decay time scale and the variance contribution, as discussed in Wu et al. (1994), can be used to pick out useful or significant P O P modes from the POP analysis. The ENSO-related mode is found to explain 32% of the total variance, with a period of 57 months. A maximum correlation of 0.99 between the imaginary part and the real part of the POP coefficient occurs at a lag of a quarter of the period (15 months), as expected. The spatial patterns of the P O P mode are shown in Fig. 3.6. The POP technique interprets the pattern as evolving cyclically in time: Q P -»• - Q -> - P -> Q ... The characteristics ofthe POP patterns are very similar to those found by Balmaseda et al. (1994), Latif et al (1993), and Barnett et al (1993). In the peak phase of an ENSO warm event (pattern P), the warm water present in the equatorial central and eastern Pacific ocean yields the warm SST and HC anomalies in the region. A strong zonal HC gradient at the equatorial central Pacific weakens the upwelhng there and intensifies the warm Kelvin waves propagating eastward. From the warm SST, atmospheric convection occurs, resulting in a convergence of mass and heat in the atmosphere on both sides of the equator, which enhances the oceanic upwelhng in the off-equatorial regions and induces the cold westward propagating Rossby waves. The warm eastward propagating Kelvin waves bring warm waters to the eastern Pacific ocean to further intensify the anomahes; while the westward propagating Rossby waves are reflected at the western boundary as cold Kelvin waves to propagate eastward. When the cold Kelvin waves arrive at the central and eastern Pacific ocean to weaken the warm anomahes in the region, the transition phase of ENSO (warm-to-cold phase) appears as pattern —Q. The zonal wind stress anomahes are closely related to the information propagation in the ocean, which can be demonstrated using Fig. 3.6e and Fig. 3.6f. During the peak Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 60 30N 15N 0 15S 30S 150E 180E 150W 120W 90W C v4i rf.V ° ° 0 5 _ •'• / " : 150E 180E 150W 120W 90W e 150E 180E 150W 120W 90W d 150E 180E 150W 120W 90W f 150E 180E 150W 120W 90W 30N 15N 0 15S 30S -0. — ~N, — ^ \ ^ \ \ c \ ^ —" ^ _ • . '-o.ol D.-005-- /" — — 150E 180E 150W 120W 90W Figure 3.6: Spatial patterns of the ENSO-related POP mode for N H C M . Patterns P (real) are shown on the left, and patterns —Q (imaginary) appear on the right. SST anomalies are in panels (a) and (b), HC anomahes are in panels (c) and (d), and zonal wind stress are in panels (e) and (f). C.i.= 0.005°C for (a),(b),(c) and (d); C.i.= 0.005 N rn"2 for (e) and (f). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 61 phase of ENSO, large westerly wind anomahes prevail over the equatorial central Pacific associated with the eastward propagating equatorial Kelvin wave; whereas during the transition phase ( — Q), the westerly wind anomahes in the equatorial central Pacific are much weaker, and easterly wind anomahes dominate over much of the tropic Pacific. The wind stress anomahes in turn act on the ocean. In summary, the coupled model with nonlinear atmosphere can produce a realistic climatology; and major features of ENSO variabihty in the tropic Pacific can be weU captured. The ENSO dynamics shown in the coupling simulation can be explained by the delayed-action oscillator. 3.3.2 Some differences between N H C M and L H C M The same analysis procedure is also performed for L H C M . The L H C M chmatology is basically in good agreement with the N H C M chmatology described above, i.e. the L H C M also produces a realistic mean structure. In this subsection, we will emphasize some differences between the L H C M and the N H C M in terms of their chmatology and ENSO features. Some studies have shown that ENSO osciUations in the H C M can be sensitive to small changes in the construction of the atmosphere model. For example, Syu and Neelin (2000) found noticeable sensitivity in the simulated ENSO cycles using a slightly shorter or longer data period in estimating the SVD model, in particular for the coupled period. However, as mentioned in Chapter 3.2, an ensemble average over individual atmospheric models derived from the cross-validation scheme (Chapter 2) should to some extent alleviate the sensitivity of ENSO characteristics to the data record. One difference is shown in the simulation of the equatorial cold tongue of the eastern Pacific. Fig. 3.7 is the L H C M SST climatology and its difference with the N H C M ch-matology. The L H C M produces a warmer SST climatology at central equatorial Pacific than observations. Fig. 3.7c shows that the N H C M is warmer than L H C M in the western Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 62 equatorial Pacific and cooler in the central-eastern equatorial Pacific. Some differences are also found in the subsurface ocean temperature at the depth of around 100m. Fig. 3.8 shows the errors of subsurface temperature climatology of L H C M and N H C M verified against the idealized subsurface temperature. As shown in Fig. 3.8, the errors of the L H C M are generally larger than those of the N H C M , reaching about a factor of 2 in the equatorial eastern Pacific. However, as the errors are generally small, these differences between the climatology mean of the two models are not very significant. The spatial distribution of the SST differences between N H C M and L H C M is in notable contrast to that between the uncoupled ocean responses to the two atmospheric models shown in Chapter 2. For example, the SST difference between the uncoupled ocean response to the linear atmosphere and that to the the nonlinear atmosphere is remarkable in the off-equatorial areas and the equatorial eastern Pacific (Fig. 2.15). A notable difference between L H C M and N H C M occurs in the upper-ocean thermal structure. Fig. 3.9 shows the Hovmoller diagrams of heat content anomalies along the equator for N H C M and L H C M during year 31 to year 60 in the standard coupling run of 100 years. For comparison, the time-longitude diagram of the heat content anomalies taken from a control integration from 1961-1990 is also presented (Fig. 3.9b). The most pronounced feature in Fig. 3.9a and Fig. 3.9c is the regular oscillation structure in equa-torial Pacific, especially in the eastern Pacific. Fig. 3.9a and Fig. 3.9c also demonstrate some differences between L H C M and N H C M : (1) L H C M has a longer oscillation period than N H C M , with about 5 years for N H C M and 7-8 years for L H C M ; (2) the western Pacific has a much stronger oscillation in L H C M than in N H C M ; and (3) For the L H C M , often when there is a strong anomaly in the eastern and central equatorial Pacific, there is an anomaly of comparable strength but opposite sign on the western side, which is not observed for ENSO. Compared with Fig. 3.9b, L H C M does not capture ENSO features as well as N H C M . These differences between N H C M and L H C M in terms of ENSO can also Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 63 Figure 3.7: (a) The SST chmatology of L H C M from the run of 97 years (c i . = 1°C), (b) The SST chmatology of L H C M minus the observed SST chmatology during 1961-1990 (c i . = 0.5°C) and (c) The SST climatology of N H C M minus the SST climatology of L H C M (c i . = 0.1°C). Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 64 Figure 3.8: The model error of subsurface ocean temperature at the depth of around 100m verified against the idealised subsurface temperature for (a) L H C M and (b) N H C M . Negative contours shown as dashed curves denote that the subsurface model temperatures are lower than the 'idealized data', c i . = 0.02°C. Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 65 be seen in the structure of SST anomahes (Fig. 3.10). The structure of SST anomalies in N H C M shows anomalies starting in the central equatorial region, spreading eastward and westward (Fig. 3.10a), with a period of about 4-5 years. The lack of interannual signals in the west in N H C M is consistent with observations. In contrast, the L H C M SST anomalies are able to propagate to the western Pacific, with a period of about 7-8 years (Fig. 3.10c). A P O P analysis was performed for the L H C M simulation as was done for the N H C M . Fig. 3.11 is the leading POP mode, which explains 29% of the total variance. This mode has a period of 87 months, much longer than the period of the N H C M and realistic ENSO oscillations. During the peak phase, the equatorial SST anomahes change sign at 170°W (Fig. 3.11a), in contrast to those for the N H C M , where anomalies of the same sign reach 150°E (Fig. 3.6a). During the transition phase, the SST and HC anomahes are very weak in the eastern equatorial Pacific unlike the structures found in N H C M (Fig. 3.6). For the — Q pattern, the easterly anomalies in the off-equatorial areas in L H C M are slightly weaker than the anomalies found in N H C M . Similarly, during this transition phase, the westerly anomalies in the equatorial central Pacific do not withdrawn so quickly as in N H C M , with the strong westerly anomahes still dominating the equatorial central Pacific in the L H C M . The POP analysis with observed data shows that during the transition phase, the whole domain of the central equatorial Pacific is almost dominated by easterly anomalies (Latif and Fhigel, 1991). The withdrawal rate of the westerly anomalies in the model may influence the period of the oscillation— the slow withdrawal may explain why L H C M has a longer period than N H C M . An important feature of the ENSO oscillation is its phase-locking, i.e. the peaks tend to occur during a particular season. The phase-locking might be mainly due to the non-linear interactions between the seasonal cycle and self-sustained interannual osciUations. The seasonal cycle plays a central role and often dominates the strength of phase-locking Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 66 Figure 3.9: Time-longitude diagrams of ocean heat content anomalies along the equator from (a) N H C M , (b) control run forced with observed FSU and (c) L H C M . C.i.= 0.1°C. Fig. 3.9b was smoothed by a 5-point running mean in time. Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 67 Figure 3.10: Time-longitude diagrams of SST anomalies along the equator from (a) N H C M , (b) control run forced with observed FSU and (c) L H C M . C.i.= 1°C. Fig. 3.10b was smoothed by a 5-point running mean in time. Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 68 150E 180E 150W 120W 90W C 150E 180E 150W 120W 90W e 150E 180E 150W 120W 90W d 150E 180E 150W 120W 90W 30N 15N 0 15S 30S 150E 180E 150W 120W 90W f 150E 180E 150W 120W 90W Figure 3.11: Spatial patterns of the ENSO-related POP mode for L H C M . Patterns P (real) are shown on the left, and patterns —Q (imaginary) appear on the right. SST anomahes are in panels (a) and (b), HC anomahes in panels (c) and (d), and zonal wind stress in panels (e) and f. C.i.= 0.005 in all panels. Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 69 (Chang et al. 1995). Fig. 3.12 is the histogram of the number of warm ENSO events per month of the calendar year according to the time series of SST anomalies averaged over the NIN03 (150°W-90 oW,5 oN- 5°S) region for both coupled models. While annual phase-locking exists in both models, the phase-locking in L H C M appears to be too strong and about 3 months too late in comparison with observations, implying that the seasonal cycle is over-emphasized in L H C M . In contrast, the more scattered phase-locking be-haviour and the timing of the phase-locking found in N H C M (with warm events peaking mainly in late fall and winter, and never in summer) are more realistic ( Syu and Neelin 1998). The M E M (maximum entropy method) analysis of the NIN03 index derived from the N H C M simulation exhibits double significant spectral peaks at 53-month period and 28-month period, compared with only a single significant peak at 85-month period in the L H C M , revealing that the N H C M oscillations are more complicated than simply the sinusoidal oscillations found in L H C M . 3.4 Sensi t ivi ty experiments To fully explore nonlinear interactions between the seasonal cycle and interannual oscil-lations in the two coupled models, a series of experiments are conducted by varying the coupling parameter 8. Each experiment consists of a three-year initial spin-up followed by a 47-year simulation. The experiments show that the oscillatory behaviour of both L H C M and N H C M has a sensitive dependence on the parameter, with L H C M being slightly more sensitive. In the first experiment, we decrease 8 to 0.83, a value equivalent to setting the adjusting factor in section 3.2 to unity for the standard coupling runs. A damped oscillation with a period of 55 months was found in the N H C M , compared with a similar oscillation with a 75-month period in the L H C M . As the coupling strength 8 is gradually increased starting from the standard coupling, the oscillatory behaviour of the Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 70 2 1 0 -1 -2 -3 20 40 60 80 100 120 1 2 3 4 5 6 7 8 9 1011 12 1 2 3 4 5 6 7 8 9 1011 12 Figure 3.12: The NINO 3 index found in N H C M and L H C M in the standard case, (a) The time series of NINO 3 index, with solid curve for N H C M and dot-dash curve for L H C M . A histogram of the number of warm events peaking in Jan. (month 1) to Dec. (month 12) for (b) N H C M , and (c) L H C M . Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 71 two coupled models are basically similar to those obtained in the last section, i.e. reg-ularity and phase-locking. As the coupling parameter increases, the period is gradually shortened. As the coupling parameter 8 is further increased up to a specific critical value, the two models finally produce a 2 year quasi-periodic oscillation, with the peaks always occurring at the same months of the calendar year, as in many other ENSO models (Syu and Neelin 1995; Chang et al 1995, 1996). The critical value is 1.5 for L H C M and 1.68 for N H C M . The lower critical value for the L H C M results from the L H C M being more sensitive to 8 than the N H C M . Frequency locking has been noted in a number of other ENSO models (Anderson and McCreary 1985; Battisti 1988; Syu et al. 1995). Chang et al (1995, 1996) and Tziperman et al.(1994) investigated the relationship between such a frequency-locking behaviour and the transition to chaos in several ENSO coupled models. They found that a further increase of the coupling strength led to chaotic behaviour of the model ENSO. Interestingly, a further increase of the coupling parameter 8 in our models did not exhibit the typical behaviour of a chaotic system. Instead a further increase of the coupling parameter beyond the critical value leads to a climate drift for N H C M (8 = 2.0) and for L H C M (8 = 1.66). The climate drift at high coupling is likely associated with the spurious effects of flux correction in the coupled system as discussed in Neelin and Dijkstra (1995). Further analyses for the N H C M and L H C M simulations with the critical values are given in Fig. 3.13. A mode-locked state with a quasi-periodic cycle of 2 year can be found for both N H C M and L H C M in Fig. 3.13. The L H C M has a very strong phase-locking to January; but N H C M shows a more realistic broader phase locking. In the M E M analyses, the single sharp peak is clearly shown, and a broad spectrum, a necessary indicator of chaos, does not appear for either model. In the construction of the phase plane, we still use the methods commonly used in ENSO studies, i.e. decorrelated time-step K for Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 72 the time lag arid 2 dimensions SST(i) and SST(£ -f- K) for dynamical dimensions, to be consistent with the discussions of other ENSO models. Fig. 3.13e and Fig. 3.13f are the phase planes reconstructed with K—W months, which also clearly show the limit cycle solution for the both models. As we have not covered the entire parameter regime in our search, there might exist some specific parameters which can produce chaotic behaviour associated with ENSO in our models, but recent works do show that the random forcing may play a central role in ENSO evolution ( Blanke et al. 1997; Eckert and Latif 1997). The irregular interannual oscillation simulated by coupled CGCMs (Philander et al. 1992) does not appear to be low-order chaos. The nonlinear time series analysis performed by Chang et al. (1996) suggests that stochastic processes, rather than chaotic dynamics, are likely to be a major source of ENSO irregularity in CGCMs and in nature. In summary, the characteristics of the N H C M (with an N N atmosphere) and the L H C M (with an L R atmosphere) in terms of ENSO oscillations have been investigated in this chapter. A series of sensitivity experiments were carried out to further examine the dynamic behaviours of the two coupled models. The results show that the L H C M has a longer ENSO oscillation period (87 months), than the more realistic period found in the N H C M (57 months), and a very narrow phase-locking to the annual cycle compared to the more realistic phase-locking in the N H C M . Chapter 3. Hybrid coupled models of the tropical Pacific — Interannual variability 73 Figure 3.13: Analyses of the NIN03 index for N H C M and L H C M at the critical coupling parameter S. N H C M are shown on the left (5 = 1.68), and L H C M appear on the right (5 = 1.5). (a) and (b) are the histograms of the number of warm events peaking in month 1 to 12 (i.e. January to December); (c) and (d) are the power spectra; (e) and (f) are the phase planes reconstructed from the time series. Chapter 4 Hybrid coupled models of the tropical Pacific — ENSO prediction 4.1 Introduction In chapter 2 and chapter 3, We have developed two hybrid coupled models, and examined their ENSO behaviour. Although the dynamic response of the atmosphere and ocean in the tropical Pacific is quasi-linear (Chapter 2), the coupled interactions between the two components seem more complicated than the relations extracted from data and from uncoupled model responses (Chapter 3). In this chapter, we compare the N H C M and L H C M in ENSO prediction. This chapter is scheduled as follow: Section 4.2 briefly describes the prediction strat-egy. Section 4.3 displays the prediction skills with the 'stress-only' initiahzation scheme for the N H C M and L H C M . Section 4.4 shows the prediction experiments with the alter-native 'stress-from-ocean' initialization scheme for the two coupled models and, gives the prediction results of several typical ENSO episodes. 4.2 The strategy of prediction experiments For both N H C M and L H C M , two kinds of prediction experiments were carried out in this study, 'hindcasts' from 1964-1990 and 'forecasts' from 1990-1998. More precisely, during the 'forecast' period, the latest forecasts were made in 1998, but the forecasts were into 1999, while during the 'hindcast' period, the forecasts were into 1991. To aUeviate the artificial skills, the cross-validated atmospheric models described in Chapter 2 were used 74 Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 75 in the hindcasts experiments, i.e. the atmospheric models were developed excluding data from the hindcast period. For example, when the period 1970-1980 is being hindcasted, the atmospheric models used were only built from the data of 1964-1969 and 1981-1990. Thus, no data from the hindcast period have entered the coupled model through L R or N N reconstructed winds in any hindcast experiments. For the forecast experiments of 1990-1998, the atmospheric models are the ensemble of individual cross-validated models built for the hindcast experiments from 1964-1990, so data from 1990-1998 were never used in model building. There is a slight difference with the atmospheric models used in Chapter 3. The atmosphere models with upper ocean heat contents as predictors were used here for all the prediction experiments because they yield slightly better skills than the atmospheric models with SST as predictors (Chapter 2). The dynamical behavior and interannual variability of the coupled models with upper ocean heat contents as predictors are generally very similar to these described in Chapter 3. An important factor affecting the ENSO oscillation and model predictability is the coupling coefficient S (Chapter 3), as the change in the coupling strength can easily modify the amplitudes of the coupled oscillation, and exert a strong influence on the frequency locking behavior (Chang et al. 1996; Chapter 3). After many sensitivity experiments, we found that the best hindcast and forecast skills for L H C M and N H C M are attained neither from the standard coupling with regular ENSO oscillations as in Chapter 3, nor from stronger coupling than the standard coupling. The best hindcast and forecast skills occurred for 'damped' coupling, with S between 0.7 to 0.8. We thus chose 8 to be 0.7 for all hindcast and forecast experiments presented in this study for both N H C M and L H C M . The seasonal cycle, used to obtain the surface sea temperature anomalies (SSTA) for evaluating the model skills, was calculated from the hindcast experiments during 1964-1990. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 76 4.3 Prediction experiments with 'stress-only' initializations The ocean model was forced with the observed FSU (Florida State University) wind stress (Goldenberg and O'Brien 1981) from 1961 to 1998. From April , 1964 onward, at three months intervals (1 January, 1 April , 1 July, 1 October), the conditions of this control run were selected as the initial conditions for the N H C M and L H C M , which were run for 15 months. We will refer to this type of initialization as 'stress-only' initialization, to be contrasted with a different initialization in the next section. For each H C M , a total of 104 runs were made starting from April 1964 to January 1990 (the hindcast period), and 35 runs starting from April 1990 to October 1998 (the forecast period). Fig. 4.1 shows the N H C M and L H C M predictive skills for the NIN03 index for the hindcast period and the forecast period. In the first 2 months, the correlation skills of the two HCMs are below persistence (the predictions at all lead times being simply the initial observed SSTA), due to initialization shock caused by the initial conditions taken from the control run not being consistent with the governing equations of the coupled model. The two HCMs outperform persistence at lead times greater than 3 months, with the skills of the HCMs decreasing gradually up to a lead time of 6 months, beyond which the skills are relatively stable up to a lead time of 15 months. For N H C M , the correlation is around 0.5 up to a lead time of 14 months for the hindcast period and greater than 0.5 at all lead times up to 15 months for the forecast period (Fig. 4.1). Compared to N H C M , L H C M has the slightly higher skills during the hindcast period but slightly lower skills during the forecast period. Both coupled model predictions have weak amplitudes, in particular at long lead time. In order to alleviate the weak amplitude, a systematic error correction is adopted for the predictions at each lead time by scaling: S S T A S S T A ™ = ^±L<Tobs, (4.1) CT Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 77 a b Lead time in months Lead time in months Figure 4.1: Correlation skills for the predictions of SST anomahes in the NIN03 area as a function of the forecast lead time for (a) the hindcast period from 1964-1990, and (b) the forecast period from 1990-1998, by N H C M (solid), L H C M (dashed) and persistence (dot-dashed). Both HCMs use the 'stress-only' initiahzation scheme. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 78 where r is the lead time, and SSTA is model prediction and S S T A n e u ; is the corrected prediction. crT, the standard deviation of SSTA at lead time r , is obtained from the 104 hindcast cases. The 0~obs is the standard deviation of the observed SSTA from the hindcast period of 1964-1990. The systematic error correction does not change any correlation scores but effectively improves the R M S E skills for the prediction of ENSO events. Fig. 4.2 and Fig. 4.3 show the time series of predicted NIN03 SSTA by the N H C M at lead times of 3 , 6 and 9 months (solid line). The predicted NIN03 SSTA shown here has removed the systematic errors using Eq. 4.1. The observed NIN03 SSTA are also shown (dash line) for comparison. The N H C M predicted the low-frequency interannual evolution ofthe SST anomalies reasonably well, especially at short lead time. The model successfully predicted all typical E l Nino and La Nino events, such as the 1972/73, 1982/83, 1986/87, 1997 E l Ninos and the 1988, 1998 La Nihos, up to a lead time of 9 months. For lead time greater than 9 months, the model has also useful predictive skills and can to some extents predict these typical ENSO events in amphtude and phase (not shown). But for the two strongest warm events 82/83 and 97, the model suffered the same problem as many prediction models, i.e. weak amphtude prediction for 82/83 and lagging behind in the phase prediction for 97 (e.g. Syu and Neehn 2000; Fliigel and Chang 1998; Saunders et al. 1998). The model also produced a spurious warm event for 74/75 (a good case to test model skill in ENSO community) at the prediction of 9 month lead time, like many coupled models (Wu et al. 1994). In the observation, there seemed to be an incipient E l Nino in the middle of 1974, but it aborted later in 1974. However the model predicted the aborted event well at short lead time. The time series of predicted NIN03 SSTA with L H C M shows almost the same features as N H C M . Fig. 4.4 is its prediction for the forecast period. Compared with the N H C M , the amphtude prediction in L H C M seemed weaker than that in N H C M . Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 79 Figure 4.2: Predicted NIN03 SST anomalies (°C) by the N H C M from 1964-1990 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for observations. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 80 Figure 4.3: Predicted NIN03 SST anomalies (°C) by the N H C M from 1990-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Sohd line is for the predicted values and dash line for observations. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 81 Figure 4.4: Predicted NIN03 SST anomalies (°C) by the L H C M from 1990-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for observations. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 82 Decadal dependence of ENSO predictability has been found in many ENSO forecast models (e.g. Balmaseda et al. 1995; Chen et al. 1997; Flugel and Chang 1998). While all models tended to have very good forecast skiUs in the 1980s, they suffered low skiUs in the 1990s— even with an improved initiahzation strategy (Chen et al. 1997). A clear contrast in terms of the characteristics of the interannual variability between 1980s and 1990s also exists in many observations such as sea level pressure, SST, low-level zonal wind, and subsurface ocean heat content anomahes in the Pacific (Kleeman et al. 1996; Latif et al. 1997; J i . et al. 1996). It has been suggested that the mechanism of decadal dependence of predictability is associated with the decadal changes in the mean state leading to decadal variability of ENSO (e.g., Wang 1995; Zhang et al. 1997). Several recent studies suggested that the changes in mean state is probably caused by (1) the remote response in the tropical atmosphere to the mid-latitude decadal osciUations (Wallace et al. 1998); (2) anthropogenic global warming (Trenberth and Hoar 1997); (3) the interaction between tropical and extratropical oceans to produce decadal variations in the mean thermocline structure of the tropical Pacific (Gu and Philander 1997); and (4) amplification of the decadal signal through atmospheric high-frequency variability (Kirtman and Schopf 1998). Fig. 4.5 shows the decadal dependence of predictabihty in the two coupled models. The obvious decadal changes of predictability appear only for L H C M . For N H C M , except for low skill in the 1970s which is probably associated with poor data quahty and low signal-to-noise ratio (Chen at al. 1997), prediction skills did not have significant changes among the 1960s, 1980s and 1990s. In particular, in L H C M , there is a large drop in skills from the 1980s to the 1990s, which is absent in N H C M . What could be the cause for this difference? Using a neural network model for nonhnear canonical correlation analysis ( N L C C A , Hsieh 2000), Hsieh (2001) found that the relation between the tropical Pacific SST and Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 83 Figure 4.5: Predictive (correlation) skills of the NIN03 SSTA from (a) L H C M and (b) N H C M , for the 1990s (solid), 1980s (circled), 1970s (dashed), and the 1960s (dot-dashed), as a function of the lead time. Both coupled models use the 'stress-only' initialization scheme. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 84 the sea level pressure has become more nonlinear with time. We apphed the N L C C A of Hsieh (2001) to study the relation between the first 3 principal components (PCs) (i.e. the E O F time coefficients) of the observed zonal wind stress anomahes and the first 3 PCs of the observed SSTA. P C analysis (i.e. E O F analysis) was first performed for the period 1964-1999 for the zonal wind stress anomahes and the SSTA separately. The N L C C A was then applied to the PCs in the 1980s, then to the PCs in the 1990s. The first N L C C A mode for the 1980s is shown in the PC-space of the zonal wind stress anomalies (Fig. 4.6), and in the PC-space of the SSTA (Fig. 4.7). The classical C C A solutions are shown as straight hnes in Figs. 4.6 and 4.7 for comparison. As a measure of the degree of nonhnearity in the relation, the ratio R between the M S E (mean square error) of the N L C C A mode to the M S E of the C C A mode can be calculated. A small R would indicate a nonhnear relation, while R 1 would indicate an essentially hnear relation. R is 0.941 for the zonal wind stress, and 0.917 for the SSTA in the 1980s. In contrast, the N L C C A mode 1 is shown in Figs. 4.8 and 4.9 for the 1990s. It is evident the relation in the 1990s has become more nonhnear. R is 0.810 for the zonal wind stress and 0.507 for the SSTA, both considerably lower than the corresponding values for the 1980s. Using the nonlinear technique of N L C C A , we have found that indeed the relation between the zonal wind stress and the SSTA has become more nonlinear in the 1990s than in the 1980s. This could explain why L H C M , which did well in the 1980s, fared poorly in the 1990s, and why N H C M , which is capable of simulating the nonhnear relation between the wind stress and the surface ocean conditions, did not suffer a loss of forecast skills as the relation turned more nonhnear for the 1990s. These decadal comparisons suffer from the fact that a decade is far too short for statisical significance. One can only interpret the N L C C A findings as being at least consistent with our finding that the L H C M forecast skills appeared lower than those of the N H C M during the 1990s. As in other models (e.g. Balmaseda et al. 1995; Zebiak and Cane 1987), the two Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 85 CM O CL 0.1 0.05 0 -0.05 -0.1 (a) -0.15 -0.2 -0.1 0 PC1 0.1 0.2 -0.21 -0.15 -0.1 -0.05 0 PC2 0.05 0.1 Figure 4.6: Nonlinear C C A (NLCCA) mode 1 for the 1980s shown in the P C (principal component) space of the zonal wind stress anomalies. The N L C C A analyzed the relation between the first three PCs of the zonal wind stress anomalies and the first three PCs of the SSTA. The data are shown as dots, and their projection onto the N L C C A mode 1 are the circles, and the classical C C A solution is shown as a straight line for comparison. The solutions are shown in (a) the PC1-PC2 space, (b) the PC1-PC3 space, (c) the PC2-PC3 space, and (d) a 3-D view. Following the principle of parsimony, the architecture of the N L C C A was set for minimal nonlinearity, i.e. only 2 hidden neurons used, and the weight penalty parameter was 0.1 in all cost functions (Hsieh 2001). Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction Figure 4.7: N L C C A mode 1 for the 1980s in the P C space of the SSTA. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 87 Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 88 Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 89 0.9 0.3 0.2 0.9 1 4 7 10 13 Lead time in months 0.2 1 4 7 10 13 Lead time in months Figure 4.10: Seasonal dependence of the predictive skills for the NIN03 SSTA during 1964-1998 from (a) L H C M and (b) N H C M . Both coupled models use the stress-only initialization scheme initialized on January 1 (dot-dashed), April 1 (solid), July 1 (dotted) and October 1 (dashed). coupled models show the seasonal dependence of predictability as shown in Fig. 4.10, where the correlation is calculated with all 139 cases during both the hindcast and forecast periods. Clearly, the 'spring barrier' appears for both models, i.e. the correlation skill of the predictions decreases markedly as the predictions pass through northern spring. The two models share a similar characteristic of seasonal dependence of predict ability. The cause of seasonal dependence of predictability is still not clear. One explanation is that spring has the minimum signal-to-noise ratio, i.e. the signal variance is minimum (Xue et al. 1994). However this is only true for SSTA. For other observed quantities such as wind stress and heat content, their variances are not minimum in spring (Fig. 4.11). Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 90 Figure 4.11: Seasonal variation of variances of NIN03 index for (a) observed SST, (b) observed HC , (c) observed zonal (solid) and meridional (dashed) wind stress and (d) model SST taken from control run forcing with FSU winds. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 91 4.4 Prediction experiments with 'stress-from-ocean' initialization The above prediction experiments are all based on the initialization scheme which only uses the surface wind stress information and ignores the model feedback from the ocean to the atmosphere. Although it has shown considerable success in establishing initial conditions for our models, we will further explore a different initialization scheme, which takes into account the feedback from the ocean to the atmosphere, by using ocean states, including subsurface ocean states. This 'stress-from-ocean' initialization scheme uses the output from the ocean model's control run forced by the observed wind stress to reconstruct the atmosphere wind stress which is then provided to initialize the ocean model for predictions. It thus involves two steps: (1) the ocean model is forced with the observed F S U wind stress to obtain the model ocean states; (2) the upper ocean heat contents from (1) are input into the empirical atmospheric models (LR or NN) to obtain the model wind stress. This model wind stress then drives the ocean model, until the start of a forecast when the H C M is run for 15 months. As in section 4.3, 104 hindcast cases and 35 forecast cases were carried out with initialization on January 1, April 1, July 1 and October 1. As upper ocean heat contents are used to derive the wind stress and only the first 3 E O F modes were used in the atmospheric models for estimating the zonal or meridional wind stress, this initialization scheme largely removes the weather noise in the initial conditions. Such a 'stress-from-ocean' initialization has been used in Barnett et al. (1993) and Syu and Neelin (2000). The correlation skills of the predictions by the stress-only initialization scheme is compared with those of the stress-from-ocean initialization scheme, during the hindcast period from 1964-1990 (Fig. 4.12) and the forecast from 1990-1998 (Fig. 4.13) for L H C M and N H C M . In Fig. 4.12, the prediction skills of the stress-from-ocean initializations are Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 92 almost identical to those of the stress-only initializations for L H C M during the hindcast period. However, the stress-from-ocean initiahzations have higher prediction skills at lead times greater than 6 months for N H C M during the same period. For the forecast period of 1990-1998, the skill improvement by the stress-from-ocean initializations can be observed for both models, slightly more so for L H C M (Fig. 4.13). Thus it appears that the stress-from-ocean initialization is better than the stress-only initiahzation. The time series of predicted NIN03 SSTA under the stress-from-ocean initiahzation for lead times of 3, 6 and 9 months (sohd line) are shown in Fig. 4.14 for N H C M and in Fig. 4.15 for L H C M . As the standard deviation of the predicted values are generally smaller than the observed ones, especially at longer lead times, the predicted values have been scaled to match the observed standard deviations in Figs. 4.14 and 4.15. A comparison of predicted skills between L H C M and N H C M (Fig. 4.16) reveals that the skills are very similar during the hindcast period, but N H C M has higher skills than L H C M during the forecast period. This is again consistent with the earlier finding by N L C C A that the relation between zonal wind stress and surface ocean condition has become more nonlinear in the 1990s, which would put L H C M in a disadvantage against N H C M . 4.4.1 Case study In previous sections, we have presented the general skills of N H C M and L H C M with two different initialization schemes. Here we will examine the model predictive capability for typical ENSO events by case studies. Three most renowned episodes of ENSO, the 1982/83 and 1997 E l Nino and 1998 La Nina, were chosen for study. As the stress-from-ocean initialization generates slightly better skill than the stress-only initialization for both N H C M and L H C M , the cases were thus selected from the stress-from-ocean initialization experiments. Each case starts from the 1 January initial conditions and the Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 93 Figure 4.12: Contrast of the predictive skills for the NIN03 SSTA during 1964-1990 between the stress-only initialization scheme (sohd hne) and the stress-from-ocean ini-tiahzation scheme (dashed line) for (a) L H C M and (b) N H C M . Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 94 Figure 4.13: Contrast of the predictive skills for the NIN03 SSTA during 1990-1998 between the stress-only initialization scheme (solid line) and the stress-from-ocean ini-tialization scheme (dashed line) for (a) L H C M and (b) N H C M . Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 95 Figure 4.14: Predicted NIN03 SST anomalies (°C) by N H C M using the stress-from-ocean initialization scheme during 1964-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for the observations. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 96 Figure 4.15: Predicted NIN03 SST anomalies (°C) by L H C M using the stress-from-ocean initiahzation scheme during 1964-1998 at (a) 3-month, (b) 6-month and (c) 9-month lead times. Solid line is for the predicted values and dash line for the observations. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 97 Figure 4.16: Correlation skills of the predicted NIN03 SSTA for (a) the hindcasts period from 1964-1990, and (b) the forecast period from 1990-1998, by N H C M (sol id) , L H C M (dashed). Both H C M models used the stress-from-ocean initialization scheme. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 98 forecast period is 15 months. The predicted value at each grid point is scaled to match the observed standard deviation at that grid. The evolution of the predicted monthly SSTA along the equator are shown in time-longitude sections. Fig. 4.17 shows the 1982/83 E l Nifio. Compared with the observed SSTA, N H C M and L H C M both successfuUy predicted the peak of this event about a year ahead. A restart of warming at March 1983 (i.e. month 15) in the observations were also well captured by the two models about 15 months in advance. But both models predicted the start of the warming about two months too early. The predicted warming was too narrowly confined around 120°W, and the equatorial western Pacific also had spurious cooling predicted. For the 1997 E l Nino (Fig. 4.18), the two coupled models predicted weU the start time, mature time, and the amplitude for this strongest event in this century. But they both have later decay and end time than the observations. The predicted warmings were too narrowly confined around 120°W, and entirely missing near the eastern boundary, when compared with observations. Again, the predictions by both N H C M and L H C M had spurious cooling in the west. In the 1998 La Nino (Fig. 4.19), N H C M displays a much better prediction than L H C M . This event appeared at June 1998 and reached its mature stage at around December 1998. The maximum anomahes reach -2°C. These features were basically predicted by N H C M except the start time, which was 2 months late. In contrast, L H C M predicted the start time 5-6 months later than observed. In the equatorial Pacific just east of the dateline, the models performed badly, missing the the transition from warm to cold conditions, especially for L H C M . In summary, the prediction skills for the two coupled models, N H C M and L H C M , with two different initialization schemes have been studied in this chapter. The ensemble of 104 'hindcasts' experiments from 1964-1990 and 35 'forecasts' experiments were performed Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 99 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Longitude Longitude Longitude Figure 4.17: Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1982. The observed SSTA from January 1982 to March 1983 is shown in the middle panel. The ordinate ranging from 1 to 15 is the calendar month from January 1982 to March 1983. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 100 1 L : u J I. \ > '. J 1 l i : i / i1 . i I 1 | _ J j i_, L L , J 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W 150E 180E 150W 120W 90W Longitude Longitude Longitude Figure 4.18: Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1997. The evolution of the observed SSTA from January 1997 to March 1998 is shown in the middle panel. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 101 Figure 4.19: Time-longitude diagrams of predicted SSTA along the equator from N H C M (left panel) and L H C M (right panel) initialized on 1 January 1998. The evolution of observed SSTA from January 1998 to March 1999 is shown in the middle panel. Chapter 4. Hybrid coupled models of the tropical Pacific — ENSO prediction 102 to examine the prediction skills of the two coupled models. Three typical ENSO episodes were chosen to further show the predictive ability of the two coupled models. The stress-from-ocean initiahzation scheme generally has better predictive skills than the stress-only scheme. The main difference in the forecast skills between the N H C M and the L H C M occurs in the 1990s, where the N H C M has better skills. The N L C C A was used to further examine the possible reasons to cause the prediction skill differences between the two coupled models in the 1990s. Chapter 5 Impact of data assimilation on ENSO simulations and predictions 5.1 Introduction In Chapter 4, we examined ENSO predictability in two hybrid coupled models. The results show that the initialization has an important impact on model skills. ENSO prediction is an initial value problem, and the further evolution of the system depends highly on the initial state from which it started. This is at least valid over a certain period, which is called the predictability time. Although ENSO is a phenomenon of the ocean-atmosphere coupled system, its evolution is determined mainly by the ocean. Therefore, the initialization can be considered as an initial value problem for the ocean. Thus, a very important task in carrying out an ENSO forecast is to determine the oceanic initial state as accurately as possible. Knowing that neither the model nor the observations are perfect, the goal of data assimilation is to estimate optimally the past and present states with an imperfect model and noise observations. As discussed in Chapter 1, techniques for data assimilation and initialization in the ENSO field are still at a relatively early stage of development. We are at present unclear on a basic issue, i.e. which observations are most useful in initialising ENSO predictions and how do the prediction skills differ for assimilating different types of data? It seems that we need a systematic comparison of assimilating different types of data on ENSO predictability based on an identical assimilation system and an identical forecast model. Another key issue is how to assimilate subsurface heat content data into an intermediate 103 Chapter 5. Impact of data assimilation on ENSO simulations and predictions 104 complexity layer model which has low vertical resolution. Compared with the model prognostic variables such as SST, the distribution of upper ocean heat content anomahes dominate ENSO characteristics and predictabihty. As an integral of the thermal infor-mation in upper ocean, the upper ocean heat content anomahes extracted from data are less sensitive to missing values in datasets than the model prognostic variables. The problem of missing values often occurs in subsurface data sets (e.g. TAO and X B T ) . A n ideal method to assimilate heat content data into a numerical model is to use the adjoint 4-D Var approach, which can define a cost function involving the heat content. The heat content correction can be transfered to the corrections of model prognostic variables directly. However, expensive computation cost and complicated coding hmit the application of the 4-D Var approach. Therefore, our main focus in this chapter is on: 1) how to assimilate observations, especiaUy subsurface heat content data into a numerical model with a simple 3D Var assimilation scheme; and 2) which data type is the most useful in initializing the ocean model to improve the predictive skiUs? This chapter is structured as follows: Section 5.2 briefly describes the assimilation system. Section 5.3 introduces the data and proposes strategies for assimilating sea surface height anomalies and upper ocean heat content into the coupled model. The impact of assimilating different types of data on ENSO simulation and prediction are presented in Sections 5.4 and 5.5, respectively. 5.2 The assimilation system The data assimilation system is the 3D Var assimilation system described by Derber and Rosati (1989). It has been used in operational ENSO forecasts at N C E P (National Center of Environment Prediction) and at G F D L (Geophysical Fluid Dynamics Laboratory) (Ji Chapter 5. Impact of data assimilation on ENSO simulations and predictions 105 et al. 1998; Rosati et al. 1997). The 3D Var minimizes the cost function / = ^ T T E - 1 T + ^ ( D ( T ) - T 0 ) T F - 1 ( D ( T ) - T o ) (5.1) The superscript T represents the transpose, T is an N component vector containing the correction to the guess field (the guess field being generated by the model before assimilating the latest data), E is an estimate of the N x N first guess error covari-ance matrix, D is a simple bilinear operator interpolating from the model grid to the observation stations, T 0 is an M component vector containing the difference between the observations and the interpolated guess field, and F is an estimate of the M X M observational error covariance matrix, where N and M denote the number of model grid points and the number of the observational stations respectively. The first term of the functional ( T T E _ 1 T ) is a measure of the fit of the corrected field to the guess field, while the second term measures the fit of the corrected field to the observations. The result is a weighted average of the guess field (which contains information from earlier periods) and the observations. The estimated spatial error covariance matrices, E and F determines the spatial structure and amplitude of the correlation field. An ideal method to estimate the matrices is to use the Kalman filter (Gelb 1974). However, extremely high computation cost for the Kalman filter prohibits its use on coupled models. Here a simple Gaussian function is used to determine E while ignoring the vertical correlations as in Derber and Rosati (1989), i.e. the horizontal covariance between any two points is given by where (f> is the latitude of the grid point, r the distance between any two points, o=0.01, and b = 570 km, as in Derber and Rosati (1989). For the observational error covariance matrix F, we assume that the observational errors are not correlated, i.e, the matrix Chapter 5. Impact of data assimilation on ENSO simulations and predictions 106 F has only diagonal elements which are the observational error variances. For SST, surface wind stress and SSHA, each element of the diagonal error covariance matrix is constant, with the value (0.5°C) 2 for SST, (0.25 N m" 2 ) 2 for wind stress and (15 cm) 2 for SSHA. For the HC anomaly, the error variances assigned are (0.25°C) 2 in the western Pacific and (0.5°C) 2 in the eastern Pacific due to fewer observational stations in the east. These values are obtained empirically from sensitivity experiments, and following recommendations from previous studies (e.g. J i at al, 1995; Hao and Ghil, 1994). The function / in (5.1) is minimized using a preconditioned conjugate gradient al-gorithm (GiU et al. 1981). The preconditioning in the algorithm is supphed by the E matrix, which allows the solution to be found without directly inverting the E matrix. Further details about the algorithm can be found in Derber and Rosati (1989) and GiU et al. (1981). The assimilation system also includes data quahty control. The initial gross check is to ensure that the data used was not over land. This is implemented by one land-ocean mask in which the grid weights are set to 1 over ocean and 0 over land. Before the guess field and observations enter the assimilation system, they were both mask-filtered. Further quality control is to remove observations with large errors. When the magnitude difference between the guess field and the observations is greater than the 25% of the magnitude of the guess field, the observations are ehminated. This might lose effective corrections for poor simulations in a few model grid points, but ensures the ocean model is graduaUy adjusted, without problems arising from the imbalance between the velocity and density fields. The data are inserted daily into the model. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 107 Table 5.1: Information on the datasets used Data SST SSHA Wind stress H C A Domain in space 1°E-359°E, 122.25°E-71.25°W, 120°E-70°W 30°W-30°W, 89°S-89°N 35°S-45°N 30°S-30°N 60°S-60°N Resolution in space 2° x 2° 1.5° x 1° 2° x 2° 5° x 2° Domain in time 1950-2000 1980-1999 1961-1999 1955-1999 Resolution in time monthly monthly monthly monthly 5.3 Methodology for assimilating different types of data 5.3.1 Data The data used in this study are the monthly mean SST obtained from COADS , the monthly mean surface wind stress from FSU, the monthly mean sea surface height anomaly (SSHA) from the reanalysis data set of N C E P and the monthly 400-m depth-averaged heat content anomahes (HCA) from the Joint Environmental Data Analysis Center at Scripps Institution. Table 5.1 describes the datasets in some detail. SST and wind stress can be directly assimilated into the model because of their very realistic climatology and because they are prognostic variables in the model. With SSHA and H C A , the situation is more complicated. 5.3.2 SSHA observations As SSH is only a diagnostic variable in the model and so the SSHA observations cannot be directly assimilated into the ocean model. Fischer et al. (1997) used a statistical relation between the SSHA and the leading P C of the temperature vertical profile to inject the SSHA information into the temperature profile of their ocean model. We will use the SSHA to adjust the model interface (H1A) between the top layer and the second layer. The interface anomahes are mainly dominated by the lowest-order baroclinic mode Chapter 5. Impact of data assimilation on ENSO simulations and predictions 108 and approximately reflect the thermocline perturbations. A simple statistical relation between model SSHA and H1A can be found with E O F analysis. Fig. 5.1 compares the PCs of first two E O F modes for the observed SSHA, and the SSHA and the H1A from the control model run. The model SSHA was diagnosed by computing the sea surface pressure in the control run. The first two E O F modes accounted for a total of 67%, 64% and 61% of the variance for the observed SSHA, the model SSHA, and the model H1A, respectively. There is close agreement among these time series (Fig. 5.1). Fig. 5.2 shows the first two E O F spatial modes of the model SSHA and the H1A, also showing a good agreement between the spatial patterns. The almost opposite spatial structure between Figs. 5.2a and 5.2c and between Figs. 5.2b and 5.2d actually reflect the out-of-phase variations between the sea surface height and the thermocline due to baroclinic motion. So, SSHA can approximately be assumed to be negatively proportional to the H1A, i.e. ( S S H A ) i = - ^ ( H l A ) , , (5.3) where i indicates the zth observational station, and the proportional coefficient LI is simply taken to be the ratio of the standard deviation of the SSHA to that of the H1A at the ith station. With (5.3), the observed SSHA can be converted to HI A , a prognostic model variable. Gill (1982) dynamically derived the same relation as (5.3) for a 2-layer shallow-water model without external forcing, where \i is determined by the model vertical density profile. 5.3.3 HCA observations A n important component of the coupled ocean-atmosphere system on the interannual timescale is the upper ocean heat content, which is the source of memory for the system. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 109 Figure 5.1: The time series of (a) the first principal component (PC) (i.e. E O F time coefficient) and (b) the second P C , for the observed SSHA (sohd hne), the SSHA (dotted hne) from the control run (with the ocean model forced by the F S U wind stress), and the first interface H1A (dash line) from the control run. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 110 Figure 5.2: The spatial patterns of the first two E O F modes from the control model run for the SSHA (a and b) and for the HI A (c and d). Solid hnes indicate positive contours, dashed hnes, negative contours, and dotted hnes, the zero contour. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 111 Many studies have shown (e.g. Zebiak and Cane 1987; Rosati et al. 1997) that the fluc-tuations in upper-ocean heat content are both systematic and significant in the evolution of ENSO. According to (2.4), the model H C A can be computed at each model step; how-ever, assimilating the H C A data is not easy as H C A is not a prognostic model variable. When Kleeman et al. (1995) and Moore and Anderson (1989) assimilated subsurface temperature into their ocean model, they used a linear empirical relation to transfer the corrections in the subsurface thermal fields to the corrections of the model prognostic variables. We will use nonlinear NNs to find the empirical relations between H C A and the four prognostic variables— SST, H1A, the second layer temperature anomalies (T2A) and the second layer thickness anomalies (H2A). The procedure to model these relations is very similar to the construction of atmo-spheric models in Chapter 2: (1) An E O F analysis is applied to the control model H C A and each of the four prognostic variables during the period 1964-1990. The first three modes are retained to construct the empirical equation as they explained over 60% of the variance for all variables. (2) An ensemble of neural networks (25 NNs used here) is trained where the inputs to each N N are the 3 PCs of the H C A and the output is one of the PCs of a prognostic variable. Three hidden neurons are used in each N N . The final relation is the ensemble average over 25 NNs. The cross-validation scheme described in Chapter 2.3.2 is used to obtain each empirical equation of the 25 NNs to ensure that no training data are used for testing the prediction skills. A total of 12 ensembles of NNs are needed to reconstruct the 3 leading PCs of the 4 prognostic variables from H C A . (3) As the observed heat content is the integral of the temperatures over the upper 400 meters, while the model HC is defined as the integral of the temperatures over the first two layers (275 meters), the observed H C A values are scaled to match the standard deviation of the model H C A before being assimilated. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 112 0.1 0.05 0 CM o Q. -0.05 < 1-CO -0.1 CO -0.15 -0.2 -0.25 -0.2 0.15 0.1 0.05 O 0 Q. < £ -0.05 -0.1 -0.15 0 0.1 HCA PC1 -0.2 -0.2 -0.1 0 HCA PC1 0.1 0.2 Figure 5.3: Scatter plot of the first P C of the H C A versus (a) the second P C of the SSTA, and (b) the first PC of the H2A, the thickness anomahes of the second layer. The observations are shown as dots, the N N relation is indicated by the small circles and the linear regression by the straight line. Fig. 5.3 depicts some of the empirical relations between the first H C A P C and the PCs of the prognostic variables, showing that the nonhnear NNs better represent the relations than hnear regression. The sum of squared error from the N N simulations are about 15% and 6% lower than those from the linear simulation in Fig. 5.3a and b respectively. Next, we perform a consistency check: H C A from the control run was used to estimate the prognostic variables, which were then used to recompute H C A using (2.4). Fig. 5.4 shows a time-longitude plot of the H C A along equator during 1980-1998, taken from observations, control run and recomputed from the prognostic variables. As shown in Fig. 5.4c, the recomputed H C A agreed well with the original control run values, and both agreeing fairly weU with the observed values. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 113 8 4 8 2 0 -ay o 1 2 0 E 1 6 0 E 1 6 0 W 1 2 0 W 8 0 W L o n g i t u d e 9 6 9 4 92 CD E 9 0 8 8 8 6 8 4 82 r 1 6 0 E 1 6 0 W 1 2 0 W 8 0 W L o n g i t u d e 9 6 9 4 92 CD E 9 0 8 8 8 6 8 4 82 1 6 0 E 1 6 0 W 1 2 0 W 8 0 W L o n g i t u d e Figure 5.4: Time-longitude plots of the ocean heat content anomalies along the equator during 1980-1998 from (a) the observations and (b) the control run. (c) The H C A calculated from the prognostic variables estimated by the cross-validated ensemble N N model using (b) as predictors. Values prior to 1980 are not shown, as the quality of the observed HC is poor. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 114 5.3.4 Assimilation strategies Because two different classes of data, i.e. oceanic data and wind stress data, are as-similated into the numerical model, two different strategies need to be applied. For the oceanic data (SST, SSHA and H C A ) , the assimilation is performed only for ocean model. As the SSHA data starts from 1980 only and the H C A data has relatively poor quahty prior to 1980, aU assimilation experiments are therefore set to start from 1980. As in the control run, the ocean model is forced with observed F S U wind stress, but monthly oceanic data are assimilated during the 'hindcast' period, January 1980-December 1989, and the 'forecast' period, January 1990-October 1998. Predictions are then made by running the coupled model forward in time from the established ocean initial conditions for 15 months. For the wind stress data, the assimilation is carried out for the hybrid coupled model. The observed FSU wind stress are assimilated into the N H C M for 6 months before letting the N H C M to proceed with the 15-month prediction run. The assimilation domain is confined to the equatorial belt, 15°S - 15°N. Decadal dependence of ENSO predictability has been found in many ENSO forecast models (e.g. Balmaseda et al. 1995; Chen et al. 1997; Fhigel and Chang 1998). While aU models generally have very good forecast skiUs in the 1980s, they tended to have lower skills in the 1990s— even with an improved initialization strategy (Chen et al. 1997). This might be caused by decadal changes in the mean state that could lead to the decadal variability of ENSO (e.g., Wang 1995; Zhang et al. 1997). We will study the skill differences between the 1980s and the 1990s. 5.4 Impact of data assimilation on ENSO simulation First, we will examine the impact of data assimilation on ENSO simulation in the 1980s and 1990s. The ocean data was continuously assimilated into the ocean model during Chapter 5. Impact of data assimilation on ENSO simulations and predictions 115 1980-1989 and 1990-1998 so that a continuous assimilated ocean field from 1980-1989 and 1990-1998 can be obtained. In contrast, for the wind stress data, the assimilation was carried out by repeatedly integrating the coupled model for 6 months prior to each forecast, hence there is not a continuous assimilated field. So our discussions are restricted to ocean data assimilation in this section. As the assimilated ocean fields are also the initial fields of the predictions to be discussed in the next subsection, the following analyses will actually explore the features of the initial fields for the predictions. The correlations between the observed tropical Pacific SSTA and the modeled SSTA with or without data assimilation are shown in Fig. 5.5 and Fig. 5.6 for the 1980s and 1990s, respectively. For the 1980s, the simulations with data assimilations are generally better than that without data assimilation (Fig. 5.5a). The most pronounced improve-ment in the tropic Pacific SSTA simulation is from SST assimilation (Fig. 5.5b), with high correlation (over 0.8) between observations and simulations for the whole assimi-lation domain. Since the model simulations are compared to the same data that were assimilated, it is not surprising that such a high correlation was achieved. For H C A assimilation, the improvement is less than with SST assimilation but is stiU considerable (Fig. 5.5d). The assimilation of SSHA brought httle improvement to the simulation. The correlation skills from the 1990s (Fig. 5.6) share similar features to those in the 1980s, i.e. SST assimilation has the best skills, followed by H C A assimilation, while SSHA assimilation contributed almost nothing to improve the simulation. Comparing the 1990s (Fig. 5.6) with the 1980s (Fig. 5.5), one notes that the cor-relation skills of SSTA simulation dropped in the 1990s for all cases with and without data assimilation. The skill improvements from data assimilation in 1990s is less obvious than in the 1980s. For example, the 0.8 correlation skills with SST assimilation are only confined to 10°S-5°N in the 1990s, in contrast to the 1980s, where the 0.8 correlation can be found over the whole assimilation domain (15°S-15°N). In the 1980s, the H C A Chapter 5. Impact of data assimilation on ENSO simulations and predictions 116 Figure 5.5: Correlations between the observed tropical Pacific SSTA and the SSTA from the ocean model forced by the FSU wind during 1980-1989, (a) without data assimilation, (b) with SSTA assimilation, (c) with SSHA assimilation, and (d) with H C A assimilation. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 117 Figure 5.6: Correlations between the observed tropical Pacific SSTA and the SSTA from the ocean model forced by the F S U wind during 1990-1998, (a) without data assimilation, (b) with SST assimilation, (c) with SSHA assimilation, and (d) with H C A assimilation. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 118 assimilation seems to have a remote impact, as correlation up to 0.6 can be found outside the assimilation belt (Fig. 5.5d). In contrast, such remote impact cannot be found in the 1990s for all ocean data assimilation. Fig. 5.7 shows the evolution of the observed and modeled SSTA in the Ni£io3 region. The correlations between the modeled and observed SSTA in the 1980s are 0.73 for the run without data assimilation, 0.93 for SST assimilation, 0.74 for SSHA assimilation and 0.79 for H C A assimilation, in contrast with the values of 0.76, 0.91, 0.76 and 0.77 respectively, in the 1990s. While all assimilation runs captured the major observed ENSO events, the H C A assimilation has the better capability for generating the anomalous extremes than SSHA assimilation and the run without assimilation although it sometimes exaggerates the extrema. The common deficiencies for all models (except that with SST assimilation) are the delay in simulating the occurrence of the 1982-83 warm event (Balmadesa et al. 1994), and the relatively weak simulated amplitude for the 1997 event. But the phase delay in 1982 and the weak amplitude in 1997 were much alleviated with H C A assimilation. In summary, assimilating several types of data can all benefit the ENSO (Niho3 SSTA) simulations in the 1980s, with the best simulation skills appearing in SST as-similation, followed by H C A and then SSHA. However, in 1990s, their contributions to ENSO simulation seem very minor except for the SST assimilation. 5.5 The impact of data assimilation on ENSO prediction Next we examine the predictions by the N H C M . A total of 40 prediction runs were made from April 1964 to January 1990 (the 'hindcast' period) and 35 prediction runs from Apri l 1990 to October 1998 (the 'forecast' period), starting at three months intervals (1 January, 1 April , 1 July, 1 October), and continued for 15 months. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 119 a without assimilation S S T assimilation co co CO O 1980 1985 1990 1995 2000 Year c S S H A assimilation 1980 1985 1990 1995 2000 Year co co CO O 1980 1985 1990 1995 2000 Year H C A assimilation 1980 1985 1990 1995 2000 Year Figure 5.7: Time evolution of observed and modeled Nino3 SSTA during January 1980-December 1997. Solid line denotes observed SSTA in all figures and dash hne is from the ocean model run (a) without assimilation, (b) with SST assimilation, (c) with SSHA assimilation, and (d) with H C A assimilation. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 120 Fig. 5.8a shows the N H C M predictive skills with and without data assimilation during 1980-1990, where the predicted Niho3 SSTA is compared against the observed values. Compared with the persistence skill (not shown), all predictions beat persistence at lead times greater than 2 months. In particular, SST assimilation beat persistence at all lead times. Generally, except for the predictions with SST assimilation, all suffer an initialization shock in the first 2 months, leading to lower skill than persistence. The prediction skills with the assimilations of SSHA and wind stress varies in a similar way with lead time, i.e. they rapidly decline with lead time, reaching a minimum at 9 months, beyond which their skills rebound and stabilize till a lead time of 14 months— in contrast to the predictive skills with the assimilations of SST, which basically simply decline with lead time. The best overall predictive skills were attained by H C A assimilation, with a correlation skill of above 0.7 at a lead time of 14 months. The wind-stress assimilation generated the second best performance with correlation of over 0.6 at a lead time of 14 months. SSHA assimilation, which did not contribute any additional skill at lead times under 6 months, was the third best at long lead times (10-15 months). Finally, SST assimilation which did best at 1-2 month lead times, finished last at long lead times, even worse than the predictions without data assimilation (which have correlations of above 0.5 at lead times of 10-14 months). For 1990-1998, Fig. 5.8b shows the skill improvements with data assimilation to be considerably less impressive than in the 1980s (Fig. 5.8a). Overall, H C A assimilation remained the best. SSHA assimilation, like in the 1980s, was consistently ahead of the predictions without assimilation at all lead times. Wind stress assimilation, unlike the 1980s, failed to beat predictions without assimilation at lead times of 12 months and beyond. SST assimilation again finished last at long lead times. Overall, the predictions with H C A assimilation are significantly better than the pre-dictions without assimilation. Wind-stress assimilation is generally the second best, Chapter 5. Impact of data assimdation on ENSO simulations and predictions 121 0.4 h without assi . wind-stress assi . - S S H A assi . * S S T assi . o H C A assi . 0.9 0.8 .2 0 7 O v_ o O 0.6 5 10 Lead time in months 15 0.5 0.4 0.3 * o without assi . wind-stress ass i . S S H A assi . S S T assi . H C A ass i . 5 10 Lead time in months 15 Figure 5.8: Correlation between observed and predicted SST anomalies in the Niho3 region, as a function of lead time for (a) the 1980s and (b) the 1990s. Solid line is from the prediction without data assimilation. The correlation from the predictions with data assimilations are shown by the dash hne for the wind stress assimilation, asterisks for SST assimilation, dot-dash for SSHA assimilation and circles for H C A assimilation Chapter 5. Impact of data assimilation on ENSO simulations and predictions 122 followed by SSHA assimilation, which shows some improvements over the predictions without assimilation, but not as much as found by Chen et al. (1998) and Fischer et al. (1997). The physical basis for ENSO prediction is that a delayed oscillator operates over the equatorial western Pacific and the equatorial eastern Pacific, as discussed in many studies (e.g. Schopf and Suarez 1988; Battisti 1988), which requires a phase shift between the upper ocean heat content variation in the equatorial western Pacific (Nifio4 region) and in the equatorial eastern Pacific (Niho3 region), with the variations in Niho4 leading those in Niho3. The ENSO forecast skills depend highly on how strong the phase shift relation is and how long the lag is. Fig. 5.9 shows the correlation between the monthly mean H C A in Nifio3 and Niho4 as a function of the time lag. As shown in Fig. 5.9b, for the observed H C A , the positive correlation coefficients are significantly higher in the 1980s than in 1990s for lags longer than 10 months, implying ENSO events could be predicted more accurately at long lead times in the 1980s than in the 1990s. In Fig. 5.9a, for the model with H C A assimilation, the positive correlations at positive phase lags are much higher in the 1980s than in the 1990s, indicating that the delayed oscillator mechanism was much stronger in the 1980s. Both Figs. 5.9a and b point to the 1990s as a period of lower forecast skills, and a period where assimilating H C A did not provide as great a boost in forecast skiUs as in the 1980s. The difference in prediction skills between the 1980s and the 1990s is not due to our having split the record into a 'hindcast' period (up to 1989) and a 'forecast' period (1990 onward) when we built our empirical model— as we tested another empirical model built with the observational data of 1980-1998, which still manifested the decadal difference in assimilation skills. What could be the cause for this difference of predictive skills between in 1980s and 1990s? It has been found in many studies that very persistent anomahes have occured in a number of Pacific oceanic and atmospheric variables in Chapter 5. Impact of data assimilation on ENSO simulations and predictions 123 Figure 5.9: (a) Correlation between the monthly mean heat content anomalies in the Niho3 and Niho4 regions as a function of time lag between the two time series, with H C A from (a) the model with H C A assimilation and (b) observations. Dashed line is for the 1980s and solid line for the 1990s. A positive time lag means that the signal in Niho4 is leading the signal in Niho3. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 124 1990s, which is associated with the decadal dependence of ENSO predictability (Kleeman et al. 1996; Latif et al. 1997; Ji . et al. 1996). Perhaps, the most striking of such anomalies has been SST and wind stress, characterized by persistent large-scale positive SST anomalies in the tropical Pacific and a weakening of the trade winds. A quite large positive anomaly also exists at 55 ° N off the Canadian west coast in contrast a cooling in the North Pacific. Corresponding to the changes in the ocean, many atmosphere variables also display persistent anomalies in 1990s, e.g the strongly enhanced convection in the equatorial dateline region, enhanced ITCZ and deepened Aleutian low and consistently negative Southern Oscillation Index (SOI) etc.. However it is still an open question how the persistent large-scale anomalies of ocean and atmosphere in 1990s affect the ENSO prediction skills. The predictions of the SSTA along the equator from various assimilation experiments are plotted in time-longitude diagrams as shown in Fig. 5.10 for the 1982/83 E l Nino. The predicted values at each grid point have been scaled to match the observed standard deviation at that grid. Compared with the observed SSTA, all experiments successfully predicted the anomalous warming in the equatorial central-eastern Pacific. However, except for the experiment with SST assimilation and that with wind stress assimilation, the other experiments all had an unrealistic, cool western equatorial Pacific, especially for the H C A assimilation experiment. In terms of the timing of the ENSO start, peak and decay, the experiment with H C A assimilation seems to be best. The experiments assimilating SST and wind-stress had the warming starting earlier than observed. The experiment assimilating SST and to a lesser extent, that assimilating SSHA did not decay much even by March, 1983. Compared with observations, the warming in the models are generally weaker at peak amplitude and too narrowly confined to between 110°W-150°W (except the one assimilating SST). Fig. 5.11 shows the time-longitude diagram of predicted SSTA along the equator for Chapter 5. Impact of data assimilation on ENSO simulations and predictions 125 without assimilation S S T assimilation S S H A assimilation 160E 160W 120W 80W Longitude wind-stress assimilation 160E 160W 120W 80W Longitude H C A assimilation 160E 160W 120W 80W Longitude observation 160E 160W 120W 80W Longitude 160E 160W 120W 80W Longitude 160E 160W 120W 80W Longitude Figure 5.10: Time-longitude diagrams of predicted SSTA along the equator for several assimilation experiments for the 1982/83 warm event. The predictions are all initialized on January 1, 1982. For comparison, the observed SSTA and the SSTA from the model without data assimilation are shown. The ordinate with values 1 to 15 is the calendar month from January 1982 to March 1983. Chapter 5. Impact of data assimilation on ENSO simulations and predictions 126 the 1988/89'La Nina. A l l assimilation experiments predicted the anomalous coohng in the equatorial eastern Pacific, but the coohng was too weak in the SST assimilation ex-periment, and too strong in the SSHA and wind stress assimilation experiments. The H C A assimilation experiment had generally the best performance, as it successfully pre-dicted two separate coolings, centered around July 1988 and December 1988, found also in the observations around June 1988 and November 1988. The models all had spurious warming in the western equatorial Pacific. For the 1997 E l Nino (Fig. 5.12), all experiments predicted weU the start time for this very strong event. The best performance was attained by H C A assimilation, which not only predicted well the peak time, decay time and the amplitude, but also avoided spurious coohng in the west, as found in most of the other experiments. SSHA assimila-tion also predicted well the mature time, decay time and the amphtude but had strong spurious cooling in the west. Again, aU predicted warmings were a little too narrowly confined around 120°W, when compared with observations. In summary, a 3D var assimilation system was developed to investigate the impact of assimilating different types of data on ENSO simulations and predictions. The data ex-plored were the sea surface temperature (SST), sea surface height anomahes (SSHA), wind stress, and the upper ocean 400 meter depth-averaged heat content anomahes (HCA) during 1980-1998. Neural network (NN) techniques have been used to relate data for the diagnostic variables (e.g. HCA) to prognostic variables in the model, so that the assimilation of non-prognostic data can be performed. Overall, H C A assimilation is the best in improving the prediction skills for lead times greater than 5 months. The improvement from SST assimilation is the best only for lead time of 4 months or shorter; at longer lead times, SST assimilation degrades the predictions. Wind-stress assimilation is generally the second best for lead times 6 months or longer, SSHA assimilation gen-erally produces prediction skills only shghtly better than the experiment without data Chapter 5. Impact of data assimilation on ENSO simulations and predictions 127 without assimilation S S T assimilation S S H A assimilation 160E 160W 120W 80W Longitude wind-stress assimilation 160E 160W 120W 80W Longitude H C A assimilation T T T 160E 160W 120W 80W Longitude 160E 160W 120W 80W Longitude \ : 1 I " \ \ 1 l \ • H I I , i , \ • i ' i 1 11 i \ 160E 160W 120W 80W Longitude observation \ I \ ' I-\: -1 I 160E 160W 120W 80W Longitude F i g u r e 5.11: Same as i n F i g . 5.10 but for the 1988/89 cool event. T h e p r e d i c t i o n was i n i t i a l i z e d on J a n u a r y 1, 1988 Chapter 5. Impact of data assimilation on ENSO simulations and predictions 128 without assimilation S S T assimilation S S H A assimilation 160E 160W 120W 80W Longitude wind-stress assimilation 160E 160W 120W 80W Longitude H C A assimilation 160E 160W 120W 80W Longitude observation 160E 160W 120W 80W Longitude 160E 160W 120W 80W Longitude 160E 160W 120W 80W Longitude Figure 5.12: Same as in Fig. 5.10 but for the 1997 warm event. The prediction was initialized on January 1, 1997 Chapter 5. Impact of data assimilation on ENSO simulations and predictions assimilation. Chapter 6 Constructing hybrid coupled models via 4-D variational data assimilation 6.1 Introduction In Chapter 2, two statistical atmospheric models were constructed using hnear regression and neural network, based on a traditional procedure: assuming a function form relating predictors (oceanic HC) to predictands (wind stress), the parameters in the function are estimated by optimally fitting data. The dynamics in the ocean are ignored in constructing the atmospheric model, which might lead to inconsistencies between the atmosphere model and the ocean model in dynamical couphng, leading to prediction errors. Is it possible to construct statistical relations while satisfying the dynamical constraints? 4D variational data assimilation provides a good technique for optimization under dynamical constraints. It has become popular in assimilating data into numerical models (Daley 1991; Ghil and Malanotte-Rizzoli 1991; Bennett 1992; Navon 1998). The method is commonly used to optimally estimate model parameters or initial conditions, and is being implemented in operational weather and climate prediction models (e.g., Zhu and Navon 1999; Gauthier et al. 1999; Courtier et al. 1998). The objective of this chapter is to estimate a statistical equation subject to the constraints of the dynamical equations of the system. As the dynamical equations in meteorology and oceanography are generally nonlinear, N N is used to model the statistical relation. 130 Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl31 The chosen dynamical model is the simple Lorenz (1963) model, with 3 variables, aris-ing from the Fourier truncation of the Rayleigh-Benard equations describing atmospheric convection. In the field of data assimilation, the celebrated Lorenz model has served as a test bed for examining the properties of various data assimilation methods ( Miller and Ghil 1990; Gauthier 1992; Miller et al. 1994; Evensen 1997) as the Lorenz model shares many common features with the atmospheric circulation and climate system in terms of variability and predictability (Palmer 1993). By adjusting the model parameters which control the nonlinearity of the system, the model can be used to simulate near-regular oscillations or highly nonlinear fluctuations. This chapter is structured as follows: Section 6.2 describes the hybrid Lorenz model, while section 6.3 shows some simple experiments involving the hybrid model. Section 6.4 gives a general formulation for coupling a N N to a dynamical model via variational data assimilation. Section 6.5 uses variational assimilation to estimate the N N parameters in the hybrid Lorenz model. Section 5.6 studies not only estimating the N N parameters, but also retrieving the dynamical model parameters and initial conditions. 6.2 Hybrid neural-dynamical Lorenz model The non-dimensionalized Lorenz (1963) nonlinear system of 3 differential equations is dX dY ^- = XY- cZ. (6.3) dt where variables X, Y, and Z are related to the intensity of convective motion, and the temperature gradients in the horizontal and vertical directions respectively. The Jx =-aX + aY, (6.1) dt = -XZ + bX-Y, (6.2) Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl32 parameters a, b and c will be referred to as dynamical parameters, in contrast to the empirical parameters of the N N model. The so-caUed observed data or true data are from the integration of the Lorenz equations (6.1)-(6.3) over 15,000 time steps at a step size h of 0.001, using a fourth-order Runge-Kutta integration scheme. As this system is very sensitive to changes in the initial conditions and model parameters, two cases are studied (Fig. 6.1): The first case, called the weakly nonlinear case, with the parameters a, b and c set to 10, 28 and 8/3 respectively, and initial conditions for (X, Y, and Z) to (-9.42, -9.43, 28.3) (as in Gauthier 1992), displays near-regular oscillations with a graduaUy increasing amphtude in the devised integration period. The other case is the highly nonlinear case, with a, b and c set to 16.0, 120.1 and 4.0 respectively (as in Eisner and Tsonis 1992), and initial conditions to (22.8, 35.7, 114.9). Next we assume the third Lorenz equation is unavailable, and we must approximate it empiricaUy with an N N equation. Our hybrid model thus consists of: ^ = -aX + aY, (6.4) at dY — = -XZ + bX-Y, (6.5) at ^ = NN(Xt,Yt,Zt). (6.6) where NN is a feed-forward N N model (Hsieh and Tang 1998). Xt, Yt and Zt are the 3 input neurons (also denoted by V{, i = 1,2,3), inputting the values of X, Y and Z at time t, and the single output neuron is dZ/dt. 6.3 Simple hybrid model experiments A simple-minded approach to the missing third Lorenz equation, is to use data to get the N N (6.6), and then integrate (6.4)-(6.6), with (6.6) as a replacement for the missing Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl33 0 5000 10000 15000 0 5000 10000 15000 Time steps Time steps F i g u r e 6.1: T h e L o r e n z sys tem in tegra ted over 15000 t ime steps for the w e a k l y non l inear case (left co lumn) and the h igh ly nonl inear case (right c o l u m n ) . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation 134 Lorenz equation. The N N has X, Y, and Z as the inputs and (Zt+1 — Zt)/h as the target during training (with the optimization minimizing the M S E between the target and the model output), so the model output is approximately dZ/dt. The optimization of the N N model was determined over a training period of 3000 time steps, and the following period of 1000 time steps was used to independently test the N N model skills in estimating dZ/dt. The appropriate number of neurons in the hidden layer is chosen based on the trade-off between under-fitting (too few neurons) and over-fitting (too many neurons) as in Chapter 2. After trying models with different numbers of hidden neurons, we found that 5 hidden neurons yielded the best skills over the test period. To alleviate the instability problems associated with N N modelling (Hsieh and Tang, 1998), an ensemble of 25 identical NNs were trained with random initial weights and biases, with the outputs from the 25 NNs used to provide an ensemble mean output. The skills of the ensemble mean generally exceeded the skills of an individual member (Hsieh and Tang 1998). The ensemble mean N N simulations of dZ/dt for both the weakly nonhnear and the highly nonlinear cases during the training period and the test period are shown in Fig. 6.2, where the N N gave good simulations of dZ/dt, although the oscillation peaks were underestimated. Now that the third Lorenz equation seems to be reasonably approximated by the N N (6.6), we proceed to integrate the hybrid system (6.4)-(6.6) from perfect initial conditions and compare the output with respect to that from the true Lorenz system (6.1)-(6.3). This simple-minded hybrid approach worked reasonably well in simulating the true Lorenz system for the weakly nonlinear case (Fig. 6.3), but failed completely for the highly nonlinear case, where the hybrid model yielded X and Y values which were off by orders of magnitude. The highly nonhnear Lorenz system is known for its extremely high sensitivity to perturbations. Clearly, the N N approximation to the third Lorenz equation introduced significant errors in Z, which quickly caused the X and Y Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation 135 F i g u r e 6.2: T h e ensemble m e a n of 25 N N s (dot-dash curve) for a p p r o x i m a t i n g the t h i r d L o r e n z equa t ion (sol id curve) for the weak ly nonl inear case (left co lumn) a n d the h i g h l y non l inear case (right c o l u m n ) . Top panels (a) and (b) are for the t r a i n i n g p e r i o d of 3000 t i m e steps, and b o t t o m panels (c) and (d) for test p e r i o d of 1000 t i m e steps. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationliQ 0 1000 2000 3000 0 100 200 300 400 500 Time steps Time steps Figure 6.3: The hybrid model (dot-dash curve) and the true Lorenz model (solid curve) integrated from identical initial conditions for the weakly nonhnear case (left column) and the highly nonlinear case (right column). values to deviate far off course. The defect of this simple approach lies in the fact that the optimization of the N N was achieved without imposing the dynamical constraints of (6.4) and (6.5). A better approach is variational assimilation, where the dynamical constraints are imposed. 6.4 The hybrid model using variational data assimilation The scalar equations (6.4)-(6.6) can be expressed as a vector equation: Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl37 ^ = f ( v , P > 0 (6.7) where v is the vector denot ing (X,Y,Z), and p is the paramete r vec tor ( w h i c h c o u l d con t a in the parameters o f the N N or o f the d y n a m i c a l equat ions) . T h e va r i a t i ona l ass imi la t ion m e t h o d involves m i n i m i z i n g a quadra t i c cost f unc t i on J subject to m o d e l const ra int (6.7), where J is defined as J ( v , p , v 0 ) = [T D{v,t)dt (6.8) Jo D= ( v - V o b . y ' W - ^ v - V o b . ) (6.9) where the subscr ip t obs denotes the observed values, the superscr ip t tr the t ranspose, T the l e n g t h of the ass imi la t ion w indow (also ca l led t r a i n i n g pe r iod ) , Vo the i n i t i a l cond i t ions , a n d W the covariance m a t r i x of the measurement errors, assumed here to be d iagona l . I n t roduce the Lagrange func t ion L, L = J+ [Tv*tr[^-f}dt (6.10) Jo at where v*(i) is a vector of Lagrange mul t ip l i e r s (or adjoint var iables) . A f t e r i n t eg ra t ion by par ts , L becomes TT r1v*tr L = [ v * % £ + Jo[D- ( — v + v*trf)]dt (6.11) T h e first order va r i a t ion of L is SL = [v*tr8v}l + / [ V v M v - [ ( -77 - ) < r Jv + (5-)* pv'*v + (?-)trv*8p}dt (6.12) Jo at ov ap T h e h y b r i d adjoint m o d e l can be ob ta ined by le t t ing 8L/8v = 0: Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl38 v*(T) = 0. (6.14) According to (6.13) and (6.14), the formulas for computing the gradients of J with respect to p and v 0 can be obtained by differentiating (6.12) with respect to these un-knowns (Lu and Hsieh 1998): < 6 1 5 > r j ^ = - V ( 0 ) . (6.16) Details of the variational assimilation are given in Appendix B . As the above assimilation process, also referred to as the training process, is com-pletely subject to the model (6.4)-(6.6) after the initial guesses are given, it is caUed the strong continuity constraint assimilation scheme (SC). 6.5 Determining N N parameters via variational assimilation In this section, we determine the N N parameters, i.e., the weights and biases, in the hybrid model by variational data assimilation. A N N with 5 hidden neurons is used. Preliminary experiments are done in Section 6.5.1 and 6.5.2 to determine the best first guesses and assimilation window sizes, respectively, before the main experiments are run in Section 6.5.3. In this section, all experiments were performed with perfect initial conditions and dynamical parameters a and b. 6.5.1 The initial guesses Whether an assimilation process is successful depends greatly on the choice of the first guesses. To help get a good initial guess, two assimilation schemes are introduced in Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl39 addition to the SC scheme. Under SC, at each step of the integration, the initial con-ditions are the model outputs from the previous step (Fig. 6.4a). In the no continuity constraint assimilation scheme (NC), at each integration step, the initial conditions are the observed data, rather than the model output from the previous step (Fig. 6.4b). Between these two extremes, we can have the partially strong continuity constraint assimilation scheme (PSC) (Fig. 6.4c), where the assimilation window T is divided into smaUer assimilation segments (AS). Within each AS, integration at each step uses the model output from the previous step as the initial conditions. Only at the start of the AS, are the data used directly as the initial conditions. Both SC and NC can be regarded as special cases of PSC, where the AS is T in SC and 1 time step in NC. The reason for introducing NC and PSC is that using SC without good initial guesses generally does not lead to successful assimilation over a long assimilation window. Chop-ping the window into shorter AS means that the dynamic model constraints are no longer imposed over a long period, thereby making the nonhnear optimization a much easier task. Our strategy is to divide the assimilation process into 2 stages: (a) use less de-manding schemes such as NC and PSC to obtain reasonable parameters estimate, which will then be used as initial guesses in (b) the full variational assimilation under SC. For the weakly nonlinear case, the NC assimilation scheme was first applied to the hybrid model (6.4)-(6.6), with initial guesses for the N N parameters taken randomly from one of the 25 NNs in Section 6.3. The assimilation window is 3000 time steps. The N C scheme can produce good results during the training period but poor skill during the test period (Fig. 6.5). Since the AS is very short (1 time step) in the N C scheme, the initial values at each integration step vastly outweigh the importance of the dynamical constraints, hence the poor skill during the test period. However the result from the N C can provide reasonable first-guesses for the following PSC. A series of PSC assimilations with AS = 100, 200, 500 and 1000 time steps was Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlAO t. * • t t Figure 6.4: Schematic diagram for the (a) SC (Strong Continuity constraint), (b) N C (No Continuity constraint) and (c) PSC (Partial Strong Continuity constraint) assimilation schemes. The solid curves are the model trajectory, dots, the observed values and dash lines, the model errors to be minimized. Here a model trajectory starts from an observed value. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlil 24 L 2 0 1 , , , , — _ 1000 2000 3000 0 200 400 600 800 1000 Time step Time step F i g u r e 6.5: T h e N C ass imi la t ion results , w i t h the m o d e l t r a i n i n g results shown o n the left panels , a n d pred ic t ions over the test pe r iod on the r ight panels; where the sohd curves show the t rue da t a and dot-dash curves the outputs of the h y b r i d m o d e l . T h e do t -dash curves over lap the sol id curves i n the left panels. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlA2 performed with the hybrid model (in the weakly nonhnear case) to gradually improve the first guesses to be used in the stage (b) SC assimilation. A longer AS increases the importance of dynamical constraints relative to initial conditions in influencing the model outputs. The first guesses of the N N parameters in a PSC experiment are the parameters obtained from the previous PSC experiment with the shorter AS. The first PSC experiment (not shown) with AS = 100 time steps uses the results from the N C case (i.e. AS = 1 time step) for its first guess. The PSC assimilation over the training period, i.e. a window of 3000 time steps, with AS = 200 time steps (Fig. 6.6), and AS = 1000 time steps (Fig. 6.7), reveal that the test skills improved with the increase of the AS. The parameters estimated from the PSC experiment with AS= 1000 time steps will serve as the first guess for the stage (b) SC assimilation (in the next Section) for the weakly nonhnear case, where the assimilation window is also taken to be 1000 time steps. This approach of using a series of PSC assimilations to improve on the parameter estimates to be used as initial guesses in the stage (b) SC assimilation did not work for the highly nonhnear case. Unlike the weakly nonhnear case, the first guesses from PSC, if found, is only effective for a specific experimental configuration (e.g. AS, window length, etc.). Hence, for the highly nonhnear case, the parameters estimated from the N C assimilation over 3000 time steps are directly used as the first guesses for the stage (b) SC assimilation. 6.5.2 Impact of assimilation window in the highly nonlinear case The greatest difficulty in variational assimilation, namely the presence of multiple min-ima in the cost function arising from the nonhnearity of the problem, often renders the problem intractable. The number of local minima dramatically increases with the window T for a strongly nonlinear dynamical system. Gauthier (1992) and Miller et al (1994) Figure 6.6: The PSC assimilation results with AS= 200 time steps. The model training results are shown on the left, and model predictions on the right, with solid curves for the true data and dot-dash curves for the outputs of the hybrid model. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationliA Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlA5 have shown with the Lorenz model that variational assimilation is effective only for suf-ficiently short windows. In the weakly nonlinear case (Section 6.5.1), T did not have a major impact; hence, the following discussion will only focus on the highly nonhnear case. To study the impact of the window, SC assimilation experiments were performed, where T was varied between 100, 300, 500 and 800 time steps. For each T, 10 experiments were performed with the same first guess, but the assimilation periods were shifted by 100 time steps between one experiment and the next. The effect of T on the assimilation is different for training period (Fig. 6.8a) and the test period (Fig. 6.8b). With a short window (T = 100 time steps), the correlation skills were excellent for all 3 variables during training, but yielded no prediction skiUs over the following test period of 300 time steps. As T increased to 300 time steps and then 500 time steps, there were increases in the prediction skills (Fig. 6.8b), but as T increased to 800, the prediction skill dropped sharply, as the problem with local minima in the cost function worsened with increasing T. Hence in the following SC assimilation experiments, the window will be 500 time steps for the highly nonhnear case, and 1000 time steps for the weakly nonhnear case. FoUowing the training (assimilation) period, there is a test period of 300 time steps for the highly nonhnear case, and 1000 time steps for the weakly nonlinear case. 6.5.3 SC assimilation We now perform the main SC assimilation experiments, repeated 100 times, to examine the skills of the hybrid model for the weakly nonhnear and highly nonhnear cases. The choices of windows and first guesses are those found in Sections 6.5.1 and 6.5.2. The 100 experiments are implemented with their assimilation periods shifted by 100 time steps from one experiment to the other. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation!^ Var iab les Figure 6.8: Correlation skills (averaged over 10 experiments) for the variables X, Y and Z at various assimilation windows (T = 100, 300, 500 and 800 time steps) during (a) the training period, and (b) the test period. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl4.7 For the weakly nonlinear case, the average correlation between model and data (over 100 experiments) exceeds 0.96 for each of X, Y and Z during both the training period and the test period (Fig. 6.9a). The relative estimate error, R E E , i.e. (estimated value - true value) 2/ (true value)2) is very small, less than 0.004 (Fig. 6.9b). Thus, via variational assimilation, a hybrid model has successfuUy reconstructed the Lorenz model in the weakly nonlinear case. For comparison, the results of the simple hybrid model without data assimilation from Section 6.3 are also show in Fig. 6.9. The improvement due to variational data assimilation is dramatic for X and Y. In contrast, the highly nonlinear case has much lower skill. During training, the skills were good- 80% of the experiments had correlation above 0.75, and 60% had R E E under 0.05 (not shown). However, during the test period, the correlation skiU and R E E attained were much poorer (Fig. 6.10). Generally, for the highly nonlinear case, there were three classes of assimilation results found in the 100 experiments. A typical example of the class 1 results is shown in Fig. 6.11, showing successful assimilation. About 15% of the total experiments are of this class. Almost 60% of the total experiments belong to a class 2 of moderately successful assimilation, as illustrated by the example in Fig. 6.12, where during the test period, the forecast skills were fairly good for the first 200 time steps. The remaining 25-30% of the experiments belonged to a class 3 of failed forecasting, as illustrated in Fig. 6.13, where despite excellent fitting during the training period, the forecast results were totally wrong. To see why there is a drastic difference between class 1 and class 3 behaviour, we plotted the system trajectories of the two cases in the X — Y — Z space (Fig. 6.14). The Lorenz attractor for the highly nonhnear case is known to have a 'butterfly' shape with two wings. In the class 1 experiment, both wings of the Lorenz attractor were covered by the training data trajectory. In contrast, with the class 3 example, the training data Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation 148 Q Training (with assimilation) 23 testing (with assimilation) • training (without assimilation) E3 testing (without assimilation) .9> 0.8 0.6 0.4 0.2 0 X Y V a r i a b l e s Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlA9 LU LU CC 0.07 0.06 0.05 0.04 0.03 0.02 0.01 H Training testing • Training • testing (with assimilation) (with assimilation) (without assimilation) (without assimilation) Variables Figure 6.9: (a) Correlation skills and (b) R E E for the variables X, Y and Z during the training and test periods (averaged over 100 experiments), for the weakly nonlinear case, where the SC assimilation retrieves the N N parameters. The error bars show the standard deviation. Results from the simple hybrid model without data assimilation (Section 3) are also shown for comparison. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlbO a £ o .a £ 0.9 0.8 0.7 0.6 0.5 0.4 Correlation trajectory resided in one wing. The system then evolved to the other wing during the testing period. As the model had never been trained with data from the second wing, it obviously could not forecast properly in that regime, hence the disastrously poor forecast results found in Fig. 6.13. Thus the bimodality of the highly nonlinear Lorenz system created cases where the training data covered only one wing of the attractor, and the model thus trained had no capability of forecasting the transition to the other wing. One would be tempted to think that the problem could be corrected by simply using longer training periods, so that the full Lorenz attractor would have been learned by the model. Unfortunately, a longer training period greatly increases the difficulty in the nonlinear optimization due to multiple local minima in the cost function, leading eventually to the failure of the variational assimilation method. Figure 6.10: Bar charts from 100 assimilation experiments for the highly nonhnear case showing the number of experiments with (a) correlation above a particular correlation level and (b) R E E below a particular R E E level, during the test period. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl52 50 X 0 -100 i 200 0 100 200 300 400 500 e -20L 0 100 200 300 400 500 f 160 r F i g u r e 6.11: A t y p i c a l exper iment w i t h class 1 results for the h i g h l y non l inear case, showing (a) the fitting d u r i n g the t r a in ing pe r iod (left panels) a n d (b) the forecast ing d u r i n g the tes t ing p e r i o d (right panels) . Sohd curves denote the t rue value, wh i l e dot-dashed curves denote m o d e l results . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation!53 F i g u r e 6.12: Same as F i g . 6.11 but for a t y p i c a l exper iment w i t h class 2 (i.e. m o d e r a t e l y successful forecast) results . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl54 F i g u r e 6.13: Same as F i g . 6.11, but for a t y p i c a l exper iment w i t h class 3 (i.e. unsuccessful forecast) resul ts . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation!^ F i g u r e 6.14: T h e sys tem trajectories i n the X — Y — Z space d u r i n g the t r a i n i n g p e r i o d ( t h i n curves) and the tes t ing p e r i o d ( th ick curves) for (a) the class 1 expe r imen t a n d (b) the class 3 expe r imen t . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlSG 6.6 Determining dynamical parameters and initial conditions by the hybrid model In variational data assimilation, parameters and initial conditions can be jointly estimated (e.g., Le Dimet and Talagrand 1986; Tziperman et al 1989; Lu and Hsieh 1998; Zhu and Navon, 1999). I next conducted an experiment, for the weakly nonlinear case, to estimate simultaneously the dynamical parameters a and 6, and the initial conditions of X,Y and Z in the hybrid Lorenz system (6.4)-(6.6). The N N from Section 6.5.1 (i.e. PSC assimilation with AS = 1000 time steps) was used. For the initial guess, a and b and the initial conditions were all scaled down by 10% from their true values. The result of retrieving the 2 parameters and 3 initial conditions by the hybrid model with an assimilation period of 3000 steps under SC showed the retrieval to be very good, as the REEs were all of the order of 1 0 - 3 - 10~4 (not shown). Next the true values were scaled down by 20% to provide the initial guesses, and the retrieved results were still generally good (Fig. 6.15). But when the initial values were the true values scaled down by 30%, the retrieval nearly failed (not shown). Next, we had the N N parameters retrieved as well during the assimilation, i.e. we conducted 100 additional SC experiments to estimate simultaneously the N N parameters, the dynamical parameters a and b, and the initial conditions of X, Y and Z, for the weakly nonlinear case. The first guess for the initial conditions and the parameters a and b were simply the true values scaled down by 10%. The choice of the first guess of the N N model parameters, cost function and assimilation window were the same as those discussed in the previous section. The average correlation skills over the 100 experiments are shown in Fig. 6.16, while the average R E E are all around 0.05 (not shown). The correlation and R E E were generally poorer than those in Fig. 6.8, where only the N N parameters were retrieved. This agrees Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation\57 F i g u r e 6.15: T h e results f rom re t r i ev ing the d y n a m i c a l parameters a, b and the i n i t i a l condi t ions of the h y b r i d m o d e l , where the i n i t i a l guesses were the t rue values scaled d o w n by 20%. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl58 1.4 1.2 0.8 as 0.6 8 0.4 0.2 0 n Training Testing fW£W,J WW.W A - J W-JI-JI-J A-J"« J--J - . J i . i . i . / • S i ® i S ' S > S > S > ' k£wtf www W£W : v v v y • • ^ • ^ • ^ • ^ W W W www i • • V"V_"V «^  • T^i-1 ?1 iw: mm'mm* w< w< w< W i wikiv.i www ww# W W W www www WWtf-www W W W i . - • . . ••• X V a r i a b l e s F i g u r e 6.16: T h e cor re la t ion s k i l l averaged over 100 exper iments where the N N parame-ters, the d y n a m i c a l parameters a and 6, and the i n i t i a l condi t ions were r e t r i eved for the w e a k l y nonhnear case. E r r o r bars ind ica te the s t anda rd dev ia t ion . Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationl59 0.9 0.8 0.7 Cor re la t ion 0.6 0.5 0.4 with many studies which found that the assimilation quality deteriorated while retrieving more model unknowns (Thacker and Long 1988; Lu and Hsieh 1998), as the problem of multiple minima in the cost function might become worse when retrieving more model unknowns. Likewise, 100 experiments for the highly nonlinear case were also carried out to si-multaneously estimate the N N parameters, the dynamical parameters a and b and the initial conditions. The first guesses were the true values scaled down by 10%, and the first guess for the N N model, the cost function and assimilation window were again taken from Section 6.5. While the skills in the training period were as good as those retrieving the N N parameters only (not shown), the skills during the testing period (Fig. 6.17) were considerably lower than those found when only N N parameters were retrieved (Fig. 6.10). Class 1 examples like Fig. 6.11can hardly be found now. In summary, a hybrid neural-dynamical variational data assimilation procedure aimed Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilation!60 1 00 w 80 A—-CD CD CL X LU CD E 60 40 20 0.1 0.2 0.3 0.4 0.5 REE Figure 6.17: Bar charts from 100 assimilation experiments for the highly nonlinear case showing the number of experiments with (a) correlation above a particular correlation level and (b) R E E below a particular R E E level, during the test period. The N N param-eters, the dynamical parameters a and b, and the initial conditions were retrieved. Chapter 6. Constructing hybrid coupled models via 4-D variational data assimilationlQl at simulating missing dynamical equations has been formulated in this chapter. This pro-cedure was then applied to the Lorenz model, where the third dynamical equation of the model was missing and had to be simulated empirically by an N N equation. Several important issues in 4 D var assimilation were discussed with a series of numerical ex-periments. Overall the numerical experiments showed that the hybrid model based on N N and variational assimilation is successful in simulating the weakly nonlinear Lorenz model, with forecast skills similar to the original Lorenz model. For a highly nonlinear Lorenz system, the hybrid model can also produce reasonable prediction skills for most experiments, although the assimilation basically failed for a fair portion of the experi-ments. Chapter 7 Summary and Discussion In this thesis, two hybrid coupled models were developed by coupling two different at-mospheric models to a 6-layer nonlinear dynamical ocean model. Neural network and linear regression were used to construct the two atmospheric models— one nonlinear and the other linear. A variety of numerical experiments were carried out to investigate the two coupled models in terms of ENSO oscillations and ENSO predictions. To explore the effects of assimilating different types of data for initializing the model for ENSO pre-dictions and to further improve model prediction skills, a three-dimensional variational system was developed to assimilate various observational data into the coupled model. The possibility of a hybrid coupled model with '4D' variational data assimilation was tested with the simple Lorenz model. 1. Empirical atmospheric models and the oceanic responses to them In Chapter 2, the possibility of using a nonlinear empirical atmospheric model for hybrid coupled modelling was investigated by developing a neural network (NN) model for predicting the contemporaneous wind stress field from the ocean state, and comparing the N N model with a linear regression (LR) model. Upper ocean heat content (HC) from an ocean model was found to be a better predictor of the wind stress than the (observed or modelled) SST. The N N model generally had slightly better skills in predicting the contemporaneous wind stress than the simple L R model in the off-equatorial tropical Pacific and in the eastern equatorial Pacific, mainly through better predictions of the second and third wind stress modes. 162 Chapter 7. Summary and Discussion 163 When the N N and L R model produced wind stresses were used to drive the ocean model, slightly better SST skills were found in the off-equatorial tropical Pacific and in the eastern equatorial Pacific when the N N winds were used instead of the L R winds. Better skills for the model HC were found in the western and central equatorial Pacific when the N N winds were used instead of the L R winds. Because changes in the HC in the western equatorial region can lead to SST anomalies in eastern Pacific, the potential skill improvement for HC by N N could be of interest. There are several possible reasons why N N failed to show more significant improve-ments over L R in the equatorial Pacific for wind stress and SST. The first (and the most probable) is that the relation between the surface ocean and the atmosphere in the equatorial Pacific over the seasonal time scale is basically linear, with nonlin-ear processes playing only minor roles. The linear assumption was also supported by Xue et al. (1994) and Tang et al. (2000). Another reason why N N failed to show more significant improvements over L R in the equatorial Pacific for wind stress and SST is that the data records are not long enough. N N is more capable than L R but the data requirement is also considerably higher. To extract more than the linear rules from the data, longer records of good quality data are needed. As Barnett et al. (1993) built L R models for individual calendar months, we also tested empirical atmospheric models for individual sea-sons. Although the seasonal approach can include the effect of the different basic states in different seasons of the year, we did not find the overall cross-validated anomaly skills better than the model developed using all the data - the tradeoff being that by dividing the data record into seasons, we have even less data to construct each seasonal model. 2. E N S O oscillations in the two hybrid coupled models Chapter 7. Summary and Discussion 164 The characteristics of the N H C M (with an N N atmosphere) and the L H C M (with an L R atmosphere) in terms of ENSO oscillations were investigated in Chapter 3. The N H C M can simulate a realistic chmatology and the seasonal cycle, and captures the major features of ENSO variability. The L H C M has a longer ENSO oscillation period than the more realistic period found in the N H C M , and a very narrow phase-locking to the annual cycle compared to the more realistic phase-locking in the N H C M . Model sensitivity studies demonstrated that both the N H C M and the L H C M are sensitive to the coupling strength of the system, with the L H C M being more sen-sitive. Depending on the changes in the coupling parameter, both coupled models can change from damped oscillations to periodic oscillations. With the increase of the coupling strength, the seasonal cycles have more effect on the inherent os-cillation of the coupled models. As the parameter reaches the critical value, both models are locked in a biennial frequency oscillations peaking in winter. The chaos mechanism has been proposed in many ENSO models. Our sensitiv-ity studies do not exhibit chaotic behaviour for either coupled model. Instead, the coupled models will shift to another climate regime as the coupling parame-ters increase beyond their critical values. There are two possibilities why the two coupled models did not exhibit chaotic behaviour: (1) There might exist some crit-ical parameters which we failed to find, which would produce chaotic behaviour. (2) Chaotic behaviour might not play a central role in ENSO oscillations, instead stochastic forcing may be the source of ENSO irregularity. Several recent works do show that the random forcing may play a central role in ENSO evolution (Blanke et al. 1997; Eckert and Latif 1997). The irregular interannual oscillation simulated by coupled CGCMs (Philander et al. 1992) does not appear to be low-order chaos. Chapter 7. Summary and Discussion 165 The nonlinear time series analysis performed by Chang et al. (1996) suggests that stochastic processes, rather than chaotic dynamics, are likely to be a major source of ENSO irregularity in CGCMs and in nature. 3. ENSO predictions in the two hybrid coupled models In Chapter 4, the prediction skills for the two coupled models, N H C M and L H C M , with two different initialization schemes were studied. The ensemble of 104 'hind-casts' experiments from 1964-1990 and 35 'forecasts' experiments from 1990-1998 shows that both models yielded useful predictive skills for the prediction of the equatorial east Pacific SST anomalies up to 1 year lead time. Two initialization schemes were tested: The stress-from-ocean initialization scheme, which considers the ocean-feedback in the initial conditions, generally has better predictive skills than the stress-only scheme, where the ocean is simply forced by the observed wind stress. The stress-from-ocean scheme also manifests less decadal variability in the forecast skills than the stress-only scheme. From 1964-1998, with the stress-from-ocean initialization scheme, the forecast correlation skill at 12-month lead time for the NIN03 sea surface temperature anomalies (SSTA) is about 0.55 for the N H C M and 0.50-0.55 for the L H C M . The main difference in the forecast skills between the N H C M and the L H C M occurs in the 1990s, where the N H C M has better skills. A nonlinear canonical correlation analysis of the wind stress and the SSTA shows that the relation between the two was indeed more nonlinear in the 1990s than in the 1980s, which would give the N H C M an advantage over the L H C M in the 1990s. The predictions for 3 typical ENSO episodes in 1982/83, 1997 and 1998 shows that both models can fairly predict these events in terms of their mature time and strength of anomalies around 1 year in advance. But both models displayed weakness in predicting the start time and end time, especially for the L H C M . Chapter 7. Summary and Discussion 166 4. The impact of data assimilation on ENSO simulations and predictions The impact of assimilating different types of data on ENSO simulations and pre-dictions has been studied using a 3D Var assimilation scheme. The data explored were the sea surface temperature (SST), sea surface height anomalies (SSHA), wind stress, and the upper ocean 400 meter depth-averaged heat content anomalies (HCA) during 1980-1998. Neural network (NN) techniques have been used to relate data for diagnostic variables (e.g. HCA) to prognostic variables in the model, so that the assimilation of non-prognostic data can be performed. In terms of ENSO simulation (SSTA simulation), the most pronounced improve-ment in correlation skills resulted from SST assimilation as expected, whereas the skill improvement from H C A assimilation was less significant than SST assimila-tion but was still considerable. SSHA assimilation did not significantly improve the SSTA simulation. ENSO predictability is strongly influenced by the phase shift relation between the variations of the upper ocean heat content in the eastern equatorial Pacific (Nifio3 region) and in the western equatorial Pacific (Niho4 region). Compared with the phase shift relation in the 1980s, both the observed H C A and the assimilated model H C A contain less pronounced phase relations between the Niho3 and Nifio4 regions in the 1990s, leading to less significant improvements to ENSO prediction with data assimilation in the 1990s than in 1980s. Overall, H C A assimilation is the best in improving the prediction skills for lead times greater than 5 months. The improvement from SST assimilation is the best only for lead time of 4 months or shorter; at longer lead times, SST assimilation degrades the predictions. Wind-stress assimilation is generally the second best for lead times 6 months or longer, though it degrades the prediction skills for lead times Chapter 7. Summary and Discussion 167 longer than 12 months in 1990s. SSHA assimilation generally produces prediction skills only slightly better than the experiment without data assimilation. 5. Hybrid coupled models via 4-D variational data assimilation Proceeding from 3-D var to 4-D var, a hybrid neural-dynamical variational data assimilation procedure aimed at simulating missing dynamical equations was for-mulated in Chapter 6. This procedure was then applied to the Lorenz model, where the third dynamical equation of the model was missing and had to be simulated empirically by an N N equation. First guesses in variational data assimilation can be vital to retrieving the model pa-rameters and initial conditions successfully. In order to get reasonable first guesses, the N C (No Continuity constraint) and PSC (Partial Strong Continuity constraint) schemes were proposed in this study. Experiments show that such a treatment of choosing first guesses is effective in the weakly nonlinear case, which allows the cost function to be successfully optimized in almost all experiments under the SC (Strong Continuity constraint). For the highly nonlinear case, the N C scheme pro-vides reasonable first guesses for the SC experiments. The length of the assimilation window is very significant for the highly nonlinear case. When the window is either too short or too long, the prediction skills are poor. This strict demand for an appropriate assimilating window results from a delicate balance between the need for more data (longer window) and the avoidance of severe multiple minima during optimization with the highly nonlinear model, which calls for a short window. Our numerical experiments showed that the hybrid model based on N N and vari-ational assimilation is successful in simulating the weakly nonlinear Lorenz model, with forecast skills similar to the original Lorenz model. For a highly nonlinear Chapter 7. Summary and Discussion 168 Lorenz system, the hybrid model can also produce reasonable skills of simulations and predictions for most experiments, although the assimilation basically failed for a fair portion of the experiments. Deterministic optimization algorithms often have difficulties reaching the global minimum in the cost function, due to the presence of local minima. It has been shown that the variational assimilation with deter-ministic optimization algorithms often loses its power for strongly nonhnear models (Miller et al. 1994; Evensen 1997). Stochastic optimization methods such as sim-ulated annealing may help (Kruger 1993; Hannachi and Legras 1995). However, even if the global minimum is found, the data in the assimilation window might still only consist of data from one wing of the butterfly-shaped attractor, whereas the data in the forecasting test period could he in the second wing, thereby resulting in failed forecasts. Alternatively, other assimilation approaches, e.g. the Extended Kalman filter might be more powerful in deahng with strongly nonhnear systems (Miller et al. 1994; Evensen 1997; Houtekamer and Mitchell, 1998). This thesis has initiated the fusion of neural network techniques into dynamical mod-els. Using an N N model for the atmosphere, it has produced the first H C M with a nonhnear empirical atmospheric component, and showed that the nonhnear atmosphere could have advantages over a linear atmosphere in modelling ENSO variability and in ENSO prediction. This study has introduced N N for the assimilation of non-prognostic variables (e.g. upper ocean heat content) by using N N to relate the non-prognostic vari-able to prognostic variables, which are then assimilated into the model. An operational ENSO prediction system has been developed and is starting to be used to issue routine ENSO prediction (Tang and Hsieh, 2001d). The forecast correlation skills of the N H C M model without data assimilation are comparable to those attained by the pure N N model for lead time of 12 months or less, but surpass the pure N N model skiUs for lead time Chapter 7. Summary and Discussion 169 exceeding 12 months— e.g. the N H C M with the stress-from-ocean initiahzation scheme has a correlation skill of 0.49 at a lead time of 15 months versus a correlation of 0.38 for the pure N N (Tang et al. 2000). With the best data (HC) assimilated, the prediction system generates noticeably better skills than the pure N N model at long lead times, and compares very favorably with the most influential prediction models in the ENSO community (Chen et al. 1995, Chen et al 1998, J i . et al. 1998, Syu and Neehn 2000). Table 7.1 is a comparison of correlation skills between the observed SST NIN03 index and that from the Lamont model, described in Chen et al (1995, 1998) and the N H C M with HC assimilation for several different periods. The Lamont model has been used routinely in real time for ENSO forecast since 1986. In particular, it was one of the few real-time forecast models that predicted the onset and evolution of the 1991-1992 warm ENSO phase correctly. The Lamont model has been widely used in prediction and predictability studies (Latif et al. 1998). Table 7.1: Correlation between the observed SST NIN03 index and that from Lamont model, described in Chen et al (1995, 1998) arid the N H C M with HC assimilation Time period Model Lead times 3 months 6 months 9 months 12 months 15 months 1992-1997 Lamont 0.91 0.78 0.40 0.39 0.60 N H C M 0.80 0.74 0.66 0.65 0.62 1988-1995 Lamont 0.69 0.60 0.50 0.49 0.49 N H C M 0.71 0.69 0.57 0.54 0.53 1980-1987 Lamont 0.90 0.86 0.85 0.80 0.80 N H C M 0.78 0.77 0.75 0.75 0.73 As shown in Table 7.1, the N H C M with HC assimilation has better performance in 1990s than the Lamont model, in particular for lead time greater than 9 months. But in the 1980s, the Lamont model has higher correlation skills than the N H C M . However, it should be noted that such a comparison is not completely fair since different samples Chapter 7. Summary and Discussion 170 were used in calculating the correlation skills for the Lamont model and the N H C M . While the fuU 4-D var H C M is beyond the scope of this thesis, a neural-dynamical hybrid approach under 4-D var has been developed to study the simple Lorenz system. Chapter 8 References Anderson D. L. T., and McCreary, J. P., 1985: Slowly propagating disturbances in a coupled ocean-atmosphere model. J. Atmos. Sci. 4%i 615-629 —, and J . Willebrand (eds), 1989: Oceanic Circulation Models: Combining Data and Dynamics, Kluwer Academic Publishers, 287-302. Balmaseda. M . A. , D. L. T. Anderson, and M . K . Davey, 1994: ENSO prediction using a dynamical ocean model coupled to statistical atmospheres, Tellus, 46(A) 4, 497-511 —, M . K . Davey, and D. L. T. Anderson, 1995: Decadal and seasonal Dependence of ENSO prediction skill. J. Climate, 8, 2705-2715 Barnett, T., and R. Preisendorfer, 1987: Origins and levels of monthly forecast skill for United States surface air-temperature determined by canonical correlation analysis, Mon. Wea. Rev., 115, 1825-1850. —, M . Latif, N . E . Graham, M . Fliigel, S. Pazan and W. White, 1993: ENSO and ENSO related predictability, Par t i : Prediction of equatorial sea surface temperature with a hybrid coupled ocean-atmosphere model. J. Climate, 6, 1545-1566 Battisti, D. S., 1988: Dynamics and thermodynamics of a warming event in a coupled tropical atmosphere-ocean model. J. Atmos. Sci., 45, 2889-2919. 171 Chapter 8. References 172 —, and A. C. Hirst, 1989: Interannual variability in the tropical atmosphere/ocean system: Influence of the basic status, ocean geometry and nonhnearity. J. Atmos. Sci., 46, 1687-1712. —, and E. S. Sarachik, 1996: Understanding and predicting ENSO. U.S. Report to International Union of Geodesy and Geophysics: Contributions in Oceanography, American Geophysical Union, 1-14. Bennett, A . F. , 1992: Inverse Methods in Physical oceanography. Cambridge University Press, 346pp. Bishop, C. M . , 1995: Neural network for pattern recognition. Clarendon Press . Oxford, 482pp. Bjerknes, J . , 1969: Atmospheric teleconnections from the equatorial Pacific. Mon. Wea. Rev., 97, 163-172. Blank, B . , J . D. Neelin, and D. Gutzler, 1997: Estimating the effect of stochastic wind stress forcing on ENSO irregularity, J. Climate, 10, 1473-1486. Chang, P., L. J i , B . Wang, and T. L i , 1995: Interaction between the seasonal cycle and E l Nino-Southern Oscillation in an Intermediate Coupled Ocean-Atmosphere Model, J. Atmos. Sci., 52, 2353-2372. —, —, H . L i , and M . Fliigel, 1996: Chaotic dynamics versus stochastic processes in E l Nino-Southern Oscillation in Coupled Ocean-Atmosphere Models, Physica D., 98, 301-320. Chen, D., M . Cane, E. Zebiak and A. Kaplan, 1998: The impact of sea level data assimilation on the Lamont model prediction of the 1997/98 E l Nino. Geophys. Res. Letter, 25, 2837-2840. Chapter 8. References 173 —, S. E . Zebiak and M . A. Cane, 1997: Initialization and Predictability of a coupled ENSO forecast model. Mon. Wea. Rev., 125, 773-788. —, —, A . J. Busalacchi and M . A . Cane, 1995: An improved procedure for E l Nino forecasting: Implications for predictability. Science, 269, 1699-1902. Courtier P., E . Anderson, W. Heckley, J. Pailleux, D. Vasiljevic, M . Hamrud, A . Hollingsworth, E . Rabier, and M . Fisher, 1998: The E C M W F implementation of three-dimensional variational assimilation (3D-VAR)-I-Formulation. Quart. J. Roy. Met. Soc, 124, 1783-1807. Daley, R., 1991: Atmospheric data analysis. Cambridge Univ. Press., 457pp. Davey, M . K . , S. Ineson, and M . A . Balmaseda, 1994: Simulation and hindcasts of tropical Pacific Ocean interannual variability. Tellus, 46A 433-447 Delecluse, P., M . K . Davey, Y . Kitamura, S. G. H . Philander, M . Suarez and L. Bengts-son, 1998: Coupled general circulation modelling of the tropical Pacific, J. Geophys. Res., 103, 14357-14373. Derber, J. and A . Rosati, 1989: A global oceanic data assimilation system. Phys. Oceanogr., 19, 1333-1347. Eckert, C. and M . Latif, 1997: Predictability of a stochastic forced hybrid coupled model of E l Nino. J. Climate , 10, 1488-1504. Eisner, J. B . , and A. A. Tsonis, 1992: Nonlinear prediction, chaos and noise. Bull. Amer. Meteor. Soc, 73, 49-60. Evensen, G., 1997: Advanced Data Assimilation for strongly nonlinear dynamics. Mon. Wea. Rev., 125, 1342-1354 Chapter 8. References 174 Fischer, M . , M . Latif, M . Fliigel and M . J i , 1997: The impact of data assimilation on ENSO simulations and predictions. Mon. Wea. Rev, 125, 819-829. Fliigel, M . and P. Chang, 1998: Does the predictability of ENSO depend on the seasonal cycle?. J. Atmos. Sci., 55, 3230-3243. Gauthier, P., 1992: Chaos and quadric-dimensional data assimilation: a study based on the Lorenz model. Tellus, 44A, 2-17. —, C. Charette, L. Fillion, P. Koclas and S. Laroche, 1999: Implementation of a 3D variational data assimilation system at the Canadian Meteorological Centre. Part I: The global analysis. Atmosphere-Ocean, 37, 103-156. Gelb, A . , 1974: Applied optimal estimation. The M . I. T. Press, 374pp. Ghil, M . and P. Malanotte-Rizzoli, 1991: Data Assimilation in Meteorology and Oceanog-raphy. Advances in Geophysics, 33, 141-265. Giering, R., and Kaminski,T., 1998: Recipies for Adjoint Code Construction, ACM TOMS, 24(4), 437-474. Gi l l ,A. E . , 1980: Some simple solutions for heat induced tropical circulation. Q. J. Roy. Met. Sci., 106, 447-462. —, 1982: Atmosphere-Ocean Dynamics. Academic Press, 662pp. GiU, P. E. , W. Murray and M . H. Wright, 1981: Practical Optimization, Academic Press, 401pp. Goldenberg, B . and J . J . O'Brien, 1981: Time and space variability of tropical Pacific wind stress, Mon. Wea. Rev., 109, 1190-1207 Chapter 8. References 175 Goswami, B . N . , and J . Shukla, 1991: Predictability of a coupled ocean-atmosphere model, J. Climate , 4, 3-22. Gu, D., and S. G. H . Philander, 1997: Interdecasdal climate fluctuation that depend on exchanges between the tropics and the extra-tropics. Science, 275, 805-807. Hannachi, A . , and B . Legras, 1995: Simulated annealing and weather regimes classifi-cation. Tellus, 47A, 955-973. Hao, Z., M . Ghil, 1994: Data assimilation in a simple tropical ocean model with wind stress errors. J. Phys. Oceanogr., 24, 2111-2128. Holton, J.R., 1992: An Introduction to Dynamical Meteorology. Academic Press, 511pp. Hirst, A . C., 1986: Unstable and damped equatorial models in simple coupled ocean -atmosphere models, J. Climate , 3, 1254-1281. Hoskins, B . J . , and R. P. Pearce, 1983: Large-scale dynamical processes in the atmo-sphere. Academic Press, 397pp. Houtekamer, P. L. , and H . L. Mitchell, 1998: Data assimilation using an ensemble Kalman filter technique. Mon. Wea. Rev., 126, 796-811. Hsieh, W. W, 2000: Nonlinear canonical correlation analysis by neural networks. Neural Networks, 13: 1095-1105. —, 2001: Nonlinear canonical correlation analysis of the tropical Pacific climate variabil-ity using a neural network approach. J. Climate, (accepted 2000/10/27). [preprint download from http://www.ocgy.ubc.ca/ william/pubs.html] —, and B . Tang, 1998: Applying neural network models to prediction and analysis in meteorology and oceanography. Bull. Amer. Meteor.Soc, 79, 1855-1870. Chapter 8. References 176 J i , M . , D. W. Behringer, and A. Leetmaa, 1998: An improved coupled model for ENSO prediction and implications for ocean initialization-PART II-The coupled model. Mon. Wea. Rev., 126, 1022-1034. —, A. Leetmaa and J . Derber, 1995: An Ocean Analysis for seasonal to interannual climate studies. Mon. Wea. Rev., 123, 460-481. —, A. Leetmaa and V . E . Kousky, 1996: Coupled model prediction of ENSO during the 1980s and 1990s at the National Centers for Environmental Prediction. J. Climate, 9, 3105-3120. Jin, F-F. , J . D. Neehn, and M . Ghil, 1998: E l Nino/Southern Oscillation and the annual cycle: subharmonic frequency-locking and aperiodicity. Physica D, 98, 442-465. —, — and —, 1994: ENSO on the devil staircase. Science , 264, 70-72. Kang, I. and J . Kug, 2000: An E l Nino prediction system using an intermediate ocean and a statistical atmosphere. Geophys. Res. Letter, 27, 1167-1170. Kirtman, B . , and P. S. Schopf, 1998: Decadal variabihty in ENSO predictabihty and prediction. J.Climate, 11, 2804-2822. —, and S. E. Zebiak, 1997: ENSO simulation and prediction with a hybrid coupled model. Mon. Wea. Rev., 125, 2620-2641. Kleeman, R., 1993: On the dependence of hindcast skiU in a coupled ocean-atmosphere model on ocean thermodynamics. J. Climate, 6, 2012-2033. —, A. Moore, and N . R. Smith, 1995: Assimilation of subsurface thermal data into a simple ocean model for the initiahzation of an intermediate tropical coupled ocean-atmosphere forecast model. Mon. Wea. Rev., 123, 3103-3113. Chapter 8. References 177 Kleeman, R., R. A . Colman, N . R. Smith, and S. B . Power, 1996: A recent change in the mean state of the Pacific basin climate: Observational evidence and atmospheric and oceanic responses. J. Geophys. Res., 101, 20483-20499. Kruger, J . , 1993: Simulated Annealing - A Tool for Data Assimilation into an Almost Steady Model State. J. Phys. Oceanogr., 23, 679-688 Latif, M . , D. Anderson, T. Barnett, M . Cane, R. Kleeman, A . Leetmaa, J . O'Brien, A . Rosati, and E. Schneider, 1998: A review of the predictability and prediction of ENSO. / . Geophys. Res., 103, 14375-14393. —, A . Sterl, E . Marier-Reimer, and M . Junge, 1993: Structure and predictability of the E l niho/Southern Oscillation Phenomenon in a Coupled Ocean-Atmosphere General Circulation Model. J. Climate, 6, 700-708. —, J . Biercamp, H . von Storch, M . J . McPhaden, and E. Kirk, 1990: Simulation of ENSO-related surface wind anomalies with an atmospheric G C M forced by observed SST. J. Climate, 3, 509-521. —, M . Fliigel, 1991: An Investigation of Short-Range climate predictability in the Tropical Pacific. J. Geophys. Res,96, 2661-2673. —, N . E . Graham, 1992: How much predictive skills is contained in the thermal structure of an Oceanic GCM? J. Phys. Oceanogr., 22, 951-962 —, A. Villwock, 1990: Interannual variability as simulated in coupled ocean-atmosphere model, J. Mar. Syst., 1, 51-60 —, R. Kleeman, and C. Eckert, 1997: Greenhouse warming, decadal variability or ENSO? A n attempt to understand the anomalous 1990s. J. Climate, 10, 2221-2239. Chapter 8. References 178 Le Dimet, F. , and 0. Talagrand, 1986: Variational algorithms for analysis and assimi-lation of meteorology observations: theoretical aspects. Tellus, 38A, 97-110. Lorenc, A . , 1981: A global three-dimensional multivariate statistical interpolation scheme. Mon. Wea. Rev., 109, 701-720. Lorenz, E. N . , 1963: Deterministic nonperiodic flow. J. Atmos. Sci., 20, 130-141. Lu, J . , and W. W. Hsieh, 1998: On determining initial conditions and parameters in a simple coupled model by adjoint data assimilation. Tellus, 50A, 534-544. Masters, T., 1995: Advanced Algorithms for neural networks—A C + + source book. John Wiley & Sons, Inc., 431pp. Miller, A . J . , T. P. Barnett, and N . E. Graham, 1993: A comparison of some tropical ocean models: Hindcast Skill and E l Nirio Evolution. J. Phys. Oceanoger., 23, 1557-1591. Miller, R. N . , M . Ghil and F. Gauthiez, 1994: Advanced data assimilation in strongly nonhnear dynamical systems. J. Atmos. Sci., 51, 1037-1056. —, and —, 1990: Data assimilation in strongly nonlinear current system. Proc. Int. Assimilation of Observation in Meteorology and Oceanography, pp 93-98, 9-13 July 1990, Clermont-Ferrand, France. Monahan, A . H . , 2000: Nonhnear principal component analysis by neural networks: theory and application to the Lorenz system. / . Climate, 13, 821-835. Moore, A . M . and D. L. T. Anderson, 1989: The assimilation of X B T data into a layer model of the tropical Pacific ocean. Dynamics of Atmospheres and Oceans, 13, 441-464. Chapter 8. References 179 Navon, I. M . , X . Zhou, J . Derber and J . Sela, 1992: Variational data assimilation with an adiabatic version of the N M C spectral model. Mon. Wea. Rev., 120, 1433-1446. —, 1998: Practical and theoretical aspects of adjoint parameter estimation and identi-fiability in meteorology and oceanography. Dynamics of Atmospheres and Oceans, 27 (1-4), 55-79. Neelin, J . D., 1990: A hybrid coupled general circulation model for E l Nino studies, J. Atmos. Sci., 5, 674-693. —, D. S. Battisti, A . C. Hirst, F-F. Jin, Y . Wakata, T. Yamagata and S. Zebiak, 1998: ENSO theory, J. Geophys. Res., 103, 14261-14287. —, and H . A . Dijkstra, 1995: Ocean-atmosphere interaction and the tropical chmatol-ogy. Part I: The dangers of flux correction. Climate Dynamics, 12, 101-112. Oberhuber, J . M , 1988: An atlas based on COADS dataset: the budgets of heat, buoy-ancy and turbulence kinetic energy at the surface of the global ocean, Max-Plank-Institute fur Meteorologie, Report No. 15 Palmer, T. N . , 1993: Extended-Range atmospheric prediction and the Lorenz model. Bull. Amer. Meteor. Soc, 74(1), 49-65. —, and D. L. T. Anderson, 1994: The prospects for seasonal forecasting -A review paper, Q. J. Roy. Met. Sci., 7, 755-793. Pares-Sierra, A. , and J . J . O'Brien, 1989: The seasonal and the interannual variability of the Cahfornia current system: A numerical model, J. Geophys. Res., 94, 3159-3180. Philander, S. G., 1990: E l Nino, La Nina, and the Southern OsciUation. Academic Press, 293pp. Chapter 8. References 180 —, R. C. Pacanowski, N .C . Lau and M . J . Nath, 1992: Simulation of ENSO with a global atmospheric G C M . J. Climate, 5, 308-329 Press, W. H . , S. A . Teukolsky, W. T. Vetterling, and B . P. Flannery: 1992: Numerical recipes in Fortran, Second edition, Cambridge University Press, 963pp. Rasmusson, E. M . , and T. H . Carpenter, 1982: Variations in tropical sea surface tem-perature and surface wind fields associated with the Southern Oscillation/El Nifio. Mon. Wea. Rev., 110, 354-384. Rosati,A, K . Miyakoda and R. Gudgel, 1997: The impact of ocean initial conditions on ENSO forecasting with a coupled model. Mon. Weather Rev., 125, 754-772. Saunders, A . , M . Chil and D. Neelin, 1998: Forecasts of Nifio 3 SST anomalies and SOI based on singular spectrum analysis combined with the maximum entropy method. Forecast Bulletin, 7, 41-43. Schopf, P. S, and M . J . Suarez, 1988: Vacillations in a coupled ocean-atmosphere model. J. Atmos. Set., 45, 549-566 Smith, T. M . , R. W. Reynolds, R. E. Livezey, and D. C. Stokes 1996: Reconstruction of Historical Sea Surface Temperatures Using Empirical Orthogonal Functions. J. Climate, 9, 1403-1420. Slutz, R. et al 1985: The Comprehensive ocean-atmosphere data set release.I: Climate Research Program, E R L / N O A A , Boulder,CO, 39pp. Suarez, M . J . , and P. S. Schopf, 1988: A delayed action oscillator for ENSO. Atmos. Sci., 45, 3283-3287. Chapter 8. References 181 Syu, H . H . , J . D. Neelin, and D. Gutzler, 1995: Seasonal and Interannual variability in a hybrid coupled G C M , J. Climate, 9, 2121-2143. —, and —, 1998: Simulation and prediction of E l Nino and the variation of ENSO phase locking, Ninth Conference on Interaction ofthe sea and atmosphere, Phoenix, 11-16 January 1998, 115-119. — and —, 2000: ENSO in a hybrid coupled model. Part II: Prediction with piggyback data assimilation. Clime Dyn., 16, 35-48. Talagrand, 0., 1991: The use of adjoint equations in numerical modelling of the at-mospheric circulation. In Andreas Griewank and George F. Corliess, editors, Au-. tomatic Differentiation of Algorithms: Theory, Implementation and Application. 169-180, SIAM. Tang, B . , W. W. Hsieh, A . H. Monahan, and F. T. Tangang, 2000: Skill compari-sions between neural network and canonical correlation analysis in predicting the equatorial Pacific sea surface temperatures, J. Climate, 13, 287-293. Tang, Y , 2001: Hybrid coupled models of the tropical Pacific - Interannual variabihty, submitted to Climate Dynamics. —, and Hsieh, W, 2001d: An ENSO operational forecast model with data assimilation, submitted to Experimental Long-Lead Forecast Bulletins. —, and —, 2001c: Impact of data assimilation on ENSO simulations and predictions, submitted to J of Climate. —, and —, 2001b: Hybrid coupled models of the tropical Pacific — ENSO predictability, submitted to Climate Dynamics. Chapter 8. References 182 — and —, 2001a: Coupling neural networks to incomplete dynamical systems via vari-ational data assimilation, Mon. Wea. Rev. 129, 818-834. Tang, Y . , and W. W. Hsieh, B. Tang and K . Haines, 2001: A neural network atmospheric model for hybrid coupled modelling, Climate Dynamics, 17, 445-455. Tangang, F . T., W. W. Hsieh, and B . Tang, 1998a: Forecasting regional sea surface temperatures of the tropic by neural network models, with wind stress and sea level pressure as predictors. J. Geophys. Res., 103, 7511-7522 , B . Tang, and A . H . Monahan, and W.W. Hsieh, 1998b: Forecasting ENSO events: a neural network-extended E O F approach. J. Climate, 11, 29-41 Trenberth, K . , W. Branstator, D. Karoly, A. Kumar, N . Lau and C. Ropelewski, 1998: Progress during T O G A in understanding and modelling global teleconnection asso-ciated with tropical sea surface temperatures. / . Geophys. Res., 103, 14291-14324. —, and T.J Hoar, 1997: E l niho and climate change. Geophys. Res. Lett., 24, 3057-3060. Thacker, W. C. and R. B . Long, 1988: Fitting dynamics to data. J. Geophys. Res., 93, 1227-1240. Tziperman, E . , M . A. Cane, and S. E. Zebiak, 1995: Irregularity and Locking to the seasonal cycle in an ENSO prediction model as explained by the Quasi-periodicity route to chaos, J. Atmos. Sci.,52, 293-306 —, and W. C. Thacker, 1989: An optimal control/adjoint equations approach to study-ing the oceanic general circulation. J. Phys. Oceanogr., 19, 1471-1485. Chapter 8. References 183 von Storch, H . , and F. W. Zwiers, 1999: Statistical Analysis in Climate Research, Cam-bridge, Cambridge Univ. Pr., 484pp. Wallace,J., E . Rasmusson, T. Mitchell, V.Kousky, E. Sarachik and von Storch, 1998: On the structure and evolution of ENSO-related climate variability in the tropical Pacific: Lessons from T O G A . J. Geophys. Res., 103, 14241-14259. Wang, B . , 1995: Interdecadal changes in E l Nifio onset in the last four decadals. J. Climate, 8, 267-285. White, W. B . , and S. E. Pazan, 1987: Hindcast/forecast of ENSO events based upon the redistribution of observed and model heat content in the western tropical Pacific, 1964-86. J. Phys. Oceanoger., 17, 264-280 Wu, D-H. , D. L. T. Anderson, and M . K . Davey, 1994: ENSO prediction experiments using a simple ocean-atmosphere model. Tellus, 49A, 464-480 Xue, Y . , M . A. Cane, S. E. Zebiak and M . B . Blumenthal, 1994: On the prediction of ENSO: a study with a low order markov model. Tellus, 4@A, 512-528 Zhang, Y . , J . M . Wallace, and D. S. Battisti, 1997: ENSO-like interdecadal variability. J. Climate, 2, 1004-1020. Zebiak, S. E . , 1985: Tropical atmosphere-ocean interaction and the E l Nino-Southern Oscillation phenomenon. Ph.D thesis, M.I.T., 261pp. —, 1986: Atmospheric convergence feedback in a simple model for E l Nino. Mon. Wea. Rev., 114, 1263-1271. —, 1989: Oceanic heat content variability and E l Nino cycles. / . Phys. Oceanogr, 19, 475-486 Chapter 8. References 184 —, and Cane, M . , 1987: A model for E l Nino-Southern Oscillation. Mon. Wea. Rev., 115, 2262-2278 Zhu, Y . , and I. M . Navon, 1999: Impact of Key Parameters Estimation on the perfor-mance of the F S U spectral model using the full physics adjoint. Mon. Wea. Rev., 127, 1497-1517. Appendix A Ocean model formulation Like the two-layer model of Balmaseda et al (1994), this 6-layer model includes explicit nonlinear thermodynamics for all layers, accounting for horizontal advection, vertical heat transport, diffusion and, in the first layer, surface heat flux. Besides the physical processes considered in the two-layer model, the 6-layer model includes a viscosity between these layers to account for wind driven Ekman effects, simple convective mixing between two adjacent layers due to upwelling and horizontal advection. The model equations are as follows: — + ui— + vi— + / k A him = {-VPi + M i + — + ^ V 2 ( ^ U i ) + V i , ( A . l ) ot dx dy J po po (wt+Uii+v%)h< = D > + d V % ' <A2» where layer subscript i = 1,...,N. The symbols u and v are the zonal and meridional components of velocity respectively, p is the pressure, h is the layer thickness, r is the wind stress, T is temperature, / is the Coriolis parameter given by / = By, p0 is a constant density, and u, d and n are the coefficients of momentum, thickness and heat diffusivity. V 2 is horizontal Laplacian. The terms D{,Mi,Hi are the mass, momentum and heat exchanges between layers and are discussed below. The term Qa is the surface heat flux determined as for the 2-layer model and V i is a vertical viscosity used if the Ekman layer is resolved. V i is parameterized as 185 Appendix A. Ocean model formulation 186 V i = Aj(ui_i - U i ) + Ai+1(ui+1 - ui), where Ai+l < A{ The pressure gradients are ^VJH = gaV^hj(Ti-TN+l) + ^ - g a ^ + f^h^JVTt (A.4) where g and a are the gravitational constant and the thermal expansion coefficients respectively. The pressure at the sea surface is ( N \ N Ho ~ zZ hj] -9J2 Pihh (A-5) J'=I / i=i and Ho is the total depth of the model (a constant), and pj is the density of the j t h layer. Di is based on the upwelhng, ee,-, and downwelling, e^, velocities between layers i and i + 1 which are parameterised as follows: eei = (Hi - hi)2/teiHi, if h{ < Hi else eei = 0, (A-6) tdi = — {Hi — hi)2/tdiHi, if /ii > Hi else = 0, where Hi is the nominal layer depth towards which entrainment or detrainment tend to force layer i, and tei and t<n are the entrainment and detrainment time scales respectively. The total vertical velocity may therefore be expressed as e'i = 2^ e " ^ l e " ^ 2^€di ~ The thickness of the ith layer can be changed by mass exchanges at both its top and bottom surfaces, with layer i — 1 above and i + 1 below. The total mass exchanged in eq. A.2 is given by Di = e'i-Appendix A. Ocean model formulation 187 Table A . l : Values of the parameters used in the 6-layer ocean model Layer H(m) todays) id(days) Hinit(m) Tinit(°C) A 1 100 1 500 100 26 0.75 2 175 150 500*150 175 16 1 3 250 150 ' 500*150 250 13 1 4 320 150 500*150 320 10 1 5 400 150 500*150 400 8 1 6 500 150 500*150 500 6 1 Note that within each layer the parameters only control entrainment or detrainment with the layer below. Any entrainment or detrainment with the layer above is determined by the thickness of the layer above. The momentum exchanges for layer i are M i = <.ui+1 - c ; _ l U , for e\ > 0, > 0, (A.7) M i = e> i + 1 - e U u ^ for e\ > 0, t[_x < 0, M i = e > - e^Ui -x for <• < 0, e\_x < 0, M i = ejm - e ^ U i for e[ < 0, e\_x > 0. The heat transport term related to upwelling and downwelling is Hi = -e'MTi - Ti+1)/hi - eUAi-iCZU - Ti)/hi for t\ > 0, < 0, (A.8) Hi = 0 for e[ < 0, > 0, where \ is a parameter controlling the effective vertical temperature gradient, taken to be 1. Table A . l lists the parameters used in the model. Appendix B Discrete form of Lorenz hybrid model for coding The computer code is directly generated from the discrete form of the hybrid model. The 4th- order Runge-Kutta is used to discretize the vector equation (6.7): , k i n . k 2 n k 3 n k 4 n . . v n + 1 = v n + - ^ - r - r + ^ - + - r , (B. l) k l n = fcf(vn>P>tn), (B.2) k 2 n = W(v„ + 0 .5k l n > p ,< n + 0.5J»), (B.3) k 3 n = hf(vn + 0.5k2 n, p, tn + 0.5h), (B.4) k 4 n = fcf(vn + k 3 n > p , *„ + /»)> (B.5) where the subscript n indicates the time level, f = ( / l , / 2 , / 3 ) t p , and f l , f2 and f3 are scalar functions: fl = -aX + aY, (B.6) fl = -XZ + bX - Y, (B.7) /3 = Wj[tanh{WljX + w2jY + w3jZ + bj)} + b, (B.8) j j=l,2. . M (M=5, the number of hidden neurons). Eq. B . l can be written as v n + i = v n + K n (B.9) The Lagrange function of Eq. (6.10) is discretized as N-l L = J+Y1 < t r (v n + 1 - v n - K n ) , (B.10) n=0 188 Appendix B. Discrete form of Lorenz hybrid model for coding 189 J = E(v» - ) t r W - 1 ( v „ - ), (B.ll) n=l where v* is the vector of Lagrange multiphers and its dimension is the same that of v . The equations for the best fit are obtained by requiring that all derivatives of L vanish: OL dL = 0, (B.12) = 0, (B.13) | = »• <B-"> where the partial derivative with respect to a vector indicates a partial derivative with respect to each ofthe vector components. Eq. (B.12) gives nothing more than the (B.9); whereas (B.10) gives ! = t 1 ( - l f K ' < B - i 5 > The gradients from (B.15) are then provided to a quasi-Newton method to seek the optimal parameters. The equations describing the trajectory ofthe Lagrange multiphers v*, i.e., the adjoint equations, can be obtained from d L/d v n and Eq. (B.10) (also see Anderson and Willebrand, 1989): 2 £ W - 1 ^ - v f ) + - < - ( ^ ) v ; = 0, n = 1 , N , (B.16) n=l a V " w. = ( B 1 7 ) = 0. (B.18) The SKn/<9vn and <9Kn/<9p can be obtained by chain-rule, e.g. dkn ldkln l d k 2 n < 9 k l n 1 <9k3„ dk2n dkln 1 <9k4n <9k3n dk2n dkln 5vn 6 <9vn 3dkln <9vn 3<9k2„<9kl n <9v„ 6 dk3n dk2n dkln <9vn Appendix B. Discrete form of Lorenz hybrid model for coding 190 Thus, the cycling procedure for computing the best fit of the hybrid model (6.4)-(6.6) to Lorenz model (6.l)-(6.3) is: (i) Given initial guess values for the unknown parameters, using (B.9) to step the hybrid model forward in time from 1 to N to compute a first approximation to the Lorenz model. (ii) Step equations (B.16) - (B.18) backward in time from N to 1 to compute approximate values of the Lagrange multipliers. (iii) Evaluate the gradient of the cost function with respect to these parameters using Eq. (B.15)) (and/or Eq. (B.17) if initial conditions need to be retrieved). (iv) Use a descent algorithm to find the minimum of the cost function in the direction opposite to the gradient by a hne search. (v) Take results of the line search as the next guesses for the parameters and repeat the process starting from (i) until convergence criteria are satisfied. The quasi-Newton method of Broyden-Fletcher-Goldfarb-Shanno (BFGS) (Press et al. 1992) is used in the above procedure for searching the optimal solutions. Fig. B . l is the variations of the normalized cost function and normahzed gradient with the number of iterations for a weakly nonlinear case (Fig. 6.7) and a strongly nonlinear case (Fig. 6.11). The BFGS algorithm has a memory requirement of 0 (N 2 ) where N is the number of parameters involved in the optimization. In our case, the memory is 183 K B . For problems with larger N , the conjugate gradient method with a memory requirement of O(N) could be more suitable. A typical case here requires about 50 iteration for convergence, using about 3 minutes of C P U time on an SGI Origin 200 server. In practice, the adjoint codes were not developed from Eqs. (B.16)-(B.18), but by the T A M C (Tangent and Adjoint Model Compiler, Giering and Kaminski, 1996) via the tangent hnear equations. These equations are derived from Eq. (B.9): o v n + 1 = A n J v n (B.19) Appendix B. Discrete form of Lorenz hybrid model for coding 191 5 10 15 20 25 30 10 20 30 40 Number of iterations Number of iterations F i g u r e B . l : V a r i a t i o n s of the no rma l i zed cost func t ion (J/ Jo) (sohd hne) a n d n o r m a l i z e d gradient (||g|| ||go||_1) (dash hne) w i t h the number of i te ra t ions for (a) the w e a k l y non-l inear case of F i g . 6.7 and (b) the s t rongly nonhnear case of F i g . 6.11. T h e be t te r i n i t i a l guesses i n (a) meant tha t the cost func t ion and its gradients h a d re l a t ive ly less d is tance to d rop to convergence t h a n i n (b). Appendix B. Discrete form of Lorenz hybrid model for coding 192 where the m a t r i x A n is the tangent l inear m o d e l at t ime-step n , and A n = I + ^ p (B.20) w i t h I the the i den t i t y m a t r i x . T h e adjointness has been checked us ing the r e l a t i on be tween the tangent- l inear code A and its adjoint code A* (Ta l ag rand 1991; N a v o n et a l . 1992): {a, Ab} = {A*a,b} (B.21) where {,} denotes the inner p roduc t , and a and b a rb i t r a ry vectors . Append ix C Summary of Exper iments Table C l : Summary of main experiments and results Experiments Results Locations Ocean model control run Fig 2.2 and Table 2.1 Chap. 2.2.2 Choice of number of hidden neuron 3 Chap. 2.3.2 Construction of two atmospheric models Table 2.3 Chap. 2.4.1 Response of ocean to two atmos. models Figs. 2.11-2.16 Chap. 2.5 L H C M and N H C M standard coupling runs (5=1) Figs. 3.2-3.12 Chap. 3.3 L H C M and N H C M coupling runs (£=0.83) damped oscillation Chap. 3.4 L H C M and N H C M coupling runs (<f=1.5 and 1.68) Fig . 3.13 Chap. 3.4 Predictions with 'stress-only' initialization Figs 4.1-4.5 Chap. 4.3 Predictions with 'stress-from-ocean' initialization Figs L 1.12-4.16 Chap. 4.4 H C A assimilation experiment for ENSO simulation Figs 5.5-5.7 Chap. 5.4 SSHA assimilation experiment for ENSO simulation Figs 5.5-5.7 Chap. 5.4 SST assimilation experiment for ENSO simulation Figs 5.5-5.7 Chap. 5.4 Wind assimilation experiment for ENSO simulation Figs 5.5-5.7 Chap. 5.4 H C A assimilation experiment for ENSO prediction Fig 5.8 Figs 5.12-5.14 Chap. 5.5 SSHA assimilation experiment for ENSO prediction Fig 5.8 Figs 5.12-5.14 Chap. 5.5 SST assimilation experiment for ENSO prediction Fig 5.8 Figs 5.12-5.14 Chap. 5.5 Wind assimilation experiment for ENSO prediction Fig 5.8 Figs 5.12-5.14 Chap. 5.5 Simple Lorenz hybrid model experiments Figs. 6.2-6.3 Chap. 6.4 Initial guess experiments for 4D Var assimilation Figs. 6.5-6.7 Chap. 6.5.1 Assimilation window experiments Fig. 6.8 Chap. 6.5.2 SC assimilation experiments (100 cases) Figs. 6.9-6.13 Chap. 6.5.3 Experiment for determining all parameters Fig. i 5.16-6.17 Chap. 6.6 193 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0053123/manifest

Comment

Related Items