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 guassian 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"))
Here are a few links that can help you decide what distribution to use in your analysis:
When to use gamma GLMs?
