A new model for polymer silylation is proposed that takes into account three mechanisms : The relaxation of the polymer during silylation, a diffusivity of the silylating agent which increases with resist matrix expansion, as well as a local reaction rate retardation due to stress induced by resist swelling. These mechanisms capture changes in profile shape during silylation. The resulting nonlinear differential equations are solved two-dimensional (2D) using the finite element method and the influence of the main parameters is studied as a function of feature size and feature type, illustrating the effects of initial aerial image and processing conditions on the final resist image. The simulation profiles reproduce the silylation depth reduction with feature size and abrupt encroachment heretofore unexplained in various experimental studies.