R: Mixed ANOVA works with the aov{stats} function and fails with the anova_test{rstatix} function. Why?

I can implement a rather complicated mixed ANOVA in R with the aov{stats} or the lmer{lme4} functions, but not with the anova_test{rstatix} nor the Anova{car} functions. My question follows: How can I make the anova_test{rstatix} or Anova{car} functions work to access Greenhouse-Geisser corrections?

I have a factorial, incomplete blocks experiment with repeated measures. I replicated factorial nitrogen (N), phosphorus (P) and potassium (K) addition four times, blocking control, +NP, +NK and +PK treatments versus +N, +P, +K and +NPK treatments within each replicate. I do not evaluate the three-way N:P:K interaction because blocking reduces degrees of freedom. I determined (and standardized) soil P concentration (P.conc) at five depths. Depth is the repeated measures variable. I retain three depths to illustrate my problem. The data follow:

z = data.frame(plot=c(1:7, 9:15, 17:24, 26:30, 32:36))
z$P = factor(ifelse(z$plot %in% c(1,3,5,7,13:15,17,20:22,24,29,30,33,35), 1, 0))
z$K = factor(ifelse(z$plot %in% c(2,5,7,9,10,13,17:21,27,29,32:34), 1, 0))
z$N = factor(ifelse(z$plot %in% c(3:5,9:11,15,17,20,22,23,27:29,34,35), 1, 0))
z$block = factor(ifelse(z$plot %in% c(1,2,4,5,9,12,13,15,19,20,23,24,28:30,32), 1, 2))
z$rep = ifelse(z$plot %in% c(1:7,10), 1, 0)
z$rep = ifelse(z$plot %in% c(9,11:15,17,18), 2, z$rep)
z$rep = ifelse(z$plot %in% c(19:24,26,27), 3, z$rep)
z$rep = factor(ifelse(z$plot %in% c(28:30,32:36), 4, z$rep))
z$plot = factor(z$plot)
z = rbind(z,z,z)
z$depth = factor(rep(c("0-5 cm", "10-20 cm", "50-100 cm"), each=32))
z$P.conc = c(1.59, -0.34, 1.78, -0.32, 1.76, -0.29, 3.42, -0.44, -0.51, -0.48, -0.56,
             3.55, 1.87, 0.69, 2.40, -0.17, -0.08, 2.03, 1.62, 1.16, -0.59, 1.91, -0.17,
             0.00, -0.46, 1.28, 1.86, -0.12, 1.44, -0.50, 2.21, -0.45, 0.81, -0.72, -0.10,
            -0.53, 0.55, -0.84, 1.95, -0.90, -0.75, -0.81, -0.68, 1.02, 0.63, 0.10, 0.13,
            -0.78, -0.70, 0.34, 0.36, 0.71, -0.86, 1.02, -0.94, -0.26, -0.67, 0.75, 0.48,
            -0.77, 0.51, -0.44, 1.00, -0.81, -0.37, -0.85, -0.72, -0.68, -0.62, -0.60, 0.54,
            -0.94, -0.88, -0.85, -0.77, -0.56, -0.81, -0.66, -0.62, -0.73, -0.73, -0.52, -0.54,
            -0.57, -1.12, -0.97, -0.80, -0.50, -0.98, -0.56, -0.63, -0.77, -0.33, -0.96, -0.83,
            -0.93)

I can use the aov{stats} function to analyze these data as follows:

m.aov = aov(P.conc ~ (N*P + N*K + P*K + rep/block)*depth + Error(plot/depth), data=z)
summary(m.aov)

I can use the lmer{lme4} function to analyze these data as follows:

library(lme4)
m.lmer = lmer(P.conc ~ (N*P + N*K + P*K + rep/block)*depth + (1|plot), data=z)           
anova(m.lmer)

Happily, these analyses with the aov{stats} and lmer{lme4} functions produce identical degrees of freedom, sums of squares, mean squares and F-values.
I then used the mauchly.test{stats} function to evaluate the sphericity assumption as follows:

temp = z[1:32, c("N", "P", "K", "rep", "block")]
mdl.mtrx = model.matrix(~ N*P + N*K + P*K + rep/block, temp)
temp = cbind(z[as.character(z$depth)=="0-5 cm","P.conc"],
             z[as.character(z$depth)=="10-20 cm","P.conc"],
             z[as.character(z$depth)=="50-100 cm","P.conc"])
m.lm = lm(temp ~ mdl.mtrx)
mauchly.test(m.lm, X=~1)

Sadly, the sphericity assumption is rejected. The p-value is ~10^-6 when I retain all five depths.
So, to access Greenhouse-Geisser corrections, I attempted to implement the anova_test{rstatix} function as follows:

library(rstatix)
anova_test(data=z, formula=P.conc~(N*P+N*K+P*K+rep/block)*depth + Error(plot/depth))

This application of the anova_test{rstatix} function produces the following error message:

Error in solve.default(wcrossprod(model.matrix(mod), w = wts)) :
Lapack routine dgesv: system is exactly singular: U[27,27] = 0

So, I next attempted to use the Anova{stats} function as follows:

# m.lm was defined above to implement the mauchly.test{stats} function
ws.factors = unique(z$depth)
m.Anova = Anova(m.lm, idata = data.frame(ws.factors), idesign = ~ws.factors, type=3)

This application of the Anova{stats} function produces the following the following error message:

Error in solve.default(crossprod(model.matrix(mod))) :
Lapack routine dgesv: system is exactly singular: U[2,2] = 0

The two error messages are very similar. This is unsurprising because anova_test{rstatix} is a wrapper around Anova{car} and aov{stats}. Both error messages state “… system is exactly singular …”.

But, the analysis works with the aov{stats} function so we know the system is not singular. The problem must lie with my implementation of the anova_test{rstatix} and Anova{car} functions. What have I done wrong?