This paper presents a quantitative validation of Scrape-Off Layer Plasma Simulation-ITER (SOLPS-ITER) simulations against the TCV-X21 reference case and provides insights into the neutral dynamics and ionization source distribution in this scenario. TCV-X21 is a well-diagnosed diverted L-mode sheath-limited plasma scenario in both toroidal field directions, designed specifically for the validation of turbulence codes (Oliveira, Body et al 2022 Nucl. Fusion 62 096001). Five new, neutrals-related observables are added here to the extensive, publicly available TCV-X21 dataset. These are three deuterium Balmer lines in the divertor and neutral pressure measurements in the common and private flux regions. The quantitative agreement metric used in the validation is combined with the conjugate gradient method to approach the SOLPS-ITER input parameters that return the best overall agreement with the experiment. A proof-of-principle test of this method results in a modest improvement in the level-of-agreement; the shortcomings impacting the result and how to improve the methodology are discussed. Alternatively, a scan of the particle and heat diffusion coefficients shows an improvement of 10.4% in the level-of-agreement, approximately twice as high as that achieved by the gradient method. This result is found for an increased transport coefficient compared to what is usually used for TCV L-mode plasmas. The simulations further indicate that similar to 65% of the total ionization occurs in the SOL, from which similar to 70% in the divertor regions, despite being a sheath-limited regime, motivating the inclusion of self-consistent neutral models in future turbulence simulations on the path towards improved agreement with the experiment.