The design of metasurfaces and metasurface-based devices (e.g. reflective intelligent surfaces - RIS) is inherently affected by unpredictable variations of the unit-cell response, due, for instance, to fabrication tolerances, physical discontinuities, presence of bumps, etc. Such effects may dramatically affect the accuracy of the conventional design models, such as the one based on the local surface impedance (LSI). In this work, we present a statistical model resulting in a probability density function that allows defining a random local surface impedance (rLSI), that takes properly into account potential random variations derived from the von Mises distribution. To compensate for these variations, we propose a transformation kernel that matches the expected rLSI value to the design target. By using an analytical formulation and a proper set of full-wave simulations, we demonstrate the effectiveness and robustness of the method against unpredictable variations caused by the mutual coupling among the unit-cells and fabrication tolerances/errors. The proposed method paves the way to foster advancements in realistic metasurface engineering.