Pultruded fiber-reinforced polymer composites are susceptible to microstructural nonuniformity such as variability in fiber volume fraction (Vf), which can have a profound effect on process-induced residual stress. Until now, this effect of non-uniform Vf distribution has been hardly addressed in the process models. In the present study, we characterized the Vf distribution and accompanying nonuniformity in a unidirectional fiber-reinforced pultruded profile using optical light microscopy. The identified nonuniformity in Vf was subsequently implemented in a mesoscale thermal–chemical–mechanical process model, developed explicitly for the pultrusion process. In our process model, the constitutive material behavior was defined locally with respect to the corresponding fiber volume fraction value in different-sized representative volume elements. The effect of nonuniformity on the temperature and cure degree evolution, and residual stress was analyzed in depth. The results show that the nonuniformity in fiber volume fraction across the cross-section increased the absolute magnitude of the predicted residual stress, leading to a more scattered residual stress distribution. The observed Vf gradient promotes tensile residual stress at the core and compressive residual stress at the outer regions. Consequently, it is concluded that it is essential to take the effects of nonuniformity in fiber distribution into account for residual stress estimations, and the proposed numerical framework was found to be an efficient tool to study this aspect.