This paper introduces the advanced MP-PIC-VOF model tailored for dense particle-laden flows with free surface, which has been developed and extensively tested across a set of validation cases found in literature and original bulk particle water entry case. A distinctive feature of the MP-PIC method is its demonstrated ability to accurately capture the behavior of closely packed particles in a fluid, even in the absence of direct pairwise particle-particle interactions. At a closed packed limit, the MP-PIC method achieves the accurate representation of the state through the resolved mean particle velocity field and implementation of the velocity limiter in the inter-particle stress force. The new model integrates a trilinear interpolation technique, specifically adapted for unstructured hexahedral meshes, and a weighted least squares method for efficient gradient computation that operates at a sub-cell level, enabling more accurate calculation of inter-particle stress gradients. Other key contributions include the integration of hydrostatic pressure adaptation in the momentum equation and a volume-conservative alpha transport equation that ensures mass conservation during the transfer of the solid phase between distinct fluid phases. The coupling framework includes a range of coupled fluid-particle forces important for particles immersed in liquid, including a dense virtual mass force. The model's validation against experimental data and CFD-DEM-VOF results focuses on key flow parameters, specifically particle velocity, dispersion profile, and cavity evolution during bulk particle water entry. The model is shown to accurately simulate complex solid-liquid-gas interactions, demonstrating its potential for optimizing a wide range of complex industrial processes such as liquid fluidized beds, solid-liquid stirred tanks, and clarifiers.