Accurately predicting urban environmental airflows using the Reynolds-averaged Navier‒Stokes (RANS) approach presents significant challenges due to unsolved closure issues. Previous attempts to address these limitations through coefficient adjustments were constrained by isotropic assumptions, resulting in limited generalizability and reliability, especially for complex building configurations. In this study, we proposed a recalibration approach within a multi-objective optimization framework. Nondominant sorting genetic algorithm-II (NSGA-II) and computational fluid dynamics (CFD) numerical computations were utilized to optimize the closure coefficients of a nonlinear eddy viscosity (NLEV) model, with a focus on improving the accuracy of the mean velocity, turbulent kinetic energy, and airflow characteristic predictions. To ensure the robustness of the model, appropriate objective functions were carefully defined. The proposed approach was evaluated using three-dimensional (3D) building cases related to urban prototypes, and the results demonstrated the necessity of multi-objective optimization. Our findings indicated that the trade-off solution on the Pareto front consistently outperformed the baseline and single-objective optimal versions across different types of urban airflow, demonstrating superior generalizability. The simulation results of this solution exhibited closer agreement with the wind tunnel experimental data and provided more accurate distributions of the time-averaged values and airflow characteristics. The normalized root mean square errors are reduced by about 33 %, 52 %, 35 % and 9 % in velocity, and 31 %, 25 %, 24 % and 11 % in turbulent kinetic energy, respectively. These results underscore the effectiveness of the multi-objective optimization framework in addressing closure issues and improving the prediction accuracy of complex urban environmental airflows.