We consider the efficient solution of partial differential equationsfor strongly elliptic operators with constant coefficients and stochastic Dirichlet data by the boundary integral equation method. The computation of the solution's two-point correlation is well understood if the two-point correlation of the Dirichlet data is known and sufficiently smooth.Unfortunately, the problem becomes much more involved in case of rough data. We will show that the concept of the H-matrix arithmetic provides a powerful tool to cope with this problem. By employing a parametric surface representation, we end up with an H-matrix arithmetic based on balanced cluster trees. This considerably simplifies the implementation and improves the performance of the H-matrix arithmetic. Numerical experiments are provided to validate and quantify the presented methods and algorithms