Slope estimator robust to repeated values I have a hydrologic dataset that contains many repeated values (in my case, 0s), and I want to understand whether there has been a trend through time. Here is an idealized example of what the data look like:
# make sample data with linear increase
year <- seq(1, 51)
value <- seq(0, 500, 10)

# replace just over 50% of values with 0s
value[seq(1,51,2)] <- 0

plot(year, value)


Typically, I would use the non-parametric Mann-Kendall test to determine whether there is a significant change through time, and the Theil-Sen slope estimator to determine the slope of that change.
In this case, the Mann-Kendall test rejects the null hypothesis of no change through time and has a positive tau, interpreted as a statistically-significant increase through time. However, the Theil-Sen estimator returns a value of 0:
manken <- rkt::rkt(year, value)
manken$tau  # Kendall tau = 0.235
manken$sl   # p-value = 0.009
manken$B    # Theil-Sen estimator = 0

I believe this is happening because the Theil-Sen estimator returns the median slope for all pairs. Therefore, whenever there are >50% of points that have the same value, the median slope (I believe) will always be equal to 0.
Question: Are there statistical slope estimators that are better suited for data that have many repeated values?
For what it's worth, a linear model (lm(value ~ year)) also returns a significant (p=0.001) positive slope. And so does the highly scientific "eyeball test".
Edit 8/3/2020: For additional context, my real y-axis is the number of days with zero flow per year. So, this issue arises when there is a stream that flows year-round for >50% of years, but has goes dry for parts or all of some years. I'd like to know whether "dryness" (frequency/duration) is increasing.
 A: 
For additional context, my real y-axis is the number of days with zero
flow per year.

Your simulated data does not contain any uncertainty and is therefore not very useful. Also, how can a year have more than 365 days? I will simulate your dependent as a count variable, i.e., with a Poisson distribution.
I'm not an expert for count models and I haven't seen your actual data, so other distributions (such as the negative binomial) might be better for your model. It might even be necessary to use a distribution with an upper limit (if you have values close to 365 days).
I will use a binomial distribution to simulate your zero values. Again, other distributions might represent your data better.
set.seed(42)

year <- seq(1, 51)
value <- rpois(length(year), lambda = exp(year * 0.07 + 1))

# replace about 50% of values with 0s
value[as.logical(rbinom(length(value), 1, 0.5))] <- 0
mean(value == 0)
#[1] 0.5686275

DF <- data.frame(year, value)

plot(value ~ year, data = DF)

We can now fit a hurdle model. A hurdle model combines two models. The first one models if values are zero or non-zero. The second one models the non-zero values. Both of these are generalized linear models.
library(pscl)
fit <- hurdle(value ~ year, dist = "poisson", zero.dist = "binomial", data = DF)

summary(fit)
#Call:
#hurdle(formula = value ~ year, dist = "poisson", zero.dist = "binomial")
#
#Pearson residuals:
#    Min      1Q  Median      3Q     Max 
#-0.7026 -0.6698 -0.6171  1.1072  2.0128 
#
#Count model coefficients (truncated poisson with log link):
#            Estimate Std. Error z value Pr(>|z|)    
#(Intercept) 1.458876   0.151480   9.631   <2e-16 ***
#year        0.058816   0.003802  15.471   <2e-16 ***
#Zero hurdle model coefficients (binomial with logit link):
#             Estimate Std. Error z value Pr(>|z|)
#(Intercept) -0.881742   0.618411  -1.426    0.154
#year         0.003785   0.020518   0.184    0.854
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#
#Number of iterations in BFGS optimization: 11 
#Log-likelihood: -73.36 on 4 Df

As you see, the output tells us that the probability of a value being zero is independent of the year (that's how we simulated it). The count model shows a strongly significant intercept and slope (note the log link). Let's plot predictions from the count model:
curve(predict(fit, type = "count", newdata = data.frame(year = x)), add = TRUE, col = "red")


I believe hurdle models could help you but you'd need to investigate a bit more which assumptions would be sensible regarding the distributions and link functions. Of course, for this it would be helpful to have mechanistic knowledge about why non-zero values occur and what might cause the increase with time. Additional predictors would be useful.

I'd like to know whether "dryness" (frequency/duration) is increasing.

The zero model would tell you if frequency of years with zero flow depends on time (in the simulated data it does not). The count model would tell you if the number of days with zero flow in dry years ("severity" of dryness) depends on time.
Note that hurdle models are for zero-inflated data, they assume that two "processes" are involved. One controls if a value is non-zero, the other the magnitude of non-zero values. Your simulated data supports this assumption. Your real data might not be zero-inflated.
A: The Theil-Sen estimator is robust against a skewed or fat-tailed error distribution, but it still assumes that the model is still linear and has a single, well-defined slope. That is, it assumes that the specification is still of the form
$$ y = \beta_1 x + \beta_0 + \epsilon $$
But instead of assuming $\epsilon \sim \mathcal{N}(0, \sigma^2)$ it places no restriction on the distribution of $\epsilon$ except perhaps mean 0. According to Wikipedia, it is robust to up to 29% of points being changed.
However, the way you generated the data (called a mixture model) doesn't assume there is a single slope, but rather posits that there are two classes, each with a different slope, that have been randomly mixed together.
The Kendall $\rm{T}$ test still works, because it is clear the mixture is not a random order, but the slope estimation isn't applicable because the assumptions have been badly violated.
There is a standard way to fit models to data generated by mixture models: latent variable models and the Expectation-Maximization algorithm. You data could be called a "mixture of regressions." There is an R package which can handle this case called flexmix. Here is how I would use that package to fit your fake data. The k=2 parameter is telling it there are two classes, which we know a priori.
# generate 51 equally spaced points along a line
year <- seq(1, 51)
value <- seq(0, 500, 10)

# add a little bit of noise to prevent likelihood underflow
value <- value + rnorm(n=51, mean=0, sd=1)

# replace just over 50% of values with 0s
value[seq(1,51,2)] <- 0


#install.packages("flexmix")
library(flexmix)

1model <- flexmix(value ~ year, k=2)
summary(model)

plot(year, value, col = clusters(model), pch=19)
abline(parameters(model)[1:2, 1], col = "black", lty=2)
abline(parameters(model)[1:2, 2], col = "red", lty=2)


Call:
flexmix(formula = value ~ year, k = 2)

       prior size post>0 ratio
Comp.1   0.5   26     26 1.000
Comp.2   0.5   25     26 0.962

'log Lik.' -118.8863 (df=7)
AIC: 251.7726   BIC: 265.2954

The way the EM algorithm works is by guessing which class each point belongs to. It starts by assuming that each point has a 50% chance of being in each class. Then it fits a weighted regression model for each class. Then, based on the two fitted regression models, it goes back and updates the probabilities of being in each class for every point. For example, if a point was initially assumed to be equally likely to be in either class but ended up very close to the regression line for class 1 and very far from the regression line for class 2 after the first iteration, its probabilities would be updated to 80% for class 1 and 20% for class 2. This process then repeats until convergence is reached. At that point, we have a pretty good guess of which class each point came from, and two separate regression lines; because of weighting, we can imagine that each line was fit only to those points which are likely to belong to the same class.
The EM algorithm is good but not perfect. The hyperparameter k must be chosen very carefully. Although likelihood is guaranteed to increase with each iteration, the algorithm can sometimes be unstable and converge to different solutions if fit to a different random subsample of the data. In some cases, the likelihood can actually go off to infinity; this actually happens with your fake data set because all of the data lie in a perfectly straight line! (Adding a little bit of random noise fixes that problem, which is very unlikely to occur in read world data anyway.) However, if the assumptions are met it can be a very powerful technique.
