Machine learning (ML) is experiencing an immensely fascinating resurgence in a wide variety of fields. However, applying such powerful ML to construct subgrid interphase closures has been rarely reported. To this end, we develop two data-driven ML strategies (i.e., artificial neural networks and eXtreme gradient boosting) to accurately predict filtered subgrid drag corrections using big data from highly resolved simulations of gas-particle fluidization. Quantitative assessments of effects of various subgrid input markers on training prediction outputs are performed and three-marker choice is demonstrated to be the optimal one for predicting the unseen test set. We then develop a parallel data loader to integrate this predictive ML model into a computational fluid dynamic (CFD) framework. Subsequent coarse-grid simulations agree fairly well with experiments regarding the underlying hydrodynamics in bubbling and turbulent fluidized beds. The present ML approach provides easily extended ways to facilitate the development of predictive models for multiphase flows. K E Y W O R D S coarse-grid simulation, filtered two-fluid model, fine-grid simulation, fluidization, machine learning, subgrid closure model 1 | INTRODUCTION Gas-solid fluidized bed reactors (FBRs) are ubiquitous in process engineering. Unfortunately, inhomogeneous mesoscale flow structures in such reactors are incredibly complex and manifest persistent instabilities in flow properties spanning all sorts of spatiotemporal scales. 1-3 This dynamic nature is of fundamental importance in essential process features such as the interphase momentum, species, and heat transfer. 4 Development of accurate and reliable modeling methods for capturing such complex dynamics has been the grand challenge for FBR modelers. One of the broadly-practiced predictive methods is the kinetic theory of granular flow (KTGF) based two-fluid model (TFM), which treats the fluid and particle phases as interpenetrating continuums. 5 However, the TFM simulation with sufficiently fine grids is necessitated for resolving all relevant scales of flow hydrodynamics. Fullmer and Hrenya 3 indicated that a typical mesh spacing of~3.3d s is required to obtain the grid independency for simulation statistics at intermediate solid phase fraction, and the grid size less than 10d s for dilute clustering flows is demanded for computational accuracy. For dense bubbling FBRs, Wang et al. 6 recommended a grid resolution of 2-4 d s . Such high requirements (10 −5 -10 −3 m) for reactor-scale simulations (10 −1 -10 1 m) are computationally prohibitive.Alternatively, subgrid models (SGMs) for interphase closure corrections such as the energy minimization multi-scale (EMMS), [7][8][9] filtered TFM (fTFM), 2,10-15 spatially-averaged TFM (SA-TFM), 16-18 and gradient-based method, [19][20][21] have been developed to account for effects of unresolved inhomogeneous structures. SGM enables a dramatic reduction of computational overhead by employing coarse grids while still maintaining high accuracy. In the establis...