Simple prediction using Katla data alone...
model {

for (i in 2:(Nk-1)) {
t[i] <- D[i] - D[i-1]
}
for (i in 2:Nk) {
t[i] ~ dweib(r, mu)I(cens[i],)
}
r ~ dunif(0, 10)
sigma ~ dgamma(0.001, 0.001) # Jeffreys prior
mu <- 1/pow(sigma, r)
median <- sigma*pow(log(2), 1/r) # median time to event
p.erupt.1 <- 1 - exp(pow(92/sigma,r) - pow((92+1)/sigma,r))
p.erupt.5 <- 1 - exp(pow(92/sigma,r) - pow((92+5)/sigma,r))
p.erupt.10 <- 1 - exp(pow(92/sigma,r) - pow((92+10)/sigma,r))
p.erupt.50 <- 1 - exp(pow(92/sigma,r) - pow((92+50)/sigma,r))

}

Data:
list(Nk=19,
D = c(1177, 1262, 1311, 1357, 1416, 1440, 1450, 1500, 1550,
1580, 1612, 1625, 1660, 1721, 1755, 1823, 1860, 1918, NA),
cens = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,92))

Inits:
list(r=1, sigma=10)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   t[19]   106.2   16.22   0.3188   92.33   101.3   148.4   1000   5001
   median   45.23   6.158   0.1093   33.42   45.17   57.87   1000   5001
   p.erupt.1   0.07364   0.03261   5.453E-4   0.02759   0.06752   0.1525   1000   5001
   p.erupt.5   0.3161   0.1157   0.002027   0.132   0.3009   0.5782   1000   5001
   p.erupt.10   0.5286   0.1539   0.002845   0.2502   0.5212   0.8376   1000   5001
   p.erupt.50   0.9647   0.05829   0.00123   0.7927   0.9892   1.0   1000   5001
   r   2.111   0.4266   0.008172   1.376   2.085   3.028   1000   5001
   sigma   54.13   6.865   0.1215   41.9   53.85   68.91   1000   5001
   mu   7.828E-4   0.00176   3.57E-5   4.612E-6   2.49E-4   0.004547   1000   5001

Prediction using Eyjafjallajökull as a time-dependent covariate...
model {
for (i in 2:19) {
t[i] <- D[i] - D[i-1]
inactive[i] <- 1
erupt[i] <- 1
}
# Likelihood contributions from inactive periods
for (i in 2:19) {
inactive[i] ~ dbern(p.inactive[i])
p.inactive[i] <- exp(-mu*pow(t[i], r))
}
# Likelihood contributions at eruption times
for (i in 2:18) {
erupt[i] ~ dbern(q[i])
q[i] <- mu*r*pow(t[i], r-1)*exp(beta*eyja[i])
}
beta ~ dnorm(0, 0.7)
rel.risk <- exp(beta)
r ~ dunif(0, 10)
sigma ~ dgamma(0.001, 0.001)
mu <- 1/pow(sigma, r)
p.erupt.1 <- 1 - exp(mu*rel.risk*
(pow(92, r) - pow(92+1, r)))
p.erupt.50 <- 1 - exp(mu*rel.risk*
(pow(92, r) - pow(92+1, r)))*
exp(mu*
(pow(92+1, r) - pow(92+50, r)))
}

Data:
list(
D = c(1177, 1262 ,1311, 1357 , 1416 , 1440, 1450, 1500, 1550, 1580 , 1612, 1625, 1660, 1721, 1755, 1823, 1860, 1918, 2010),
eyja = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1))

Inits:
list(r=1, sigma=10, beta=0)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   beta   2.862   1.19   0.005721   0.5375   2.857   5.199   4001   216000
   p.erupt.1   0.6546   0.2942   0.00139   0.09711   0.7047   1.0   4001   216000
   p.erupt.50   0.983   0.0386   1.077E-4   0.8677   0.9981   1.0   4001   216000
   r   2.112   0.4135   9.768E-4   1.365   2.09   2.983   4001   216000
   rel.risk   35.6   61.91   0.2786   1.712   17.41   181.1   4001   216000
   sigma   54.1   6.911   0.01656   41.69   53.68   68.95   4001   216000