Correlation operators are fundamental for representing background and observation errors in variational assimilation and for localizing low-rank sample covariance matrices in ensemble-variational assimilation. Furthermore, a factored (“square-root”) form of a correlation operator can be used to generate Gaussian random fields and thus is an important tool for perturbing input parameters in ensemble methods. This presentation describes recent advances in the development of diffusion-based correlation operators and their application to global ocean data assimilation. The starting point is an implicitly formulated diffusion operator, which has been shown in previous works to support symmetric and positive smoothing kernels that are closely related to those from the Matérn correlation family. We propose an iterative algorithm based on the Chebyshev iteration to solve the sequence of linear systems defining the implicit diffusion process. We discuss its attractive features compared to other methods, in terms of efficiency, accuracy, memory cost, ease-of-implementation and parallelization properties on high-performance computers. This work has been supported by LEFE-MANU.