Computational walking simulations could facilitate the development of improved treatments for clinical conditions affecting walking ability. Since an effective treatment is likely to change a patient's foot-ground contact pattern and timing, such simulations should ideally utilize deformable foot-ground contact models tailored to the patient's foot anatomy and footwear. However, no study has reported a deformable modeling approach that can reproduce all six ground reaction quantities (expressed as three reaction force components, two center of pressure (CoP) coordinates, and a free reaction moment) for an individual subject during walking. This study proposes such an approach for use in predictive optimizations of walking. To minimize complexity, we modeled each foot as two rigid segments—a hindfoot (HF) segment and a forefoot (FF) segment—connected by a pin joint representing the toes flexion–extension axis. Ground reaction forces (GRFs) and moments acting on each segment were generated by a grid of linear springs with nonlinear damping and Coulomb friction spread across the bottom of each segment. The stiffness and damping of each spring and common friction parameter values for all springs were calibrated for both feet simultaneously via a novel three-stage optimization process that used motion capture and ground reaction data collected from a single walking trial. The sequential three-stage process involved matching (1) the vertical force component, (2) all three force components, and finally (3) all six ground reaction quantities. The calibrated model was tested using four additional walking trials excluded from calibration. With only small changes in input kinematics, the calibrated model reproduced all six ground reaction quantities closely (root mean square (RMS) errors less than 13 N for all three forces, 25 mm for anterior–posterior (AP) CoP, 8 mm for medial–lateral (ML) CoP, and 2 N·m for the free moment) for both feet in all walking trials. The largest errors in AP CoP occurred at the beginning and end of stance phase when the vertical ground reaction force (vGRF) was small. Subject-specific deformable foot-ground contact models created using this approach should enable changes in foot-ground contact pattern to be predicted accurately by gait optimization studies, which may lead to improvements in personalized rehabilitation medicine.