The young beta Pictoris system has been monitored with high-resolution optical spectrographs for decades. These observations have revealed strongly variable, stochastic absorption in the CaII H&K lines attributed to in-falling cometary bodies. Since 2003, over 9000 HARPS observations of beta Pictoris have been taken and many of these have not yet been used for exocomet studies. We aim to search these spectra for new exocomet phenomenology enabled by the long time coverage and large volume of this dataset. We systematically carry out telluric correction of the HARPS spectra using molecfit, compare multi-year observations at the wavelengths of the CaII and NaI lines and use a Bayesian fitting algorithm to extract exocomet line parameters. We explore the usage of an unbiased reference spectrum with which to calibrate the continuum, and investigate Keplerian orbital solutions to observed exocomet acceleration. We find a general absence of exocometary sodium line absorption, with only two instances of clear ~2% deep exocometary sodium out of 198 nights of observation, as well as a weaker ~1% feature that persists over 13 nights in 2004. We find that these events occur during times of exceptionally deep CaII absorption, at the same red-shift, implying that strongly CaII-evaporating exocomets also exhibit detectable levels of Na I, in spite of the vast majority of NaI being rapidly photo-ionised in close proximity to the star. We find long-lived CaII absorption in 2017 and 2018 that persists on a timescale of a year, which may be difficult to explain with the classical exocomet model. Finally, we investigate two strongly accelerating, blue-shifted exocomet features observed in 2019 that show strong and sudden departures from Keplerian motion, suggesting rapid changes to the dynamics of the exocomet cloud. We hypothesize that this is caused by the destruction of the comet nuclei shortly after their periastron passages.