In this study we investigate a covariate-based stochastic approach to parameterize unresolved turbulent processes within a standard model of the idealised, wind-driven ocean circulation. We focus on vertical instead of horizontal coarse-graining, such that we avoid the subtle difficulties of horizontal coarsegraining. The corresponding eddy forcing is uniquely defined and has a clear physical interpretation related to baroclinic instability.We propose to emulate the baroclinic eddy forcing by sampling from the conditional probability distribution functions of the eddy forcing obtained from the baroclinic reference model data. These conditional probability distribution functions are approximated here by sampling uniformly from discrete reference values. We analyze in detail the different performances of the stochastic parameterization dependent on whether the eddy forcing is conditioned on a suitable flow-dependent covariate or on a timelagged covariate or on both. The results demonstrate that our non-Gaussian, non-linear methodology is able to accurately reproduce the first four statistical moments and spatial/temporal correlations of the stream function, energetics, and enstrophy of the reference baroclinic model.