A set of model equations for water–wave propagation is derived by piecewise integration of the primitive equations of motion through two arbitrary layers. Within each layer, an independent velocity profile is derived. With two separate velocity profiles, matched at the interface of the two layers, the resulting set of equations has three free parameters, allowing for an optimization with known analytical properties of water waves. The optimized model equations show good linear wave characteristics up to kh ≈ 6, while the second–order nonlinear behaviour is captured to kh ≈ 6 as well. A numerical algorithm for solving the model equations is developed and tested against one– and two–horizontal–dimension cases. Agreement with laboratory data is excellent.