Soil temperature (ST) and its spatiotemporal variations are critical factors influencing soil energy exchange and physical processes. This study proposes a novel method for estimating ST at depth using moderate resolution imaging spectroradiometer (MODIS) time series imagery. Land surface temperature (LST) and soil heat flux (SHF) time series images are gap-filled by four daily acquisitions and modeled hourly through daily cycle modeling. Soil thermal admittance is calculated as the ratio of SHF amplitude to LST amplitude at different frequencies. Utilizing a multilayer perceptron model and soil thermal admittance images, soil thermal diffusivity and damping depth are derived. Analysis in arid regions reveals that damping depth increases with soil depth. ST is estimated hourly using a sinusoidal equation and MODIS imagery. Comparison of thermal behavior shows the highest accuracy at 5 cm depth and the lowest at 20 cm. Except for 5 cm depth, where accuracy is lowest at 12:30, the lowest accuracy at other depths occurs at 6:30. Monthly average ST variations from the surface to depth are analyzed annually. A new model is then applied to improve ST estimation through daily temperature cycle modeling. The proposed model uniquely integrates MODIS-derived LST and SHF with dynamic daily cycle parameters, such as sunrise, sunset, and day length, to simulate hourly ST. This integration distinguishes it from previous models by enabling more precise temporal reconstruction of subsurface thermal behavior. This model incorporates sunrise, sunset, day length, and local noon time for each day, enabling hourly ST estimation at 5 cm (Z1), 10 cm (Z2), and 20 cm (Z3) depths by considering damping depth, phase, and diurnal variations. At Z1, the estimated ST has a root mean square error (RMSE) < 2.8 degrees C and mean absolute error (MAE) < 2.4 degrees C. At Z2, RMSE decreases by similar to 1 degrees C at 6:30 compared to previous methods, with RMSE <2.6 degrees C and MAE < 2.8 degrees C. At Z3, RMSE is <2.3 degrees C. ST is also estimated hourly for depths of 30, 50, and 100 cm. The proposed method demonstrates higher accuracy in estimating ST compared to previous approaches, offering a robust solution for spatiotemporal ST analysis.