We have developed a numerical algorithm for computing the magnetic field distribution and the nuclear magnetic resonance (NMR) sensitivity function for arbitrary topography overlying a known 3-D conductivity structure. The magnetic vector potential is split into primary and secondary terms. The primary term is obtained using a thin-wire line integral equation that accounts for arbitrary loop shape and position. It allows the singularity of the source field to be effectively removed. The secondary potential is obtained by solving the second-order partial differential equations on an unstructured tetrahedral mesh using the finite element technique. We validate the results of applying our algorithm against an explicit infinite integral solution for circular loops on a layered earth and against the results of applying a commercial simulation tool. The spatially oscillating NMR sensitivity functions to hydrogen protons (i.e., unbounded water molecules) in the sub-surface are computed on a refined unstructured grid. We apply the numerical algorithm to a number of synthetic examples in surface NMR tomography of hydrological relevance.