Generating a digital twin of any complex system requires modeling and computational approaches that are efficient, accurate, and modular. Traditional reduced order modeling techniques are targeted at only the first two but the novel non-intrusive approach presented in this study is an attempt at taking all three into account effectively compared to their traditional counterparts. Based on dimensionality reduction using proper orthogonal decomposition (POD), we introduce a long short-term memory (LSTM) neural network architecture together with a principal interval decomposition (PID) framework as an enabler to account for localized modal deformation, which is a key element in accurate reduced order modeling of convective flows. Our applications for convection dominated systems governed by Burgers, Navier-Stokes, and Boussinesq equations demonstrate that the proposed approach yields significantly more accurate predictions than the POD-Galerkin method, and could be a key enabler towards near real-time predictions of unsteady flows. and perhaps much beyond. With the recent wave of digitization, reduced order modeling can be viewed as one of the key enablers to bring the promise of the digital twinning concept closer to reality [38]. Therefore, there is a continuous demand for the development of accurate reduced order models for complex physical phenomena. In projection-based ROMs, the most widely used technique, the discrete high-dimensional operators are projected onto a lower-dimensional space, so that the problem can be solved more efficiently in this reduced space [39][40][41][42][43][44].One of the very early-developed and well-known approaches to build this reduced space is Fourier analysis. However, it assumes universal basis functions (or modes) which have no specific relation to the physical system. On the other hand, snapshot-based model reduction techniques tailor a reduced space that best fits the problem by extracting the underlying coherent structures that controls the major dynamical evolution we are interested in. Proper orthogonal decomposition (POD) is a very popular and well-established approach extracting the modes which most contributes to the total variance [45,46]. In fluid dynamics applications, where we are mostly interested in the velocity field, those modes contain the largest amount of kinetic energy [47,48]. That is why POD is usually classified as an energy-based decomposition method. Another popular approach for model order reduction is the dynamic mode decomposition (DMD) [49][50][51][52][53][54] which generates a number of modes, each characterized by an oscillating frequency and growth/decay rate. In the present study, we are interested in the application of POD for dimensionality reduction.POD generates a set of spatial orthonormal basis functions, each containing a significant amount of total energy. To obtain a reduced representation of a system, the first few modes are selected, and the remaining are truncated assuming their contribution to the system's behavior is minimum (i.e.,...