We report on the first multi-year kilometre-scale global coupled simulations using ECMWF's Integrated Forecasting System (IFS) coupled to both the NEMO and FESOM ocean-sea ice models, as part of the H2020 Next Generation Earth Modelling Systems (nextGEMS) project. We focus mainly on an unprecedented IFS-FESOM coupled setup, with an atmospheric resolution of 4.4 km and a spatially varying ocean resolution that reaches locally below 5 km grid spacing. A shorter coupled IFS-FESOM simulation with an atmospheric resolution of 2.8 km has also been performed. A number of shortcomings in the original numerical weather prediction (NWP)-focused model configurations were identified and mitigated over several cycles collaboratively by the modelling centres, academia, and the wider nextGEMS community. The main improvements are (i) better conservation properties of the coupled model system in terms of water and energy budgets, which also benefit ECMWF's operational 9 km IFS-NEMO model; (ii) a realistic top-of-the-atmosphere (TOA) radiation balance throughout the year; (iii) improved intense precipitation characteristics; and (iv) eddy-resolving features in large parts of the mid- and high-latitude oceans (finer than 5 km grid spacing) to resolve mesoscale eddies and sea ice leads. New developments at ECMWF for a better representation of snow and land use, including a dedicated scheme for urban areas, were also tested on multi-year timescales. We provide first examples of significant advances in the realism and thus opportunities of these kilometre-scale simulations, such as a clear imprint of resolved Arctic sea ice leads on atmospheric temperature, impacts of kilometre-scale urban areas on the diurnal temperature cycle in cities, and better propagation and symmetry characteristics of the Madden-Julian Oscillation.