Abstract. We have developed a new statistical approach (M3Fusion) for combining surface ozone observations from thousands of monitoring sites around the world with the output from multiple atmospheric chemistry models to produce a global surface ozone distribution with greater accuracy than can be provided by any individual model. The ozone observations from 4766 monitoring sites were provided by the Tropospheric Ozone Assessment Report (TOAR) surface ozone database, which contains the world's largest collection of surface ozone metrics. Output from six models was provided by the participants of the Chemistry-Climate Model Initiative (CCMI) and NASA's Global Modeling and Assimilation Office (GMAO). We analyze the 6-month maximum of the maximum daily 8 h average ozone value (DMA8) for relevance to ozone health impacts. We interpolate the irregularly spaced observations onto a fine-resolution grid by using integrated nested Laplace approximations and compare the ozone field to each model in each world region. This method allows us to produce a global surface ozone field based on TOAR observations, which we then use to select the combination of global models with the greatest skill in each of eight world regions; models with greater skill in a particular region are given higher weight. This blended model product is bias corrected within 2∘ of observation locations to produce the final fused surface ozone product. We show that our fused product has an improved mean squared error compared to the simple multi-model ensemble mean, which is biased high in most regions of the world.