Use the Source, Luke. We can peek inside the ANOVA function by doing getAnywhere(anova.Mermod). The first part of that function is for comparing two different models. The anova on the fixed effects comes in the big else block in the second half:

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
        ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

object is the lmer output. We start to calculate the sum of squares in line 5: ss <- as.vector ... The code multiplies the fixed parameters (in beta) by an upper triangular matrix; then squares each term. Here is that upper triangular matrix for the irrigation example. Each row corresponds to one of the five fixed effects parameters (intercept, 3 degrees of freedom for irrigation, 1 df for variety).

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

The first row gives you the sum of squares for the intercept and the last gives you the SS for the within-fields variety effect. Rows 2-4 only involve the 3 parameters for the irrigation levels, so pre-multiplication gives you three pieces of the SS for irrigation.

These pieces are not interesting in themselves since they come from the default treatment contrast in R, but in line ss <- unlist(lapply(split .... Bates scoops up bits of sums of squares according to the number of levels and which factors they refer to. There's a lot of book-keeping going on here. We get the degrees of freedom as well (which are 3 for irrigation). Then, he gets the mean squares which show up on the printout of anova. Finally, he divides all his means squares by the within-groups residual variance, sigma(object)^2.

So what's going on? The philosophy of lmer has nothing to do with the method of moments approach used by aov. The idea in lmer is to maximize a marginal likelihood obtained by integrating out the unseen random effects. In this case, the random fertility level of each field. Chapter 2 of Pinheiro and Bates describes the ugliness of this process. The RX matrix used to get the sums of squares is their matrix $R_{00}$ from equation 2.17, page 70 of the text. This matrix is obtained from a bunch of QR decompositions upon (amongst other things) the design matrix of the random effects and $\sqrt{\sigma^2/\sigma_f^2}$, where $\sigma_f^2$ is the variance of the field effect. This is that missing factor you were asking about, but it does not enter into the solution in a transparent or simple way.

Asymptotically, the estimates of the fixed effects have distribution:

$$ \hat{\beta} \sim \mathcal{N}(\beta, \sigma^2[R_{00}^{-1}R_{00}^{-T}]) $$

which means that the components of $R_{00}\hat{\beta}$ are independently distributed. If $\beta=0$, the squares of those terms (divided by $\sigma^2$) follow a chi-squared distribution. Dividing through by the estimate of $\sigma^2$ (another chi-squared when divided by $\sigma^2$) gives an F statistic. We don't divide by the mean square error for the between-groups analysis because that has nothing to do with what's going on here. We need to compare the sums of squares obtained through $R_{00}$ by comparing them to an estimate of $\sigma^2$.

Note that you would not have got the same F statistics had the data been unbalanced. Nor would you have obtained the same F statistics if you had used ML instead of REML.

The idea behind aov is that the expected mean squares for irrigation is a function of $\sigma^2$, $\sigma_f^2$ and the irrigation effects. The expected mean square for field residuals is a function of $\sigma^2$ and $\sigma_f^2$. When the irrigation effects are 0, both of these quantities estimate the same thing and their ratio follows an F distribution.

Interestingly, Bates and Pinheiro recommend using the ANOVA over fitting two models and doing a likelihood ratio test. The latter tends to be anti-conservative.

Of course, if the data are unbalanced, it is no longer clear exactly what hypothesis the ANOVA is testing. I removed the first observation from the irrigation data and refit. Here is the $R_{00}$ matrix again:

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

As you can see, the sums of squares for the irrigation parameters now contain some of the variety effect as well.