[Abstract] Grounding systems are designed to preserve human safety and grant the integrity of equipments under fault conditions. To achieve these goals, the equivalent electrical resistance of the system must be low enough to ensure that fault currents dissipate (mainly) through the grounding electrode into the earth, while maximum potential gradients between close points on the earth surface must be kept under certain tolerances (step and touch voltages) [1,2]. In this paper we present a Boundary Element approach for the numerical computation of grounding systems. In this general framework, former intuitive widespread techniques (such as the Average Potential Method) can be identified as the result of specific choices for the test and trial functions, as well as suitable assumptions introduced in the BEM formulation to reduce computational cost. Linear and higher order elements can be used in order to increase accuracy avoiding excessive segmen- tation. On the other hand, computing time is kept under acceptable levels by means of analytical integration techniques and semi-iterative methods for solving linear equations systems. Finally, an application to a real problem is presented.