Numerical simulation technique of surface nuclear magnetic resonance (NMR) responses from stratified groundwater formation with finite conductivity is concerned with two problems. One is the integration of Hankel transform of a rapidly oscillating kernel function, and another is the volume integration of surface NMR signal induced in a receiver loop from the protons in buried aquifer. In this paper, a direct numerical integration method is applied to the oscillating kernel function integration to obtain radial and vertical components of magnetic fields from a large loop source, and the spatial characteristics of the fields are described. Due to the distribution of magnetic fields, an ununiform space interval subdivision method has been developed to estimate surface NMR response in the receiver loop in order to save computing time. The surface NMR responses from different models have been evaluated, and influences of model conductivity, aquifer depth and its thickness, and loop size have been discussed. Numerical results show that the combination of direct numerical integration and ununiform space interval subdivision is feasible to get surface NMR signal. And the conductivity of formation is the key to generate the phase of surface NMR signal, but other factors, such as groundwater depth and its thickness, loop size, can make the amplitude and phase of surface NMR response greatly different.