Next-generation optical imaging surveys will revolutionize the observations of weak gravitational lensing by galaxy clusters and provide stringent constraints on growth of structure and cosmic acceleration. In these experiments, accurate modelling of covariance matrices of cluster weak lensing plays the key role in obtaining robust measurements of the mean mass of clusters and cosmological parameters. We use a combination of analytical calculations and high-resolution N-body simulations to derive accurate covariance matrices that span from the virial regime to linear scales of the cluster-matter cross-correlation. We validate this calculation using a public ray-tracing lensing simulation and provide a software package for calculating covariance matrices for a wide range of cluster and source sample choices. We discuss the relative importance of shape noise and density fluctuations, the impact of radial bin size, and the impact of off-diagonal elements. For a weak lensing source density $n_s = 10 \mathrm{arcmin}^{-2}$, shape noise typically dominates the variance on comoving scales $r_p \lesssim 5 h^{-1} \mathrm{Mpc}$. However, for $n_s = 60 \mathrm{arcmin}^{-2}$, potentially achievable with future weak lensing experiments, density fluctuations typically dominate the variance at $r_p \gtrsim 1 h^{-1} \mathrm{Mpc}$ and remain comparable to shape noise on smaller scales.