Estimating the joint probability density of seasonal geophysical variables like Aerosol Optical Depth (AOD) and Day of Year (DOY) presents three unresolved challenges. The first of these challenges is the periodic nature of the temporal axis. The second is the physically distinct scales of the two variables. The third is the marginal inconsistency introduced by smoothing operations. Existing methods for estimating probability density do not address each of these challenges simultaneously. Here we introduce CADS (the Circular-Adaptive Density Smoother), a computationally efficient algorithm for estimating the joint probability density of two seasonal geophysical variables that simultaneously addresses each of these three challenges. CADS is evaluated on two synthetic datasets and one real observational dataset of AERONET measurements from NASA Ames ($n=74,653$) using five-fold cross-validation. CADS achieves the highest mean log-likelihood relative to other methods for estimating probability density, and it is approximately 3000 times faster than kernel density estimation. An ablation study confirms that each component of CADS contributes independently to its high performance. Finally, a novel metric for evaluating the geometric correctness of the treatment of the circular boundary of the DOY variable is introduced.