Abstract. The first aim of this work is to improve the models currently available for the simulation of ice production in turbulent seawater, by means of the development of a multiphase model able to describe all the stages of ice production, overcoming the limitation of previous attempts, mainly based on Boussinesq approximation. We consider the mixture of ice and seawater as a dense compressible fluid, and we model the behaviour of seawater by an equation of state that links seawater density to temperature, salinity and pressure. The model is able to reproduce the interaction phenomena occurring between phases when the ice volume fraction exceeds the values allowed by the Boussinesq approximation, including in the momentum equations additional terms, related to the drag force between liquid and particles, and to the particle-particle interaction force. The second aim of our work is to implement and validate a numerical solver of our model. For this purpose, the model uses a sophisticated modelling approach, typically adopted for the numerical simulation of multiphase flows of industrial interest. The multiphase model can be coupled to the Large Eddy Simulation technique. The behaviour of the governing equations in the incompressible limit is investigated by means of a low-Mach number asymptotic analysis. The divergence constraint condition for the velocity field of continuous phase can then been imposed on the zero-Mach number equations by means of a projection method. The governing equations are discretized using the finite volume method, and the performance of the multiphase model has been assessed by solving a laminar Rayleigh-Bénard convection, for both large and small ice concentration regimes. In small concentration regime, the numerical solutions have been compared with the solutions obtained by a finite difference numerical code, based on the Boussinesq approximation.