  We present here the cosmo-SLICS, a new suite of simulations specially
designed for the analysis of current and upcoming weak lensing data beyond the
standard two-point cosmic shear. We sample the $[\Omega_{\rm m}, \sigma_8, h,
w_0]$ parameter space at 25 points organised in a Latin hyper-cube, spanning a
range that contains most of the $2\sigma$ posterior distribution from ongoing
lensing surveys. At each of these nodes we evolve a pair of $N$-body
simulations in which the sampling variance is highly suppressed, and ray-trace
the volumes 800 times to further increase the effective sky coverage. We
extract a lensing covariance matrix from these pseudo-independent light-cones
and show that it closely matches a brute-force construction based on an
ensemble of 800 truly independent $N$-body runs. More precisely, a Fisher
analysis reveals that both methods yield marginalized two-dimensional
constraints that vary by less than 6% in area, a result that holds under
different survey specifications and that matches to within 15% the area
obtained from an analytical covariance calculation. Extending this comparison
with our 25 $w$CDM models, we probe the cosmology dependence of the lensing
covariance directly from numerical simulations, reproducing remarkably well the
Fisher results from the analytical models at most cosmologies. We demonstrate
that varying the cosmology at which the covariance matrix is evaluated in the
first place might have an order of magnitude greater impact on the parameter
constraints than varying the choice of covariance estimation technique. We
present a test case in which we generate fast predictions for both the lensing
signal and its associated variance with a flexible Gaussian process regression
emulator, achieving an accuracy of a few percent on the former and 10% on the
latter.
