Tropical cyclones (TCs), driven by heat exchange between the air and sea, pose a substantial risk to many communities around the world. Accurate characterization of the subsurface ocean thermal response to TC passage is crucial for accurate TC intensity forecasts and an understanding of the role that TCs play in the global climate system. However, that characterization is complicated by the high-noise ocean environment, correlations inherent in spatiotemporal data, relative scarcity of in situ observations, and the entanglement of the TC-induced signal with seasonal signals. We present a general methodological framework that addresses these difficulties, integrating existing techniques in seasonal mean field estimation, Gaussian process modeling, and nonparametric regression into an ANOVA decomposition model. Importantly, we improve upon past work by properly handling seasonality, providing rigorous uncertainty quantification, and treating time as a continuous variable, rather than producing estimates that are binned in time. This ANOVA model is estimated using in situ subsurface temperature profiles from the Argo fleet of autonomous floats through a multistep procedure, which (1) characterizes the upper-ocean seasonal shift during the TC season, (2) models the variability in the temperature observations, and (3) fits a thin-plate spline using the variability estimates to account for heteroskedasticity and correlation between the observations. This spline fit reveals the ocean thermal response to the TC passage. Through this framework, we obtain new scientific insights into the interaction between TCs and the ocean on a global scale, including a three-dimensional characterization of the near-surface and subsurface cooling along the TC storm track and the mixing-induced subsurface warming on the track’s right side.