The simultaneous estimation of multiple quantiles is essential for analyzing phenomena with complex distributional characteristics, which standard parametric assumptions can not describe, and where the interest lies in tail behavior, rather than in the mean. In this work in particular, we are interested in estimating multiple quantile surfaces, which is crucial when investigating environmental phenomena. For instance, in the analysis of particulate matter (PM10) concentrations, researchers are particularly interested in tail events, which have severe consequences for human health. Estimating simultaneously multiple quantile levels is crucial for preserving monotonicity, which ensures consistent and reliable estimates. Moreover, the resulting monotone quantile surfaces allow for a fully nonparametric reconstruction of the probability density function of the variable of interest, at any location of the domain. This work introduces a novel methodology, based on a penalized functional approach, that permits to estimate simultaneously multiple quantile surfaces. The method is developed within the framework of physics-informed models, where the minimization functional includes a regularization term based on a Partial Differential Equation, that embeds physical knowledge on the phenomenon under study. In the case of PM10 analysis, this enables for instance the integration in the statistical model of information concerning wind stream over the region, a critical factor in the dispersion of air pollutants. A second regularization term is used to enforce monotonicity of the quantiles. The resulting monotone quantile surfaces facilitates the development of data analysis tools, such as spatial exceedance probability maps, which may provide policymakers with robust insights to inform strategies for mitigating the harmful effects of air pollution.