A viscoelastic model for cellulose nanofiber (CNF) suspensions using a population balance equation for fiber floc aggregation-breakage was proposed. A Krieger-Dougherty model was used to describe the dependence of viscosity on the volume fraction of flocs, and a White-Metzner type model was used to represent the effect of viscoelasticity. The rheological properties of the present model were investigated, and this model was confirmed to describe shear thinning in viscosity and viscoelastic responses of stress growth, depending on the relaxation time. Furthermore, flows of CNF suspensions through a circular tube were analyzed using the present model. The velocity profile changes with time according to the temporal change in the distribution of the effective volume fraction of fiber flocs. The numerical prediction of the velocity profile captures the typical properties of shear-thinning fluids. The development rate of the floc size distribution depends on the velocity gradient and affects the temporal changes of both velocity and stress fields. The development of the stress field is delayed longer for a longer relaxation time. These phenomena affect the temporal change in the floc size distribution and hence also affect the growth of the velocity field through the effective volume fraction of flocs. These results indicated that the present model can describe the viscoelastic behavior of CNF suspensions in flows.