The paper demonstrates the use of the response surface method (RSM) to carry out probabilistic assessment of the bearing capacity of shallow footings seated on geosynthetic reinforced granular layers over randomly variable soft clay deposits. The method substantially reduces the number of Monte Carlo simulations required to carry out the probabilistic finite element limit analyses of the bearing capacity problem. A finite element limit analysis model based on the lower bound theorem is furnished and verified using some well-known analytical methods, and is then used to generate a large synthetic database of numerical results for bearing capacity of shallow foundations on a reinforced granular fill over randomly variably soft clay deposits. To this end, a permutation of the important influencing parameters is formed, and lower bound FELA-based limit loads are sought through optimization in MATLAB. A closed-form solution is formulated using RSM-based polynomials. The RSM equations, which are acquired from least squares regression analyses, are used to carry out probabilistic Monte Carlo simulations, and the results are presented in forms of cumulative distribution functions. Results from the probabilistic analyses are introduced into reliability-based design approach to render design loads for different reliability levels.