State-of-the-art climate models rely on bulk formulae arising from the Monin-Obukhov semi-empirical theory to estimate turbulent air-sea fluxes. The mathematical structure of those formulae implies several difficulties when trying to study the numerical properties of coupling algorithms used for practical applications. This paper introduces a methodology for building physically realistic approximations of existing bulk formulae that would also satisfy suitable mathematical properties (explicit character, regularity, differentiability). This is achieved by applying the Sobol' method to compute sensitivity indices in order to reduce the number of inputs and derive a simple metamodel for the parameterization of turbulent air-sea fluxes. Numerical results show excellent agreement between our approximations and the standard bulk formulae. In particular, single column simulations using the TOGA-Coare experiment within the LMDZ atmospheric model show negligible changes in numerical results.