Wastewater treatment plants (WWTPs) face operational challenges due to variable inflows, changing wastewater composition, and fluctuating energy prices and availability, including the risk of energy-disruptions as blackouts. To enhance resilience, we develop a grey-box probabilistic model for forecasting wastewater inflow and chemical oxygen demand (COD) within a digital twin framework. The model relies on Generalized Additive Models for Location, Scale, and Shape (GAMLSS) with Johnson SU distribution to quantify uncertainty via confidence intervals. It integrates real-time WWTP data, weather observations, and network-wide monitoring inputs. Non-linear patterns such as daily, weekly, and seasonal cycles are captured using quadratic splines and ReLU activation functions. Trained on time series data from a WWTP in southern Germany, the model provides probabilistic forecasts over a flexible time horizon. Unlike black-box approaches, our method maintains interpretability while accurately modeling complex dynamics. A key innovation is the extension of GAMLSS to predict COD composition, not just inflow. The model also supports scalability: preliminary feature selection across monitoring stations identifies key predictors and can inform optimal sensor placement. By combining probabilistic forecasting, interpretable modeling, and data-driven optimization, this framework enables adaptive, energy-efficient operations and strengthens the resilience of critical water infrastructure.