  We propose a local modification of the standard subdiffusion model by
introducing the initial Fickian diffusion, which results in a multiscale
diffusion model. The developed model resolves the incompatibility between the
nonlocal operators in subdiffusion and the local initial conditions and thus
eliminates the initial singularity of the solutions of the subdiffusion, while
retaining its heavy tail behavior away from the initial time. The
well-posedness of the model and high-order regularity estimates of its
solutions are analyzed by resolvent estimates, based on which the numerical
discretization and analysis are performed. Numerical experiments are carried
out to substantiate the theoretical findings.
