In this work, we have constructed a high dimensional neural network (NN) potential energy function for simulating palladium surface properties. The NN potential was trained with 3035 density functional theory (DFT) calculations, and was shown to be nearly as accurate as DFT in molecular simulations. Important properties including lattice constants, elastic properties and surface energies as well as transition state energies and adatom diffusion barriers were predicted by the NN and were found to be in excellent agreement with DFT results. The computational time to run the NN was compared to DFT calculation time, and we found this implementation of the NN is roughly four orders of magnitude faster than DFT. This approach is general and applicable to other systems and may have applications in modeling catalytic processes at surfaces.