Many ecological insights into the function of rivers and watersheds emerge from quantifying the flux of solutes or suspended materials in rivers. Numerous methods for flux estimation have been described, and each has its strengths and weaknesses. Currently, the largest practical challenges in flux estimation are to select among these methods and to implement or apply whichever method is chosen. To ease this process of method selection and application, we have written an R software package called loadflex that implements several of the most popular methods for flux estimation, including regressions, interpolations, and the special case of interpolation known as the period‐weighted approach. Our package also implements a lesser‐known and empirically promising approach called the “composite method,” to which we have added an algorithm for estimating prediction uncertainty. Here we describe the structure and key features of loadflex, with a special emphasis on the rationale and details of our composite method implementation. We then demonstrate the use of loadflex by fitting four different models to nitrate data from the Lamprey River in southeastern New Hampshire, where two large floods in 2006–2007 are hypothesized to have driven a long‐term shift in nitrate concentrations and fluxes from the watershed. The models each give believable estimates, and yet they yield different answers for whether and how the floods altered nitrate loads. In general, the best modeling approach for each new dataset will depend on the specific site and solute of interest, and researchers need to make an informed choice among the many possible models. Our package addresses this need by making it simple to apply and compare multiple load estimation models, ultimately allowing researchers to estimate riverine concentrations and fluxes with greater ease and accuracy.