Damage from severe ground motion occurred in the downtown area of Mashiki in the Kumamoto Prefecture during the 2016 Kumamoto earthquake, Japan; such damage was heavy in the center of the downtown area. Nonlinear site responses for the first shock and mainshock, which occurred on 14 and 16 April 2016, respectively, are important factors that explain why the area was heavily damaged. We analyzed soil nonlinearity using surface and borehole records obtained from the KiK‐net KMMH16 (Mashiki) station. From our analysis, we found that S‐wave velocity models clearly depended on the amplitude of input ground motion. We estimated the strain‐dependent shear stiffness and damping ratio to explain this S‐wave velocity dependence. We conducted equivalent linear analyses at the KMMH16 site, based on a nonlinear model. From these analyses, we concluded that our synthetic surface ground motions agreed well with the observed ones, especially for the S‐wave amplitudes and phases of the first shock and mainshock noted above. In addition, we performed the same analyses at the TMP3 site, which was actually located within one of the heavily damaged zones. The synthetic motions here also agreed with the observed ones, with differences in spectral accelerations being well explained by our analyses. Our results indicated that soil nonlinearity played a major role in causing the difference of ground motions, thus leading to the heavily damaged zone in the downtown area of Mashiki.