Error bars and split-plot design

#1
I have the following experimental design :
Sampling is carried out on three random sites (Factor Site is coded : Site1, Site2, Site3).
Within each site, sampling is carried out in a dry plot and in a wet plot (Factor Hydrology is coded : Dry/Wet)
Within each of the previous plot, sampling is carried out at two Depths (D1, D2) in triplicate.
There are a total of n = 36 samples = 3 Sites * 2 Hydrology * 2 Depths * 3 Replicates.

I use the lme function in R (lnme package) as follows :

Code:
Site<-as.factor(rep(c("Site1","Site2","Site3"),each=12))
Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),3))
Depth<-as.factor(rep(rep(c("D1","D2"),each=3),6))
Variable<-rnorm(36)


mydata<-data.frame(Site,Hydrology,Depth,Variable)

mod1<-lme(Variable~ Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth)
I want to visualize the impact of Hydrology and Depth on the dependent variable without displaying the data for each Site :

Code:
interaction.HD<-interaction(mydata$Hydrology,mydata$Depth)
calculated_mean<-tapply(mydata$Variable,interaction.HD,mean)
mean_V=as.table(matrix(calculated_mean,nrow=2))
rownames(mean_V)=c("Wet","Dry")
colnames(mean_V)=c("D1","D2")

library(gplots)
graph=barplot2(mean_V, beside = TRUE)
legend("center",c("Dry","Wet"),pch=22,pt.bg = c("red","yellow"))
Each bar on the graph represents the mean of 9 values (3 replicates per Site for each combination of Hydrology and Depth).

I would like to plot error-bars on each mean. My question is the following : What is the most correct way to calculate it?

I think that it is incorrect to calculate the Standard-Deviation or the Standard-Error on the 9 values used to calculate each mean. I think it would be more correct to calculate SD or SE on 3 values (each value corresponding to the mean per site for each combination of Hydrology and Depth).

Any advice on this question would be much appreciated.