In this study, we examine the modeling of the naturally fractured reservoirs based on the fractal, Metzler, and Raghavan anomalous diffusion models which are different from Darcy flux law. The governing equations of these models are presented for a vertical well producing at either a constant rate (CR) or constant bottomhole pressure (CBHP) production in an infinite-acting reservoir. New analytical and early- and late-time asymptotic approximate solutions were derived for the finite-wellbore problem for the Raghavan anomalous model. The approximate solutions identify the fractal and anomalous diffusion parameters that influence behaviors of the early- and late-time pressure and rate transient data. The method of Laplace transformation is used to find the solutions in dimensionless forms which are then inverted to the real-time domain using a numerical inversion algorithm. In addition, we compare the solutions for the fractal, Metzler anomalous, and Raghavan anomalous models and delineate the similarities and differences among these solutions. We show the effect of fractal parameters on diffusion models, the relationship between the fractal model, and Raghavan and Metzler anomalous models, in terms of dimensionless pressure and dimensionless pressure derivative responses in the reservoir. We show that the “dimensionless” pressure defined by Raghavan is not a dimensionless quantity unless the Euclidean dimension d=2 (or for a fractal structure with the mass fractal dimension df=2). To the best of our knowledge, in the literature, there exists no study providing a detailed review and comparing the applicability of each diffusion model on real field well-test data as well as providing new analytical solutions and early- and late-time asymptotic approximate solutions for the anomalous diffusion for both the CR and CBHP cases for the finite-wellbore problem. Hence, this study fills this gap. In addition, we apply the three diffusion models to an interference well test data from the Kizildere geothermal reservoir in Turkey by using the ensemble smoother with multiple data assimilation (ES-MDA) method. The results show that the fractal model in terms of the goodness of fit (root-mean-square error – RMSE) best matches the field test data.