For the simulation of therapeutic ultrasound applications, a method including frequency-dependent attenuation effects directly in the time domain is highly desirable. This paper describes an efficient numerical time-domain implementation of the power-law attenuation model presented by Szabo [Szabo, J. Acoust. Soc. Am. 96, 491-500 (1994)]. Simulations of therapeutic ultrasound applications are feasible in conjunction with a previously presented finite differences time-domain (FDTD) algorithm for nonlinear ultrasound propagation [Ginter et al., J. Acoust. Soc. Am. 111, 2049-2059 (2002)]. Szabo implemented the empirical frequency power-law attenuation using a causal convolutional operator directly in the time-domain equation. Though a variety of time-domain models has been published in recent years, no efficient numerical implementation has been presented so far for frequency power-law attenuation models. Solving a convolutional integral with standard time-domain techniques requires enormous computational effort and therefore often limits the application of such models to 1D problems. In contrast, the presented method is based on a recursive algorithm and requires only three time levels and a few auxiliary data to approximate the convolutional integral with high accuracy. The simulation results are validated by comparison with analytical solutions and measurements.