Optoacoustic (photoacoustic) mesoscopy offers unique capabilities in skin imaging and resolves skin features associated with detection, diagnosis and management of disease. A critical first step in the quantitative analysis of clinical optoacoustic images is to identify the skin surface in a rapid, reliable and automated manner. Nevertheless, most common edge-and surface-detection algorithms cannot reliably detect the skin surface on 3D raster-scan optoacoustic mesoscopy (RSOM) images, due to discontinuities and diffuse interfaces in the image. We present herein a novel dynamic programming approach that extracts the skin boundary as a 2D surface in one single step, as opposed to consecutive extraction of several independent 1D contours. A domain-specific energy function is introduced, taking into account the properties of volumetric optoacoustic mesoscopy images. The accuracy of the proposed method is validated on scans of the volar forearm of 19 volunteers with different skin complexions, for which the skin surface has been traced manually to provide a reference. Additionally, the robustness and the limitations of the method are demonstrated on data where the skin boundaries are low-contrast or ill-defined. The automatic skin surface detection method can improve the speed and accuracy in the analysis of quantitative features seen on RSOM images and accelerate the clinical translation of the technique. Our method can likely be extended to identify other types of surfaces in RSOM and other imaging modalities.