The computation of the mean stress in coupled poroelasticity problems discretized with low-order finite elements typically gives rise to unphysical oscillations. In contrast to pore pressure instabilities due to the violation of the inf-sup stability condition, which naturally vanish in time, the mean stress oscillations are persistent throughout the whole simulation, unless deliberately removed. Stabilization is obtained by: (i) adding the mean stress as a main unknown by a mixed discretization, and (ii) enriching the finite element approximation with the scaled Laplacian of displacements, which is evaluated by the Physical Influence Scheme (PIS). It can be proved that the proposed approach is equivalent to introduce a minimization constraint on the gradient of the mean stress function. The resulting discrete algebraic problem is characterized by a double saddle-point matrix, which requires special preconditioning tools to accelerate convergence with Krylov subspace methods. A key step is the calculation of the scaling stabilization parameter ℎ by solving an optimization problem at the elemental level. To avoid the high computational cost associated with this procedure, a machine learning model is proposed to predict the optimal ℎ as a function of the elemental shape and size. The benefit of combining PIS with a machine-learning model for the computation of ℎ is that the stabilization scheme does not rely on any type of heuristic or user-specified tuning parameter, as often required by other stabilization methods. The results show that the proposed stabilization strategy can effectively remove mean stress oscillations in coupled poroelasticity. The calculation of ℎ is shown to be critical for the quality of the stabilization, with the machine learning-based approach providing an otpimal compromise between numerical diffusion and solution accuracy.