In the following the numerical solution of the kinematic boundary condition of the free surface is described. Eq.(5) is discretized using a finite difference method, with the node points at the vertex points of the free surface boundary of the grid. Since the bulk flow is discretized using a finite volume method, with the variables stored in the center of each grid cell, the velocities has to be extrapolated from the interior of the grid to the vertex of the grid point at the free surface. The discretized two dimensional version of Eq.(5) is written

(14)

where , the free surface is located at *j*=*j*_{max} , and is the vertical distance the surface at *x*=*x*_{i} moves at time step *n* .
Eq.(15) is solved using a three stage Runge Kutta routine.

(15)

where is the term within the brackets in Eq.(15), solved for The iteration scheme gives the new location of the free surface. The iteration of the bulk flow equations and the free surface equations is done in two independent loops. They are not coupled in a way that would prevent flux to go through the surface in every intermediate time step. However, as the solution converges the flux through the surface will go to zero, the velocity at the boundary will become parallel to the surface, and the surface will behave as a material surface. This approach follows the method presented in [6].

When the new location of the free surface is known, the grid is adjusted to fit the new boundary. This is done by vertical shifting the grid lines closest to the free surface. The number of grid lines to be moved is user defined, given in advance. The grid lines located deeper than a grid line, *j*=*j*_{fixed} remains unchanged throughout the simulation. The *x*– component of the vertex points remains unchanged in the deformed part of the grid,

(16)

whereas the

(17)

After the grid is adjusted to the new location of the free surface, the metrics of the deformed grid are computed. Thereafter, the bulk flow can be further advanced in time.