Fracture matrix interactions strongly affect anomalous heat transfer in geological sites. This study investigates effects of the interactions between fractures and rock matrix by using the method of multiple interacting continua (MINC). The MINC generates different temperature histories for varied fracture spacings. Two analytical solutions of each porous model and fracture model are used to fit the numerical results for temperature histories due to cold-water injection. The porous model is good agreement with the result for small fracture spacing, while a solution of the fracture model fits the result for large fracture spacing. The MINC yields intermediate behaviors in between a porous medium and a single fracture. A fractional heat transfer equation (fHTE) has been developed to describe anomalous thermal diffusion in a fractured reservoir. The fHTE accounts for heat flux from fracture into matrix by using a temporal fractional derivative. The fHTE can capture numerical results for temperature histories with different fracture spacings. The fracture spacing has correlations to the fHTE best-fit parameters (i.e., the orders of fractional derivatives and the retardation parameters). The fHTE with varying time fractional derivatives can cover descriptions of subdiffusion, Fickian diffusion, and superdiffusion. (C) 2016 Elsevier Ltd. All rights reserved.