Terrestrial water availability, quantified by precipitation minus evapotranspiration (P-E), is essential in Earth's water cycle, whereas model simulation of P-E is still largely biased and requires a post-processing procedure. This study introduces the grid-by-grid cumulative distribution function (CDF) matching method to correct simulation bias in P-E, based on the ERA5-Land dataset and outputs from 13 selected CMIP6 global climate models. The CDF matching method has a particular advantage in preserving the trends simulated by laws of physics in climate models, and three (additive, multiplicative, and additive-multiplicative mixed) trend preservation strategies are compared in this study. The cross-validation from 1951 to 2014 indicates that all the trend preservation strategies effectively improve the simulated spatial characteristics of P-E with increased spatial correlation, enhanced sign agreement and reduced mean absolute error. Specifically, the additive strategy outperforms in improving the spatial similarity and accuracy of P-E in the humid region and global average, whereas the mixed strategy is the optimal in the hyper-arid, arid, and semi-arid regions. Furthermore, the mixed strategy has a significant advantage in preserving the signs of P-E across the globe. This study exhibits a computationally efficient statistical approach for bias correction of P-E simulation, and validates its flexible correction strategies regarding different terrestrial aridity conditions.