In this paper, we propose some new box relaxation numerical schemes on staggered grids to solve the stationary system of partial differential equations arising from the dual minimization problem associated with the total variation operator. We present in detail the numerical schemes for the scalar case and its generalization to multichannel (vectorial) images. Then, we discuss their implementation in digital image denoising. The results outperform the resolution of the dual equation based on the gradient descent approach and pave the way for more advanced numerical strategies.