A river engineer is challenged with the task of setting up an appropriate model for a certain application. The model needs to provide suitable answers to the questions asked (i.e. be effective) and needs to do this within the available time (i.e. be efficient). To set up such a model with sufficient accuracy and certainty, a modeller needs to fully understand all processes that determine the flow patterns and the flow resistance. These encapsulate both the physical processes, such as bottom friction and turbulent mixing, as well as the unwanted, ’numerical processes’, due to discretization errors and grid effects. Unfortunately, these errors can be considerably large and can greatly influence model results. To quantify the effects of numerical inaccuracies on the flow patterns and resistance (or backwater) in a river, several building blocks of the governing flow equations were analyzed. In particular for moderate resolutions, where a part of the geometrical variation in a river is captured on the grid, the influence of the momentum advection scheme and the turbulence model on the model results increases. For these modelling aspects, certain common methods were therefore analyzed concerning their accuracy, efficiency and convergence properties, for grid resolutions applied in practical engineering work. The common consensus is that the backwater in river models is dominated by bottom friction and that momentum advection only has a local effect on the water levels and flow patterns. However, in this work, it is demonstrated that the artificial backwater contribution from the momentum advection approximation can be of the same order of magnitude as the bottom friction contribution, depending on the advection scheme. First this is shown using a one-dimensional (1D) analysis and then it is verified using 1D and two-dimensional (2D) numerical experiments with a wavy bed, with emerged and submerged groynes and finally for an actual river. For each test, the backwater contribution of three basic first-order and two second-order accurate advection schemes are computed and compared. The size of this contribution is found to be largely determined by the conservation/constancy properties of the scheme and to a lesser extent by the order of the scheme. Of course, the bottom friction forms the most important contribution to the total backwater. In 2D models, the bottom friction computation is considered to be straightforward, in particular when applying the newly-developed subgrid method by Casulli and Stelling [51] and Stelling [219]. However, for three-dimensional (3D) hydrodynamic models, the computation is more complicated due to the vertical structure of the flow. Most 3D river models apply the popular σ-layering, where the grid nicely follows the bed and free-surface. At present, z-layer models are seldomly applied in river computations, because they suffer from the problem of inaccurate and discontinuous bottom shear stress representation, commonly assumed to arise due to the staircase bottom representation. At higher grid resolution, where more features of the topography are represented on the grid, a terrain-following coordinate system such as the σ-layering can result in a strong distortion of the grid. This is avoided using a z-layer discretization. Additionally, the latter discretization could be very efficient for river applications, due to the fact that excessive vertical resolution is avoided in shallow areas, such as floodplains. For this purpose, the discretized equations for the z-layer model are analyzed and the cause of the inaccuracies is clearly shown to come from the emergence of thin near-bed layers. Based on this analysis, a new method is presented that significantly reduces the errors and the grid dependency of the results. The method consists of a near-bed layer-remapping and a modified near-bed discretization of the k-ε turbulence model. The applicability of the approach is demonstrated for uniform channel flow, using a schematized 2D vertical (2DV) model and for the flow over a bottom sill using the Delft3D modelling system (Deltares [69]). Finally a new modelling strategy is presented for improving the efficiency of computationally intensive flow problems in environmental free-surface flows. The approach combines the recently developed semi-implicit subgrid method by Casulli and Stelling (Casulli [46], Casulli and Stelling [51], and Stelling [219]) with a hierarchical-grid solution strategy. The method allows the incorporation of high-resolution data on subgrid scale to obtain a more accurate and efficient hydrodynamic model. The subgrid method improves the efficiency of the hierarchical grid method by providing better solutions on coarse grids. The method is applicable to both steady and unsteady flows, but it is particularly useful in river computations with steady boundary conditions. There, the combined hierarchical grid-subgrid method reduces the computational effort to obtain a steady state with factors up to 43. For unsteady models, the method can be used for efficiently generating accurate initial conditions and further dynamic computations on high-resolution grids. Additionally, the method provides automatic insight in grid convergence. The efficiency and applicability of the method is demonstrated using a schematic test for the vortex shedding around a circular cylinder and a real-world case study on the Elbe River in Germany. @en