This thesis aims to improve the understanding of El Niño Southern Oscillation (ENSO) diversity, in its future change, modeling and predictability. How might ENSO change in the warming climate? To reach a comprehensive understanding, a set of empirical probabilistic diagnoses (EPD) is introduced to measure the ENSO behaviors as to tropical Pacific sea surface temperature (SST) climatology, annual cycle, ENSO amplitude, seasonal phase locking, diversity in peak location and propagation direction, as well as the El Niño-La Niña asymmetry in amplitude, duration and transition. This diagnosis is applied to the observations, and consistency with previous studies indicates it is valid. Analysis of 37 CMIP5 model simulations for the 20th century and the 21st century shows that, other than the projected increase in SST climatology, changes in other aspects are largely model dependent and generally within the range of natural variation. The change favoring eastward propagating El Niños is the most robust seen in the SST anomaly field. To what extent can we trust the future projection? CMIP5 models show large spreads in terms of 20th century ENSO performance. So a data-driven approach called Empirical Model Reduction (EMR) is carried out, by fitting a low-dimensional nonlinear model from the observation with a representation of memory effect and seasonality. The stochastic simulation of EMR is able to reproduce a realistic ENSO diversity statistics and a reasonable range of natural variation, thus provides an additional benchmark to evaluate the CMIP5 model biases. What are the key model components leading to a good performance to simulate and predict ENSO? Using a suite of models under the aforementioned framework of EMR, control experiments are conducted to advance the understanding of ENSO diversity, nonlinearity, seasonality and the memory effects. Nonlinearity is found necessary to reproduce the ENSO diversity features by simulating the extreme El Niños. Nonlinear models reconstruct the skewed distribution of SST anomalies and improve the prediction of the El Niño-La Niña transition. Models with periodic terms reproduce the SST seasonal phase locking but do not improve the prediction appreciably. Models representing the ENSO memory effect, based on either the recharge oscillator (multivariate model with tropical Pacific subsurface information) or the time-delayed oscillator (multilevel model with SST history information), both improve the prediction skill dramatically. Models with multiple ingredients capture several ENSO characteristics simultaneously and exhibit overall better prediction skill. In particular, models with a memory effect show an alleviated skill drop during the spring barrier and a reduced prediction timing delay. One new ENSO prediction target is to predict not only the occurrence and amplitude of El Niño (EN) but also the peak location is at the central Pacific (CP) or the eastern Pacific (EP). Many prediction models have difficulty with it, which motivates the investigation on whether such ENSO diversity has intrinsically limited predictability. Here three aspects are addressed including the source/limit of predictability, time range and uncertainty. Approaches are combined including linear inverse modeling, singular vector analysis and probabilistic measure. The results show that two similar initial conditions with western Pacific SST warming anomalies may finally develop to either CPEN or EPEN. The equatorial Pacific subsurface evolution is important to tell the final outcome. Restricted by the chaotic property, the prediction horizon appears to be ~4 months before CPEN and ~7 months before EPEN. A flavor prediction model using data's transition probabilities is introduced as a new benchmark for probabilistic prediction.