A global dataset which is composed of more than 20,000 records is used to develop an empirical nonlinear soil amplification model for crustal earthquakes. The model also includes the deep soil effect. The soil nonlinearity is formulated in terms of input rock motion and soil stiffness. The input rock motion is defined by the pseudo-spectral acceleration at rock site condition (PSA(rock)) which is also modified with between-event residual. Application of PSA(rock) simplifies the usage of the site model by diminishing the need of using the period-dependent correlation coefficients in hazard studies. The soil stiffness is expressed by a Gompertz sigmoid function which restricts the nonlinear effects at both of the very soft soil sites and very stiff soil sites. In order to surpass the effect of low magnitude and long-distant recordings on soil nonlinearity, the nonlinear site coefficients are constrained by using a limited dataset. The coefficients of linear site scaling and deep soil effect are obtained with the full database. The period average of site-variability is found to be 0.43. The sigma decreases with decreasing the soil stiffness or increasing input rock motion. After employing residual analysis, the region-dependent correction coefficients for linear site scaling are also obtained.