Implications of keeping a "low" basis dimension in GAMM

Firstly, I think you could simplify your model a bit as I don't see the point of both of these terms:

1. s(fSite, bs = "re"), and
2. s(fSite, by = fSeason, bs = "re")

Both will estimate a separate mean response for each site, with the former penalizing all sites the same regardless of season, while the latter will do the same thing but allow for different levels of penalization between wet a dry seasons. If you think sites have different random abundances between seasons, you could just use s(fSite:fSeason, bs = "re"), whih will give you a random intercept for each site by season combination.

As to the thrust of your main question...

The test for basis size is a heuristic guide, not some infallible hard-and-fast rule that you must obey. It works by comparing adjacent residuals ordered by the covariate of the focal smooth, computing some metric that measure how similar adjacent residuals are, and then computing the same metric on residuals whose order has been permuted.

When you have data that are ordered in time, and especially when you are smoothing by time, as you are here, that heuristic test is likely to be confused by autocorrelation. Indeed, wiggly smooths and strong autocorrelation are often indistinguishable mathematically.

So, the interpretation of your k.check() output could just as plausibly be the result of unmodelled autocorrelation.

I know we're supposed to increase the basis dimension if k.check indicates it's too small...

You don't have to; the test is just a guide which can indicate if there is potentially unmodelled wiggliness.

...but I think a smoother pattern is fine since my main goal is prediction.

This may work, but you are going to induce bias into your estimates if that unmodelled wiggliness is signal that you are interested in rather than noise that you aren't. And if you are doing any inference on the estimated smooths, using credible intervals, etc, these are going to biased and anti-conservative if your residuals are not conditionally i.i.d.

Can I ignore the low p-values and keep the basis dimension, k, at the default?

You can decide not to increase k, but you can't ignore the potential for a problem: you need to understand why the test rejected the null hypothesis. As I mentioned above, this could be due to unmodelled autocorrelation, but it could be due to too small a basis for the temporal components. And those two things could be the same thing.

Also, is there an issue with "using up" all the available EDF to max-out k (i.e. does setting k=15, the max, come with any issues)?

Not really; in general the main downside to fitting with larger k is the increased compute time. You can estimate implausibly wiggly effects of covariates if k is too large, so you can use k as a prior on the upper limit of expected wiggliness. But that is usually an issue with covariates that are not time. I don't think there is a a prior expectation that smooths of time should not be wiggly.

A consideration with time however is how you decompose that temporal effect. You could put it all into a very wiggly smooth (if you have the data), or you could fit a very smooth trend with the short-term dependence modelled using a correlation process in the covariance matrix.

It might also be that your trend shouldn't be smooth; the trend might be better off consider stochastic, or smooth in the sense of a random walk (via an MRF smooth).

I'm afraid of over-fitting the model.

I don't think you need to worry too much about this. When you are modelling time a wiggly trend could be just fine. As I mentioned above, it really depends how you want to model/decompose time. Do you want a simpler smooth trend? If so you can keep k low, but you need to model the autocorrelation some other way. See below.

Some further thoughts, observations

You should also be plotting residuals against predictors to understand how well the temporal smooths capture the underlying trends in the data. Or plot the observed counts and model predictions to see the same thing. You should be looking at whether you are adequately modelling your data and also trying to understand why the smooths might be being estimated as wiggly as they are. Are those oscillations in s(CYR.std) for the wet seasons justified by the data or is the smooth chasing odd years?

An alternative approach to detect unmodelled wiggliness is to fit the model

r <- resid(m)
m_2 <- gam(r ~ s(sal) + s(water_depth) + fSeason + s(CYR.std, by = fSeason) + 
    s(fSite, bs = "re") + s(fSite, by = fSeason, bs = "re"),
  method = "REML",
  family = quasi(link = "identity", var = "constant"),
  data = shrimp2)

where we are fitting the deviance residuals instead of the original response, so we can drop the offset term. We also change the family because we expect these to be distributed mean 0 with constant variance, but we aren't assuming that they have some specific conditional distribution, such as being conditionally Gaussian.

If you can, increase k on the smooths in that model above the k used originally. Also, adjust the model formulation depending on what you end up doing after considering my initial observation at the start of this answer.

When you look at summary(m_2) you should be looking to see if there is any smooth that uses more than 1 effective degree of freedom, and if there is unmodelled wiggliness. Any that do are candidates but having their k increased, when you return to the model for the observed response.

But even this model test may erroneously suggest you have too small initial basis sizes. Again because of autocorrelation. While it is OK to model this autocorrelation by smooths of time, you might prefer a different decomposition; where you include a correlation among residuals to model the short-scale temporal dependence. With the response and family you are trying to fit, this is going to be quite difficult to do in practice.

One simple GEE-like approach, although you need to provide $\rho$, would be to fit using bam() and provide an estimate of $\rho$ to argument rho. To get this estimate, you'll need to use the ACF and that assumes your data are regularly spaced in time. The trick here will be to subset the model residuals into individual time series, by site, within year, I think by the way you are decomposing time. This will yield many ACFs, each yielding an estimate of $\rho$, and you'll need to provide some average of these to rho, plus tell AR.start where each time series starts/ends as per ?bam. Note that this will fit the same AR(1) process in each site/year series, not a separate AR(1) for each, so you're only estimating a single $\rho$.

To judge the success of this, you'll then want to plot the standardized residuals, which are in the $std.rsd (or fit an ACF to those residuals, again one per time series).

If you need anything more complex that this, you might need to in practice try brms::brm() as that can fit smooths and more complex autocorrelation processes.