Large-scale river dynamics are often simulated using two-dimensional approaches based on deep-integrated Navier-Stokes equations, but capturing turbulent phenomena or flow behavior in slope-dominated river sections requires solving the full Navier-Stokes equations in three dimensions. In this study, we present a two-fluid fractional energy- and mass-conservative approach based on a standard Galerkin Finite Element (FE) formulation. A level-set method is used to track the evolution of the interface between water and air . The weak pressure discontinuity at the interface between the two fluids is handled through an enriched shape function space. To maintain the same degrees of freedom in the global system, a local condensation problem is solved within each element. Different boundary conditions are applied depending on the prevailing flow regime in the study section, ensuring appropriate treatment for both supercritical and subcritical regimes. This approach has been implemented within the Kratos Multiphysics framework, leveraging automatic terrain reconstruction procedures that facilitate the generation of realistic 3D models from Digital Elevation Models (DEM).