Single-well push-pull tracer tests are broadly employed to estimate effective parameters for solute and heat transport in aquifers. Tracer recovery curves obtained from these tests serve as inputs for solving an inverse problem to infer effective transport parameters such as porosity, thermal and solute longitudinal dispersivities, and retardation factors. However, the inherent non-uniqueness of the inverse calibration problem and associated uncertainties in field measurements create a bottleneck for multiparametric calibration. To address these challenges, we employed a computationally efficient optimization framework based on surrogate modeling via Gaussian process regression (GPR) to approximate the objective function based on six effective transport parameters to be calibrated simultaneously, which yields plausible parameter combinations. For training and model evaluation, we implemented a 1D finite-difference (FD) representation of the advection-dispersion equation for sorbing tracers featuring an adaptive explicit time stepping scheme adhering to numerical stability criteria while minimizing numerical diffusion, where an analytical radial flow field serves as input based on well hydraulic properties. The FD model includes the measured input time series of temperature and concentration as transient boundary conditions, as well as well-bore storage to accurately model push-pull test conditions. We applied this framework to push-pull tests conducted in a sandy aquifer in Horonobe (Hokkaido, Japan) using heat and three solute tracers: uranine, lithium, and iodide. The confidence intervals for field measurements were included by using repeated under identical conditions tests. The surrogate model facilitates parameter optimization by balancing the exploration of high-uncertainty regions with the exploitation of high-probability regions through a weighted probability function. The posterior parameter distribution reveals reduced uncertainty intervals for porosity and both solute and thermal dispersivities while indicating low sensitivity for the solute retardation factor. The results demonstrate the necessity of high-precision measurements for concentration and highlight the value of utilizing multiple tracers to enhance calibration accuracy under parametric and measurement uncertainty. The developed framework highlights the benefit of using machine learning techniques combined with physics-based models to efficiently address stochastic parameter optimization under parametric uncertainty. The developed framework is a useful tool that enables time-efficient stochastic evaluation of computationally expensive models and optimization of push-pull tests.