Accurate characterization of aerosol optical depth (AOD) is crucial for atmospheric science, serving as a cornerstone metric in aerosol-climate process research and radiative forcing quantification. However, satellite-based AOD observations contain data gaps due to cloud contamination and high surface reflectance. Appropriate statistical modeling of the spatio-temporal AOD process is therefore essential to filling those data gaps, so as to acquire a reliable and gap-free dataset. Owing to the high-dimensional nature of the underlying interpolation problem, traditional geostatistical approaches are computationally infeasible. In this study, a spatial mixed effect model is used to describe the AOD process, which separates the process into large-, small-, and fine-scale spatial variations, modeled through covariates and spatial basis functions. Since the dimension of the basis functions is fixed and is much smaller than the original data dimension, the method is known as fixed rank kriging (FRK), which makes the spatial prediction tractable. The empirical part of this study applies FRK on monthly and daily AOD observations from MODIS. On the monthly dataset, using the AERONET observations as the ground truth, FRK AOD is found to outperform several benchmarks, including inverse distance weighting interpolation, DeepKriging, and the MERRA-2 reanalysis. Additionally, on the daily dataset, FRK is shown to be able to effectively handle large-scale data, whereas other alternatives often fail due to computational infeasibility.