Integrating crosshole ground-penetrating radar (GPR) with seismic methods is an efficient way to reduce the uncertainty and ambiguity of data interpretation in shallow geophysical investigations. We have developed a new approach for joint full-waveform inversion (FWI) of crosshole seismic and GPR data in the frequency domain to improve the inversion results of both FWI methods. In a joint objective function, three geophysical parameters (P-wave velocity, permittivity, and conductivity) are effectively connected by three weighted cross-gradient terms that enforce the structural similarity between parameter models. Simulation of acoustic seismic and scalar electromagnetic problems is implemented using 2D finite-difference frequency-domain methods, and the inverse problems of seismic FWI and GPR FWI are solved using a matrix-free truncated Newton algorithm. The joint inversion procedure is performed in several hierarchical frequencies, and the three parameter models are sequentially inverted at each frequency. The joint FWI approach is illustrated using three numerical examples. The results indicate that the joint FWI approach can effectively enhance the structural similarity among the models, modify the structure of each model, and improve the accuracy of inversion results compared with those of individual FWI approaches. Moreover, joint inversion can reduce the trade-off between permittivity and conductivity in GPR FWI, leading to an improved conductivity model in which artifacts are significantly decreased.