You know that ice is less dense than water. Then, water that freezes will stay at the surface. Also, take in mind that water will only freeze on the surface. Then as you said, any ice that forms will break apart. But there are places on the river where these tiny pieces of ice can accumulate, forming mushy ice. There are many studies of mushy ice and ice formation on the north pole if you want to check. If the temperatures stay under the fusion point of the river water, the mushy ice will eventually become layers of ice. Once the layers form in the low-velocity stream places they will start growing and expanding. Because ice is solid, water will doge the ice, then, if the layer is thick enough water won't break it and will eventually manage to cover the whole river. Only the surface will freeze. Beneath the ice water still flows normally. You'll be surprised at how thin can the ice layer be without breaking.
To simulate this, I'll start with the Navier-stokes equations for the water motion. Then, you have to add some advection-diffusion heat transport equations. To couple them to the fluid just make the density of water a function of temperature. You could try to make viscosity also a function of temperature but that can result in some numerical instabilities. For the phase change, the enthalpy method is quite simple to implement. You know that most of the water on the river will be almost at the fusion point but won't lose enough energy to become solid. Just the water on top will. Then, you have to take into account that you have 3 mediums. Ice, water, and air. All those have different densities and heat capacitances that are important for the simulation. You'll notice that ice will become some sort of heat insulator preventing the underneath water to freeze after a certain thickness. Then the difficult part. You want the ice to move and stay on top.
So, you can have 2 approaches.
The first one is to use some boundary tracking method to separate ice from water. Then, calculate what the drag force would be and the buoyancy so the ice floats and apply that to the boundary to move it around. As you probably noticed, the problem with this approach is that the boundary has to move, and that can be tricky.
The second way is using some Brinkmann penalization for the fluid but you'll need an extra auxiliary field for the ice, and somehow, you have to calculate the forces from the fluid velocity field and pass them over to the auxiliary field. At last, you have to implement a fracture model for the ice. This model won't be that hard because ice is very well studied and there are plenty already proposed.
I hope you got a general idea of the physical phenomena. If you have any questions I'm happy to answer them.
