Analysis of incomplete block design and block-level attributes

Question

I am working in a system where experimental treatments are necessarily associated with specific blocks (a complete design is not possible) and am looking for guidance on valid model specification.

Basic background: I have data on 4 lakes. Each year I measure fish in each lake, repeated for ~10 years. Each lake has several fixed characteristics (stratified/unstratified, depth) as well as attributes that are specific to a lake but change each year (temperature degree days). Body weight is the response variable, and individual fish also have attributes like sex. I am particularly interested in how "temperature degree days" and the "stratified/unstratified" characteristic of a lake influences my response variable.

The issue I am grappling with is how to account for clustering within lakes while also modeling the effect of lake attributes themselves (the situation is exacerbated by the fact that there are only 4 lakes, so treating lake as a "random" grouping variable in multilevel model is tenuous given the presence of only 4 levels). I'd like to determine what the independent effects of my two lake level properties are on my response, but also avoid pseudoreplication by considering clustering within lakes -- however including multiple lake characteristics makes each lake fully identifiable, and so a corresponding model is rank deficient (a specific combination of stratified and depth is enough to know which lake is associated).

My preferred analysis program is R, and modeling with the package lme4. To give an idea of the problem, here is some sample data:

set.seed(1001)
dat <- data.frame(body_weight = rnorm(2400, 20, 2), 
    lake= rep(c("A", "B", "C", "D"), each=600), 
    depth= rep(c(100,150,400,550), each=600), 
    strat = rep(c("stratified", "unstratified"), each=1200), 
    sex= sample(c("M", "F"), replace=TRUE, size = 2400, 
                prob=c(0.5,0.5)), 
    year=factor(rep(rep(c(2010, 2012, 2014, 2016, 2018, 2020), 
                    each=100), 4)), 
    degree_days = rep(rnorm(24, mean=150, sd = 40), each=100))
dat$year_numeric <- as.numeric(as.character(dat$year))

I naively started with this specification:

library(lme4)

mod <- lmer(body_weight~ lake + depth + strat + sex + year_numeric + 
             degree_days + (1|year), data=dat)  

But of course this was rank deficient because "lake" is perfectly identified by my "depth" and "strat" variables. Dropping "lake" the model doesn't appear to be overtly problematic:

mod2 <- lmer(body_weight~ depth + strat + sex + year_numeric + 
              degree_days + (1|year), data=dat)       

However, I am surprised to see that specifying "lake" as part of the random effects structure does not cause singularity problems. I had originally sought to model lake as a fixed effect because it only has 4 levels, but mostly I'm just surprised that the model converges with "depth", "strat" and "lake" all included (any 2 can tell you what the third should be), as long as "lake" is a random effect. These 2 models for example:

mod3 <- lmer(body_weight~ depth + strat + sex + year_numeric + 
    degree_days + (1|year) + (1|lake), data=dat)         
mod4 <- lmer(body_weight~ depth + strat + sex + year_numeric + 
    degree_days + (1|lake/year), data=dat)

So I am wondering, are these valid model specifications for attempting to assess the effect of "degree_days" and "strat" on "body_weight"? Or am I in trouble because my continuous predictors of interest are properties of a shared grouping variable?

And if I can't include "lake" in the model, am I at risk of failing to account for real clustering (pseudoreplication) despite accounting for specific attributes that define the clusters?

So far I have followed advice from: the comments in (https://biologyforfun.wordpress.com/2015/08/31/two-little-annoying-stats-detail/#comment-324) and How to analyze this incomplete block design in R? Incomplete block design analysis / design with R