How to correct quantile deviations in glmmTMB with zero inflated, over distributed data found using the DHARMa package

Question

I am analyzing data on bird landing occurrence on ponds in proximity to the river, and whether birds with certain traits more frequently land in certain places.

I have a data frame with 19,447 observations. The observations go up to 2022 and there are 12 different ponds so the size_ha and River_prox values are repeated quite often. Each line is a summary of the observations on that day.

My data is count data that is heavily zero inflated and over dispersed. I have one pond with about 10 fold more landings than the others; also the detection frequency has increased over the years so that likely contributes to the dispersion issue as well. Therefore I am using the glmmTMB package and I have modeled with both poisson and negbin2 families to see which is better (I think negbin because of the dispersion).

The model contains the offset of pond size, since we know the bigger the pond the more birds. I also have the random effect terms of Pond_ID, Year and Day of Year (DOY). I am not sure whether I need all of these and what the best format to include them would be but I have tried many. From what I understand, Pond_ID is nested within Year since the pond sizes have changed throughout the years? Any input on that would be appreciated.

My models are as follows: Poisson:

model_full <- glmmTMB(Landed ~ Guild * River_prox + 
   offset(log(Size_ha))+(1|DOY) +  (1|Year/Pond_ID),
   ziformula= ~Guild,data=all_landings_sum_guild, 
   family="poisson")

summary()
      AIC       BIC    logLik  deviance  df.resid 
 295482.5  295646.5 -147717.2  295434.5      6855 

Negbin with a different random effect structure since I got errors with others :

model <- glmmTMB(formula = Landed ~ offset(log(Size_ha)) + 
   Guild * River_prox + 
    (1 + (1 |         Year/Pond_ID) + (1 | DOY)),
    family = nbinom2, data = all_landings_sum_guild)

I have tried multiple random effect structures and constantly get warnings and errors that I do not understand. Further, all of my residual plots look bad and everything I try from the forums either gives me errors, or does not fix the issue.

Error examples with (1|DOY)+(1|Year/Pond_ID) structure: which gets these errors :

Warning messages:
1: In finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
2: In finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old) : Model convergence problem; function evaluation limit reached without convergence (9). See vignette('troubleshooting'), help('diagnose')

I have read the trouble shooting page and tried to fix it but changing the error structure takes away the warning. The plots look essentially the same no matter what I do.

Here is the summary of the negative binomial model

summary(model)
 AIC      BIC   logLik deviance df.resid 
 50610.9  50734.0 -25287.5  50574.9     6861 

I am checking the fit using the following outputs. All of the included results are done on the negative binomial model because I think it is the best one thus far.

simulationOutput <- simulateResiduals(fittedModel = model, plot = FALSE)
plot(simulationOutput)
testQuantiles(simulationOutput)

Result:

data:  simulationOutput
p-value < 2.2e-16
alternative hypothesis: both

testZeroInflation(simulationOutput)

Results:

data:  simulationOutput
ratioObsSim = 0.042636, p-value < 2.2e-16
alternative hypothesis: two.sided

testOutliers(simulationOutput)

Result:

outliers at both margin(s) = 72, observations = 6879, p-value =
0.02487
alternative hypothesis: true probability of success is not equal to 0.007968127
95 percent confidence interval:
 0.008198258 0.013163080
sample estimates:
frequency of outliers (expected: 0.00796812749003984 ) 

                                        0.01046664

testDispersion(simulationOutput)

data:  simulationOutput
dispersion = 1.294, p-value = 0.4
alternative hypothesis: two.sided

I also did a bootstrap

outliers at both margin(s) = 71, observations = 6879, p-value =
/sup>0.22

alternative hypothesis: two.sided I have read the forums and the DHARMa and glmmTMB troubleshooting guides but I either do not understand what its asking me to do, it gives me an error, or it does not fix it. I think perhaps there is something about my data I am not understanding.