I think I am one step closer to a solution and @jkbest2 has the right idea. For Gamma distributions, I believe Distributions.jl parameterizes the Gamma distribution with shape \alpha and scale \theta while WinBUGS uses shape \alpha and mean parameter \mu = \alpha\theta. So the Gamma(0.001, 0.001) in WinBUGS is Gamma(0.001, 1) or Gamma(0.001) in Julia.
But now I am often getting warnings for numerical errors when I run this model.
@model function corr_model(x, ::Type{T}=Float64) where T
#Prior
r ~ Uniform(-0.999, 0.999)
# mu ~ MvNormal(zeros(2), 1/sqrt(0.001))
mu = Vector{T}(undef, 2)
mu[1] ~ Normal(0, 1/sqrt(0.001))
mu[2] ~ Normal(0, 1/sqrt(0.001))
lambda = Vector{T}(undef, 2)
lambda[1] ~ Gamma(0.001, 1)
lambda[2] ~ Gamma(0.001, 1)
sigma = 1 ./ sqrt.(lambda)
Σ = [1/lambda[1] r*sigma[1]*sigma[2]; r*sigma[1]*sigma[2] 1/lambda[2]]
for i = 1:size(x, 1)
x[i, :] ~ MvNormal(mu, Σ)
end
end