We study a random matching problem on closed compact $2$-dimensional Riemannian manifolds (with respect to the squared Riemannian distance), with samples of random points whose common law is absolutely continuous with respect to the volume measure with strictly positive and bounded density. We show that given two sequences of numbers $n$ and $m=m(n)$ of points, asymptotically equivalent as $n$ goes to infinity, the optimal transport plan between the two empirical measures $μ^n$ and $ν^{m}$ is quantitatively well-approximated by $\big(\mathrm{Id},\exp(\nabla h^{n})\big)_\#μ^n$ where $h^{n}$ solves a linear elliptic PDE obtained by a regularized first-order linearization of the Monge-Ampère equation. This is obtained in the case of samples of correlated random points for which a stretched exponential decay of the $α$-mixing coefficient holds and for a class of discrete-time Markov chains having a unique absolutely continuous invariant measure with respect to the volume measure.