Large scale engineering as well as natural systems, such as weather, often have high-dimensional state spaces and exhibit chaotic dynamics. To model the behavior of such systems, sets of coupled Partial Differential Equations (PDEs) are formulated and solved using high-performance computing systems. More recently, significant attention has grown toward the use of Artificial Intelligence (AI) and Machine Learning (ML) techniques, in particular, using data-driven modeling to learn fast and accurate surrogate process models trained on high-resolution data obtained from simulations, or observations of chaotic systems. Echo state networks (ESN), a family of recurrent neural network algorithms, have emerged as one of the most promising techniques to learn predictive models of chaotic dynamical systems directly from data. In spite of their success in learning chaotic dynamical systems from data, there are many open questions. Some of them are practical engineering concerns such as: how to choose training parameters (reservoir size, spectral radius, length of training sequence) for specific problems, how robust the learned models are to variations in data, and in training parameters (initialization of random weights, reservoir size, spectral radius). Others are open theoretical questions such as: why do ESNs work at all, in particular, which aspects of the underlying dynamical systems are captured by the learned reservoirs, and which factors determine the prediction horizon of the learned models. In this thesis, we study these practical and theoretical questions in the context of two models of chaotic dynamical systems, Lorenz63 and Lorenz96, which are prototypes of more complex weather models. We show that the predictive performance of the learned models is highly sensitive to initial conditions — i.e., for different training sequences all of the same lengths but with different initial states, there is considerable variation in prediction horizon from 0.1 MTU to 3.8 MTU in Lorenz63 and from 0.4 MTU to 2.8 MTU in Lorenz96. We also show that variations in the initialization of (random) input weights and (random) reservoir weights at the start of the training phase yields models with varying prediction horizon for the very same training sequence. We discuss the implications of these findings in the construction of robust ESN models for Lorenz systems. To help explain the observed variations in predictive performance with initial conditions, and to understand when and why ESNs work, we use dimensionality reduction and clustering algorithms to visualize the evolution of high-dimensional reservoir states during training and prediction. Our main finding is that, in a well-trained model, reservoir states mirror the dynamics of the chaotic system from which the data is derived. In particular, we can infer the number of dynamical components from the non-linear clustering of the reservoir states. In the context of Lorenz63, we show that the sensitivity to initial conditions stems from the locations of the initial condition relative to the two components of the underlying system.