Abstract A stochastic parameterization of ocean mesoscale eddies is constructed in order to account for the fluctuations in subgrid transport and to represent upscale turbulent cascades. Eddy-resolving simulations to derive the parameterization are performed in a quasi-geostrophic (QG) model in a double-gyre configuration. The coarse-graining of the high-resolution model is giving rise to an eddy source term which represents the turbulent Reynolds stresses. The eddy source term, its mean and fluctuations are analyzed as function of the resolved scales and external parameters. A functional form of the resolved scales, based on a representation of turbulence as a non-Newtonian viscoelastic medium and including the rate of strain, is used to describe the eddy source term mean, variance and decorrelation timescale. Probability density functions (PDFs) of the eddy source term conditional on the resolved scales are then calculated, capturing the fluctuations associated with mesoscale eddies and their impact on the mean flow. Scalings for the mean, standard deviation, skewness, and kurtosis of the conditional PDFs are provided as function of the grid size, forcing, and stratification of the coarse-resolution model. In light of these scalings, no preliminary high-resolution (QG) model runs are necessary to diagnose the subgrid forcing and the implementation of a stochastic closure based on the conditional PDFs requires in principle very little tuning.