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 termbased on a Partial Differential Equation, that embeds physical knowledge on the phenomenon under study.