Spatiotemporal (4D) neonatal cortical surface atlases with densely sampled ages are important tools for understanding the dynamic early brain development. Conventionally, after non-linear co-registration, surface atlases are constructed by simple Euclidean average of cortical attributes across different subjects, which leads to blurred folding patterns and therefore hampers the reliability and accuracy when used to register new subjects. To better preserve the sharpness and clarity of cortical folding patterns on surface atlases, we propose to compute the Wasserstein barycenter, which represents a geometrically faithful population mean under the Wasserstein distance metric, for the construction of 4D neonatal surface atlases. Comparing to the direct vertex-wise Euclidean average, the Wasserstein distance takes into account the alignment of spatial distribution of cortical attributes, thus is robust to potential registration errors during atlas building. Using this method, we constructed 4D neonatal cortical surface atlases at each week, from 39 to 44 postmenstrual weeks, based on a large-scale dataset with 764 subjects. Our 4D atlases show sharper and more geometrically faithful cortical folding patterns compared to the state-of-the-art methods, thus leading to boosted accuracy when used to align new subjects and facilitating early brain development studies.