A granular landslide impacting a river may lead to the formation of a landslide dam blocking the streamflow, and subsequently create a barrier lake. Should a barrier lake outburst, the flood may be destructive and spell disastrous consequences downstream. The last decade or so has witnessed a number of experimental and numerical investigations on barrier lake outburst flooding, whilst studies on barrier lake formation remain rare – a physically enhanced and practically viable mathematical model is still missing. Generally, barrier lake formation is characterized by multi-physical, interactive processes between water flow, multi-sized sediment transport and morphological evolution. Here, a new double layer-averaged two-phase flow model is proposed, which is an advance on existing continuum models that involve a single-phase flow assumption and presume a single sediment size, and discrete models that preclude fine grains and assume narrow grain size distributions. The proposed model is first validated against data from previous laboratory experiments of waves due to landslides impacting reservoirs and landslide dam formation over dry valleys. Then it is applied to explore the complicated mechanism and threshold for barrier lake formation. The water and grain velocities are shown to be disparate, characterizing the primary role of grains in driving water movement during subaqueous landslide motion and also demonstrating the need for a two-phase flow approach. The grain size effects are revealed, i.e., coarse grains and grain-size uniformity favour barrier lake formation. A new threshold condition is proposed for barrier lake formation, integrating the landslide-to-river momentum ratio and grain size effects. The present work facilitates a promising modelling framework for solving barrier lake formation, thereby underpinning the assessment of flood hazards due to barrier lakes.