Using GAMs to model fish density data with many "zeroes" I'm trying to fit GAM models to acoustic backscatter data (Sv_mean; measured in dB, log scale) which is proportional to the density of fish in the water column. I have four years of daily data. In the winter time there are a bunch of zeros (which in log scale are represented as -999dB). The zeroes are true zeroes. To go around having these -999, I replaced the zeroes with values 3 dB below the minimum backscatter value in the series (~ -134 dB that's roughly half the minimum fish density).
To account for autocorrelation in the data I'm including a correlation term in the gamm() function (my data is ordered and I don't have NAs). I'm also including a linear form of time (date_numeric) to account for interannual patterns and day of year (julian_day) using a cyclic cubic spline to look at seasonal patterns. I am planning to do model selection using different candidate models with a different set of environmental covariates.
I've fitted a couple of GAM models and it seems that no model is picking up the transformed zeroes. Because the zeroes have a seasonality, the normalized residuals keep having some autocorrelation.
What are possible solutions to this?
Here is the time series:

Here is how the distribution of the backscatter values (Sv_mean) look like:

Here is the example model I'm fitting:
gam_6 <- gamm(Sv_mean ~ s(bottom_temperature1, bs="cr") +  
    s(mw_salinity2, bs="cr") + s(sea_ice, bs="cr") + 
    s(nitrate_corrected, bs="cr") + s(julian_day, bs="cc", k=365)
    + s(mw_PAR2, bs="cr") + s(chla_insitu_corrected, bs="cr") + 
    s(date_numeric, bs="cr"), method="REML", data=data, 
    correlation = corAR1(form = ~ 1))

Here is the output from gam.check()

Here is the data and fitted model (you can see the transformed zeroes occurring during specific months each year):

And here are the normalized residuals with their ACF and PACF:


You could try a zero adjusted gamma (ZAGA) distribution model in the gamlss package.
ZAGA(mu, sigma,nu) distribution
takes value 0 with probability nu,
and a gamma GA(mu, sigma) distribution
with probability (1-nu).
All the ZAGA parameters, mu, sigma and nu, can be modelled using explanatory variables.
