Dear @RenQiang.Xi,
When SoilDyn is coupled to SubDyn, the linear part of the soil reaction (Ks in your nomenclature) is considered on the left-hand side of the SubDyn equations of motion (and influences the Craig-Bampton modes), and the nonlinear part of the soil reaction is applied as a force on the right-hand side of the SubDyn equations of motion. In this way, the effect of the linear reaction is not double counted. We had originally tried to include the full (linear + nonlinear) soil reaction on the right-hand side as you suggest, but the coupled solution was numerically unstable.
We had some discussion with the authors of the paper you cite and their FounDyn module, hoping to get features of FounDyn merged into SoilDyn, but that merging has not yet happened as far as I know.
Best regards,