Unsteady solute dispersion in pulsatile non-Newtonian flow in a circular tube has been investigated considering higher order moments using the Aris’s method of moments. The Luo and Kuang(1992) (K - L Model) constitutive relation is used for the non-Newtonian fluid. The axial mean concentration of solute is estimated by considering exchange, convective, dispersion coefficients. In addition, considering the 3rd and 4th central moments, the skewness, and kurtosis in the mean concentration is examined. Variations of these five moments against the absorption parameter, K-L fluid parameters, and Womersley frequency parameter are thoroughly investigated. The increase in the Womersley frequency parameter led to the occurrence of a double frequency period for the dispersion coefficient, which influenced the skewness and kurtosis coefficients that led to the deviation from usual normal distribution of the mean concentration distribution at the starting of the dispersion process. The concentration distribution is skewed towards the central region of the tube. The concentration contours also reveal that the solute convection is more in the central region of the tube. The results for Newtonian, Bingham, and Casson fluid models are retrieved as special cases. The mean concentration has the highest peak for the K-L model and the lowest for the Newtonian fluid model. It is scattered axially more for the Bingham fluid than the Casson fluid. This investigation emphasizes that higher moments provide precise information about the solute dispersion in the flow field and its deviation from the normal distribution at small time scales. • This investigation examining the unsteady axial dispersion of a solute in the non-Newtonian pulsatile flow in a cylindrical tube, considering the impact of an irreversible first-order boundary reaction. Solute dispersion in K-L model has been introduced, which gives the shear thinning and thickening behaviour of blood flow even if yield stress is zero. • The present study initiated skewness and kurtosis analysis for steady/ unsteady, non-pulsatile/ pulsatile and Newtonian/ non-Newtonian fluid flow by using method of moments. Here analytical results for all transport coefficient with skewness and kurtosis has been discussed. • The perturbation and numerical solution of coupled and non-linear equation of shear stress and momentum equation has been obtained to get the velocity profile. A new class of computationally explicit Runge–Kutta (CERK) numerical method (Yadav et. al. (2022)) has been used to solve this stiff problem. A detailed comparison of perturbation and numerical solutions is also given in this manuscript. • The solution of the concentration equation is obtained analytically adopting the Aris0 s method of moments, with fourth-order precision. In the presence of boundary reaction, the transport coefficients are affected not only by time but also by the Womersley frequency parameter, other chemical variables, the wall absorption parameter, the yield stress, the amplitude of the fluctuating pressure component, and plasma viscosity. The non-Gaussian distribution of solute is visible for small time after injection. • The mean concentration as well as the two-dimensional concentration distribution are accurately described. The results for the Bingham model, the Casson model, and the Newtonian model are derived directly from this analysis.