A non-linear lumped kinetic model (LKM) of liquid chromatography formulated for simulating the separation of multi-component mixtures using gradient elution chromatography in an overloaded column. The model incorporates parameters such as Henry constants, nonlinearity coefficients, mass-transfer coefficients, and axial-dispersion coefficients, which vary with solvent strength using the Linear solvent strength (LSS) gradient approach. Numerical solutions are obtained using a high-resolution finite volume method with a flux limiter to avoid numerical dissipations and resolve sharp discontinuities to achieve higher-order accuracy. Investigates the influence of modulator concentration variations on one, two, and three-component mixtures and compares the results with isocratic elution under fixed conditions. Additionally, analyze the effects of gradient slope, modulator concentration, and other parameters on elution profiles and the production of targeted components in multi-component gradient elution. This model relates to the probability principle and employs numerical temporal moments to enhance graphical illustrations of results during the running phase, offering insights into liquid chromatography behavior and optimization of the separation performance, sample retention, and band broadening. The proposed theoretical study for this model provides an efficient method for optimizing gradient elution chromatography processes and has implications for process optimization.