In a heterogeneous medium (usually called a rough medium) with fractured formations, the propagation of an electromagnetic (EM) field is a type of subdiffusion. Current mainstream geophysical EM data processing methods cannot be applied to data acquired on heterogeneous Earth, as they are not governed by the classic diffusion theory. To evaluate the influence of roughness on the transient EM (TEM) signal for a complex model and contribute to data inversion, we proposed a novel three-dimensional (3-D) forward modeling scheme for TEM in rough media. First, we derived the governing equation with a fractional-order time derivative for the subdiffusion of EM waves in rough media. Then, we proposed a novel time discretization using an unequal step length for the Caputo operator, which significantly reduces the total number of time steps. Finally, an implicit time-domain finite-element method using unstructured tetrahedron discretization was adopted to solve the 3-D forward problem. Furthermore, an efficient time segmentation strategy combined with parallel RHS construction was proposed to accelerate modeling. The numerical results prove that the proposed method is accurate and efficient, and will be a powerful numerical method for analyzing TEM wave propagation and processing TEM data in areas with multiscale fractures or porosity.