Sound wave propagation in a shallow water environment is complex due to e.g. the uncertainties of sound speed profile being inhomogeneous and imprecisely measured, the bottom reflections, etc. The propagation and influence of several uncertainty parameters are quantified in this paper. A four-layer model, which can approximately represent a wide range of shallow water environments, is considered; six parameters representing sound speed profile and water depth are considered as random variables. We investigate how the wave field (pressure) in this model is influenced by these uncertainties. For this purpose, the sound field is computed for different realizations of the random variables, when the medium is excited with sources whose frequencies are appropriate, for example, for marine seismic exploration applications. Since classical Monte Carlo methods require a huge sample size to converge, we use three surrogate modeling techniques (Kriging, Polynomial Chaos, and Polynomial Chaos-based Kriging). The proposed methods require much smaller sample sizes, which makes the uncertainty quantification (UQ) possible. Wavelength-to-depth ratio (lambda/d) is introduced as the key parameter that defines the degree of interaction (reflection and transmission) of the sound waves with the boundaries of the shallow water waveguide. The results show that for small and large values of lambda/d, the wave field is more sensitive to the variations of the water depth and the velocity of the bottom layer, respectively. The robustness (precision) of the surrogate models is shown to decrease for lower values of lambda/d. The proposed UQ methodology can be used for more complicated underwater environments; it is even more advantageous because it can efficiently deal with a large number of model uncertainty parameters and identify the most influential ones.