The paperpresents a comprehensive approach for full-wave modeling of grounding systems buried in a multilayer stratified ground. The dyadic spectral domain Green's function associated with the multilayer stratified ground is first derived. This Green's function is then evaluated by numerical solution of Sommerfeld's integrals along the deformed path in the complex plane. The proposed method is general as it can deal with arbitrarily shaped grounding electrodes buried in soils with any number of layers. The case of horizontal, vertical, and oblique electrodes buried in soils with different layers is analyzed. The accuracy of the proposed approach is first validated by comparing the results with those obtained byNEC-4 for a homogeneous soil, and with those obtained by finite difference time domain method for a multilayer soil.