I present synthetics of seismic wave propagation near free surface topography. The velocity-stress formulations of both the full elastic and viscoelastic wave equations are used, and I have derived exact boundary conditions for any arbitrary, smooth topography in terms of the particle velocities. Program codes are developed for 2 and 3 dimensions (2-D and 3-D) using finite-difference (F-D) methods for both spatial and temporal numerical discretizations. An 8th order F-D method is used inside the physical model space, and the spatial F-D order decreases gradually towards the free surface topography. The discretization of the medium equations along the side and bottom boundaries, the free surface topography boundary conditions, and the forward time stepping, are all by 2nd order F-D methods. The leap-frog technique is used for time stepping everywhere except for the memory variable equations in the viscoelastic cases, where an explicit version of the unconditionally stable Crank-Nicholson method is used. I show synthetics applying the schemes to isotropic 2-D and 3-D media covered by topographies that are either described by analytic expressions or by real elevation data. These data are taken from an area in South-Eastern Norway that contains the NORESS seismic receiver array. Domains up to 60 x 60 kilometers are used in 3-D simulations, and the applied sources are plane waves generated by a plane of Ricker type point sources. These sources represent earthquakes or teleseismic explosions. For 2-D simulations I have used both plane waves and point sources, since the larger models permissible in 2-D allow for point sources to represent earthquakes or teleseismic explosions quite well. For 2-D simulations I have also included examples using layered media with randomization by a 2-D von Karman function with and without apparent anisotropy. Synthetic snapshots and seismograms show Rayleigh (Rg)-waves emanating from areas of prominent topography as well as strong surface wave directivity from some topographic features. Full viscoelastic modeling with relatively low Q-values, describing near-surface sedimentary layers, exhibit intrinsic attenuation and physical dispersion of the wavefield. Results coincide with numerous observations. 3-D simulations are performed using domain decomposition parallelization implemented by Message Passing Interface (MPI).