Earth System Models (ESMs) are the primary tool for understanding the impacts of global change and several ESMs are updated on a regular basis to provide more reliable scenarios of the future. However, the confrontation of ESMs outputs to observations reveals biases that are important to correct, especially for impact applications where the absolute scale of the environmental variable is as relevant as its trends. In addition, regional impact studies need fine scale projections to devise strategic planning and management measures. Statistical downscaling provides a fast way to produce regional ocean forcing from ESMs and can additionally produce bias-corrected outputs, which are necessary for impact applications driven by or fitted to observed data, like many ecological models. Statistical downscaling can make use of different parametric distributions depending on the variables used, and generalized regression can provide a flexible approach for this purpose. We propose a multi-model approach based on non-parametric generalized regression and a suite of indicators to select a robust statistical downscaling model that can be used for projection of future scenarios. The empirical cumulative distribution of the variables to downscale is modeled, ensuring that not only the mean but also the variance and quantiles (including the minima and maxima) are properly represented, improving the prediction of extreme events and taking into account spatial autocorrelation. The approach presented here is applied to two contrasted regional case studies, the Bay of Biscay-Celtic Sea ecosystem and the Northern Peru Current ecosystem, using the Sea Surface Temperature from the IPSL-CM5A-LR ESM. The results showed that a multimodel selection approach is appropriate as individual model performance is case specific.