Error bars or standard deviation around means from zero-inflated heterogenous data?

I've constructed a plot of mean values per group (year) solely for the purpose of examining my data and have been told to always include either standard error (S.E.) or confidence intervals (C.I.'s). However, some people have informed me that S.E. doesn't make sense for data with a lot of zeros and C.I.'s are difficult to estimate from data with small sample sizes. I hope to better understand why.

I believe I understand the difference between S.E. and C.I.'s as S.E. is an estimation of the population mean from a sample mean, in my case, and C.I.'s are the the quantiles of the distribution the mean's might come from (or +- 1.96x S.E.). Given this, why wouldn't S.E. or C.I.'s work for data with some or all of these characteristics:

1. Small sample size per group (<10) *Wet season has enough data per group, but the other (dry season) does not.
2. Many zeros (75% - 96% per group)
3. Heterogeneity of variance (Does this matter too?)

Can anyone provide a simple example (possibly in R)? Is standard deviation around each mean my only remaining option? I have made a plot with standard deviation and the lower bars go below zero which don't make sense (can't have negative density, in this case).

My data:

>dput(wet_pivot)
structure(list(WYR = c(2010, 2011, 2012, 2013, 2014, 2015, 2016,
2017, 2018, 2019, 2020, 2021, 2022), CYR = c(2009, 2010, 2011,
2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021),
Season = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L), .Label = c("DRY", "WET"), class = "factor"),
N = c(59L, 63L, 69L, 70L, 72L, 71L, 71L, 72L, 71L, 68L, 70L,
48L, 72L), n_mean = c(0.00696806934430411, 0.000649730847004026,
0.00288256551918419, 0.01141088388474, 0.000536174103147671,
0.00349584646220785, 0.000482925207291882, 0.00245359625194744,
0.00292096956686587, 0.00252817293686805, 0.00196286772014134,
0.00501799463867351, 0.00132244297252478), n_median = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), sd = c(0.030946706350869,
0.00248965525641742, 0.0100973832581282, 0.051577934580242,
0.00331468784320076, 0.0266064084754242, 0.00212505905295283,
0.00675243933898364, 0.0119729983336735, 0.00639785127193391,
0.00930625647382774, 0.0136275258272549, 0.00543420856675111
), se = c(0.00402891799826298, 0.000313667078988821, 0.00121558209746373,
0.0061647423020683, 0.000390639708573979, 0.00315759975690469,
0.000252198110662322, 0.000795782607691024, 0.00142093348159893,
0.000775853428563995, 0.00111231039833223, 0.00196696392618855,
0.000640427621321956)), class = "data.frame", row.names = c(NA,
-13L))


Mean values with S.D. around each mean: