**Answer Q1:** To know if the model is correct we need more information. It depends on how you randomized the `TreatN` and `TreatH` combinations in each site (i.e. `FactorS`). If you assigned the `TreatN` and `TreatH` combinations as a completely randomized design then I would say that yes your model is correct. If you randomized the the treatment combinations as a randomized complete block design I think the model should be:   

    m2 = lmer(area ~ treatN*treatH + (1|FactorS/replicate), 
              data = data) 

because your replicate/block is nested within location.

**Answer Q2:** `lmer` fits mixed-effect models and is a type of generalized linear mixed model with a Gaussian distribution. The function `glm()` can't fit random effects.

**Answer Q3:** You can't fit the model you specified using the `glm()` function unless you treat `FactorS` as a fixed-effect; you could use `glmer()` function doing the following:

    m4 = glmer(area ~ treatN*treatH + (1|FactorS), 
               data = data, family = Gamma(link = "identity"))   

If your data follows a normal distribution you can also use the Gaussian distribtion in `glmer`, which is the same as doing the analysis in `lmer`.       

    m4 = glmer(area ~ treatN*treatH + (1|FactorS), 
               data = data, family = gaussian(link = "identity"))

(this will give you a warning saying you should just use `lmer` instead).

Here are a few links that can help you decide what distribution to use in your analysis:  
https://stats.stackexchange.com/questions/67547/when-to-use-gamma-glms