A method for simulating delay-Doppler maps (DDMs) of global navigation satellite system signals reflected from land surfaces with heterogeneous terrain is developed from first principles. The method follows previous work for ocean DDMs in the geometric optics limit of the Kirchhoff approximation. Unlike the ocean method, however, where surface heights are assumed to be random with homogeneous statistics, this method decomposes the surface heights into a deterministic part obtained from a digital elevation map (DEM) and a random part representing the residual between the surface and the DEM. The method accounts for the displacement of reflected power into bins of lower delay due to raised surface terrain. The method also provides for the modulation of the normalized bistatic radar cross section by DEM-derived surface slopes over the glistening zone of the DDM. A technique to register Cyclone Global Navigation Satellite System (CYGNSS) DDM bins in delay-Doppler space for land applications is also proposed. The DEM-based method is applied to a CYGNSS track over the Soil Moisture Sensing Controller And oPtimal Estimator (SoilSCAPE) site at Tonzi Ranch, CA, USA. The DEM-based method has potential application for spaceborne monitoring of a variety of environmental parameters.