Physical parameterizations in global atmospheric and ocean models typically include free parameters that are not theoretically or empirically constrained. New methods are required to determine the optimal parameter combinations for such models in an objective, exhaustive, yet computationally feasible manner. Here, we propose to apply computationally inexpensiveradial basis function (RBF) surrogate models to minimize a “cost”, or error, function of an atmospheric model or a physical parameterization. The RBF is iteratively updated as more input‐output pairs are obtained during the optimization. The approach is used to optimize the Eddy‐Diffusivity/Mass‐Flux (EDMF) boundary‐layer parameterization of the convective boundary‐layer in a single‐column (SCM) framework. The optimization based on surrogate models is able to identify parameter combinations that reduce the error of the untuned default setting by 41%. The probability to detect a low‐error solution increases rapidly especially over the first tens of SCM evaluations. In comparison, a quadratic polynomial model only yields an error reduction of 17% since a) high‐order parameter interactions are not accounted for and b) the construction of the polynomial is not based on an iterative sampling approach. The RBF surrogate models achieve this 17% error reduction for 40% of the polynomial model's cost. Interestingly, one of the emerging optimal settings describes a pure mass flux parameterization without eddy‐diffusivity component. A second optimal solution is characterized by a plume fractional area of 20‐30% and an eddy‐mixing time scale of ∼700 s.