[1] The probabilistic analysis by Monte Carlo simulation method (MCSM) for the transport of nonlinear reactive solutes in a three-dimensional heterogeneous porous medium is a computationally prohibitive task. For linear transport problems, the perturbation-based stochastic finite element method (SFEM) has been found to be computationally efficient with acceptable accuracy. This provides a motivation to develop the SFEM for the nonlinear reactive solute transport. In the present study, SFEM is developed for the transport of equilibrium nonlinear sorbing solutes, which follow the Langmuir-Freundlich isotherm. This method produces a second-order accurate mean and a first-order accurate standard deviation of concentration. In this study, the governing medium propertis viz. hydraulic conductivity, dispersivity, molecular diffusion, porosity, sorption, and decay coefficients are considered to vary randomly in space. The performance of SFEM is compared to MCSM for both one-and three-dimensional transport problems. The mean and the standard deviation of concentration for various test cases obtained with the SFEM compares well for the mild heterogeneity cases (standard deviation of log hydraulic conductivity less than 0.85) tested. SFEM produces a sharp front for the mean and the standard deviation of concentration while fronts obtained by MCSM are found to be dispersive. The error associated with the results obtained by SFEM is sensitive to the boundary conditions, the size of the domain, and the plume size. For a higher nonlinearity of sorption isotherm, the prediction uncertainty is higher. The pattern of the statistical moments of concentration is similar for cases with different correlation lengths of the parameters.