A major weakness in one-dimensional (1D) stellar structure and evolution modelling is the simplified treatment of convection, which leads to erroneous near-surface stratification and considerable uncertainties in predicted effective temperatures and luminosities of low-mass stars. In a series of preceding works, a novel method for coupling 1D stellar structural models with a grid of 3D surface convection simulations during stellar evolution was developed, at solar metallicity. This 1D-3D coupling method slightly shifts evolutionary tracks relative to standard calculations, meanwhile providing oscillation frequencies that agree more closely with asteroseismic observations. Here, we extend this method to model metal-poor and metal-rich F, G, Ktype stars, by implementing interpolations on-the-fly across metallicity ( -3 < [Fe / H] < 0 . 5) for mean 3D models during stellar evolution. We demonstrate quantitatively that the fundamental stellar parameters modelled within our framework are insensitive to the mixing-length parameter. A 20 per cent change in the mixing-length parameter results in evolutionary tracks with a temperature shift of less than 30 K, compared to a difference of over 200 K in standard evolution calculations. Our extension is validated against eclipsing binary systems with extremely precise observational constraints as well as stars in binaries with asteroseismic data. Using a fixed mixing-length parameter that merely governs convective heat transport in the near-adiabatic layers, the 1D-3D coupling method successfully reproduces most observational constraints for all target stars. Coupling 1D stellar evolution models with 3D simulations greatly reduces uncertainties associated with the choice of atmosphere boundary conditions and mixing-length parameters, hence offering a powerful tool for characterizing stars with seismic measurements and determining ages for globular clusters.