The use of computational fluid dynamics (CFD) has become essential for aerospace design optimization processes. The computational cost of high-fidelity CFD is often very large and can make design optimization prohibitively expensive if a large number of design evaluations are required. Reduced-order models (ROMs) are a method that can be used to mitigate this cost. ROMs are low-dimensional data-driven surrogate models that are trained using a set of computed high-fidelity simulation snapshots. Many ROMs utilize the proper orthogonal decomposition (POD), a linear subspace method for representing solution spaces. While ROMs are becoming increasingly popular, they do face some challenges in their practical use, which include maximizing accuracy for a given computational budget, the ability to generalize throughout a parameter space, and applicability to topologically dissimilar meshes. In this dissertation, algorithms are introduced to improve the performance, stability, and understanding of data-driven surrogate CFD models and their applications. As ROMs tend to use a small amount of training data, their predictive performance is highly sensitive to their choice. Algorithms to improve the data selection process for POD-based ROMs are introduced using Isomap, a versatile technique for nonlinear dimensionality reduction, resulting in significantly improved predictive performance for a given computational budget when used over a traditional and widely used statistical sampling technique. Next, ROMs using artificial neural networks, specifically convolutional autoencoders (CAEs), are introduced to address the performance limits of POD for problems that are highly nonlinear or require large amounts of training data, such as unsteady ROMs involving multiple designs. A steady ROM framework combining CAEs with Gaussian process regression (GPR) is introduced and shown to significantly outperform POD when applied to a highly nonlinear lid-driven cavity problem. Ensemble learning is also used to effectively address the issue of error propagation in unsteady ROMs, where errors made early on can accumulate and lead to large inaccuracies over long time horizons at unseen design points. Finally, field inversion and machine learning (FIML) is proposed as an an alternative to ROMs for problems that require topologically dissimilar meshes. Field inversion involves obtaining improvements to turbulence models by augmenting them with a corrective field that is obtained using gradient-based optimization. Using a machine learning model trained on local flow variables and their gradients, a data-driven turbulence model is introduced to improve the predictive capabilities of baseline turbulence models, allowing for the prediction of complex flow phenomena present in experimental results.