<p>Accurate and cost-effective surrogate modeling methods are essential for systems with stochastic response noise, where the primary challenge lies in effectively capturing the inherent uncertainty. Most existing approaches are derived from deterministic techniques and rely on expensive replications for noise estimation, with Stochastic Kriging (SK) drawing significant attention. The Nested SK (NSK) approach have addressed some of these issues by decoupling parameters in an iterative optimization process, but it is still hindered by inaccurate estimation in some complex cases. This paper introduces a novel surrogate modeling method for systems with stochastic noise, called arbitrary Polynomial Chaos Expansion Gauss-Seidel-based Stochastic Kriging (aPC-GSK), which enables accurate and robust surrogate modeling for stochastic simulations. Built on the Polynomial Chaos Kriging (PCK) framework, aPC-GSK introduces an adaptive arbitrary Polynomial Chaos Expansion (aPCE) to enhance global trend approximation without requiring prior distributional knowledge. Moreover, a Gauss-Seidel-based iterative optimization algorithm is developed to improve parameter estimation in SK modeling, addressing the limitations of NSK. Numerical experiments show that GSK significantly advances NSK in model fitting, while aPCE provides a more accurate global trend approximation. The combination of aPCE and GSK results in a highly accurate and robust surrogate model for stochastic simulation with both homoscedastic and heteroscedastic noise comparing with existing approaches. Practical case studies further validate the practicality of the proposed aPC-GSK, showcasing superior performance compared with other techniques in real-world stochastic simulations.</p>