Abstract. A numerical model, ISWFoam, for simulating internal solitary waves (ISWs) in continuously stratified, incompressible, viscous fluids is developed based on a fully three-dimensional (3D) Navier-Stokes equation using the open source code OpenFOAM. This model combines the density transport equation with the Reynolds-averaged Navier-Stokes equation with the Coriolis force, and the model discrete equation adopts the finite volume method. The k-ω SST turbulence model has also been modified accordingly to the variable density field. ISWFoam provides two initial wave generation methods to generate an ISW in continuously stratified fluids, including solving the weakly nonlinear models of the extended Korteweg–de Vries (eKdV) equation and the fully nonlinear models of the Dubreil-Jacotin-Long (DJL) equation. Grid independence tests for ISWFoam are performed, considering the accuracy and computing efficiency, the appropriate grid size of the ISW simulation is recommended to be one-one hundred and fiftieth of the characteristic length and one-twenty fifth of the ISW amplitude. Model verifications are conducted through comparisons between the simulated and experimental data for ISW propagation examples over a flat bottom section, including laboratory scale and actual ocean scale, a submerged triangular ridge, a Gaussian ridge and slope. The laboratory test results, including the ISW profile, wave breaking location, ISW arrival time, and the spatial and temporal changes in the mixture region, are well reproduced by ISWFoam. The ISWFoam model with unstructured grids and local mesh refinement can accurately simulate the generation and evolution of ISWs, the ISW breaking phenomenon and the interaction between ISWs and complex structures and topography.