We present a coupled thermal-fluid model for Bénard-Marangoni convection in a three-dimensional fluid layer. The governing equations are derived in detail for two reasons: first, we do not assume a flat free surface as commonly done; and second, we prepare for the use of flexible discretizations. The governing equations are discretized using spectral elements in space and an operator splitting approach in time. Since we are here primarily interested in steady state solutions, the focus is on the spatial discretization. The overall computational approach is very attractive to use for several reasons: (i) the solution can be expected to have a high degree of regularity, and rapid convergence can be expected; (ii) the spectral element decomposition automatically gives a convenient parameterization of the free surface that allows powerful results from differential geometry to easily be exploited; (iii) free surface deformation can readily be included; (iv) both normal and tangential stresses are conveniently accounted for through a single surface integral; (v) no differentiation of the surface tension is necessary in order to include thermocapillary effects (due to integration-by-parts twice); (vi) the geometry representation of the free surface need only be C⁰ across element boundaries even though curvature effects are included. Three-dimensional simulation results are presented, including the free surface deflection due to buoyancy and thermocapillary effects.