Creating stacked barplots beside each other from two separate columns

Question

I would like to create a barplot from two seperate column values, meanDM and meanDM1, with meanDM being stacked under meanDM1 while they are beside each other. In other words, I'd like stacked bars from both columns per year. This is a portion of the data (called df):

structure(list(lakeID = c("WE1", "WE1", "WE1", "WE1", "WE1", 
"WE1", "WE1", "WE1", "WE1", "WE1", "WE1", "WE1", "WE2", "WE2", 
"WE2", "WE2", "WE2", "WE2", "WE2", "WE2", "WE2", "WE2", "WE2", 
"WE2"), depth = c("epilimnion", "epilimnion", "epilimnion", "epilimnion", 
"hypolimnion", "hypolimnion", "hypolimnion", "hypolimnion", "metalimnion", 
"metalimnion", "metalimnion", "metalimnion", "epilimnion", "epilimnion", 
"epilimnion", "epilimnion", "hypolimnion", "hypolimnion", "hypolimnion", 
"hypolimnion", "metalimnion", "metalimnion", "metalimnion", "metalimnion"
), year = c(2017, 2018, 2019, 2020, 2017, 2018, 2019, 2020, 2017, 
2018, 2019, 2020, 2017, 2018, 2019, 2020, 2017, 2018, 2019, 2020, 
2017, 2018, 2019, 2020), N = c(16L, 34L, 15L, 35L, 1L, 9L, 5L, 
10L, 16L, 32L, 17L, 33L, 12L, 30L, 22L, 35L, 2L, 9L, 4L, 8L, 
9L, 28L, 22L, 35L), meanDM = c(0.0942303257700237, 0.0661932757564863, 
0.0438792787492709, 0.0668278055831276, 0, 0.036338386383836, 
0.0072085129156294, 0.0107874963737059, 0.0843188069646411, 0.0259432194388038, 
0.0694355219045752, 0.036678461587157, 0.200347848260264, 0.0721742248629981, 
0.144838700987705, 0.41848983842282, 0.00151069595390005, 0.0349446780961297, 
0.130617095377867, 0.017800088456605, 0.00245917385025562, 0.0527044693903904, 
0.0338059254197665, 0.0707387850979652), meanDM1 = c(0.15224547935626, 
0.225104715684883, 0.0897994121303526, 0.177881242805832, 0.0135752775485218, 
0.00373609727796816, 0.0421681308133268, 0.36185063117783, 0.0290467334036447, 
0.051002786726544, 0.0503071643296395, 0.106731388036737, 0.082208813422814, 
0.082824567032327, 0.118888180057524, 0.183962612236507, 0.0229279012476412, 
0.0306285627778002, 0.0337899368763007, 0.00186102094104106, 
0.0965735889574732, 0.047293574778759, 0.0955554969261778, 0.129723908058139
)), row.names = c(NA, -24L), class = "data.frame")

    lakeID     depth   year N    meanDM     meanDM1
1     WE1  epilimnion 2017 16 0.094230326 0.152245479
2     WE1  epilimnion 2018 34 0.066193276 0.225104716
3     WE1  epilimnion 2019 15 0.043879279 0.089799412
4     WE1  epilimnion 2020 35 0.066827806 0.177881243
5     WE1 hypolimnion 2017  1 0.000000000 0.013575278
6     WE1 hypolimnion 2018  9 0.036338386 0.003736097
7     WE1 hypolimnion 2019  5 0.007208513 0.042168131
8     WE1 hypolimnion 2020 10 0.010787496 0.361850631
9     WE1 metalimnion 2017 16 0.084318807 0.029046733
10    WE1 metalimnion 2018 32 0.025943219 0.051002787
11    WE1 metalimnion 2019 17 0.069435522 0.050307164
12    WE1 metalimnion 2020 33 0.036678462 0.106731388
13    WE2  epilimnion 2017 12 0.200347848 0.082208813
14    WE2  epilimnion 2018 30 0.072174225 0.082824567
15    WE2  epilimnion 2019 22 0.144838701 0.118888180
16    WE2  epilimnion 2020 35 0.418489838 0.183962612

I've tried the following code but the barplots havent really showed up.

barCenters <- barplot(height = df$meanDM & df$meanDM1,  beside = TRUE, las = 1, ylim = c(0, 0.8), add=TRUE, cex.names = 1, ylab = "Mean Biomass (mg-AFDM)", xlab = "Year", border = "black", axes = TRUE, legend.text = TRUE, args.legend = list(title = "Year",  x = "topright", cex = .7))

I am using the standard rstudio package and not ggplot2.