We present a model for upscaling the time-dependent effective retardation factor, R˜(t), in hierarchical porous media with multimodal reactive mineral facies. The model extends the approach by Deng et al. (Chemosphere 91(3):248–257, 2013) in which they expanded a Lagrangian-based stochastic theory presented by Rajaram (Adv Water Resour 20(4):217–230, 1997) in order to describe the scaling effect of R˜(t). They used a first-order linear approximation in deriving their model to make the derivation tractable. Importantly, the linear approximation is known to be valid only to variances of 0.2. In this note, we show that the model can be derived with a higher-order approximation, which allows for representing variances from 0.2 to 1.0. We present the derivation and use the resulting model to recalculate R˜(t) for the scenario examined by Deng et al. (2013).