Is it appropriate to approximate an unknown probability distribution with a beta distribution? I am creating a model that utilize the output probability distribution of another model, as an input. The model outputs a probability of an event occurring per year, the distribution is a probability between (between 0 and 1) and right skewed, it is defined by Mean = 0.019, Mode = 0.009, 90% confidence interval: 5th percentile = 0.0002, 95th percentile = 0.07. I want to use a beta distribution so I can perform a Bayesian update using a binomial likelihood for the update.
However, I am uncertain of 1) if it is appropriate to use a beta distribution (as opposed to a lognormal) and how to prove if it is appropriate or not, and 2) how to determine the parameters of the beta distribution from the available information.
In relation to 2) I have looked at several articles on cross validate which seem related but I cannot quite work out how to implement:
Calculate the confidence interval for the mean of a beta distribution
Calculating the parameters of a Beta distribution using the mean and variance
This one seems the closest and easiest to implement: Selecting alpha and beta parameters for a beta distribution, based on a mode and a 95% credible interval
In the model I am creating it is possible to define a reparametrized beta distribution according to Mean and standard deviation, so I thought that I would be able to solve equations with the above mentioned values to obtain the standard deviation and define the beta distribution but now I am not sure if this is possible.

A beta distribution or (if a single beta does not do an adequate job) a mixture of beta distributions is quite a sensible choice for quantifying beliefs/knowledge/information about the likely values of a probability (while a log-normal is a less natural choice, because it would consider it possible for the probability to be above 1).
As discussed e.g. by Diaconis and Ylvisaker (1985) a mixture of conjugate priors can approximate any (practically relevant) distribution you want arbitrarily well. The nice thing is that you can even still do conjugate updating with such a mixture: You just do conjugate updating of each of the mixture components and the only tricky bit is how to update the weights of each mixture component. That gets done by multiplying the old weight with the probability mass function of the posterior predictive distribution (for a binomial outcome given the posterior mixture component for the parameter after conjugate updating), which gets you a beta-binomial pmf. Then you normalize the weights. Et voilà, your posterior.
E.g. the RBesT R package provides convenient facilities for doing this:
library(RBesT)
library(ggplot2)
library(patchwork)

# Sample some data from a normal and then use the 
# inverse-logit function to get probabilities.
samples = plogis(rnorm(n=10000, mean=-1, sd=0.5))
hist(samples, breaks=50)

# Find a 3-component mixture to approximate the samples
emfit1 = mixfit(samples, type = "beta", Nc=3)
# For plotting purposes we turn this into a mixture object
# if we wanted to continue the analysis we could just use the emfit object
mixture1 = mixbeta(comp1=emfit1[1:3,1], 
                   comp2=emfit1[1:3,2], 
                   comp3=emfit1[1:3,3], param="ab")
p1 = plot(mixture1, fun=dmix, prob=1)

# Posterior updating with that mixture prior
# and observed data of 10 successes out of 100 tries.
mixture2 = postmix(mixture1, n=100, r=10)
print(mixture2)
p2 = plot(mixture2, fun=dmix, prob=1-1e-13)

(p1 + ggtitle("Prior mixture)) + (p2 + ggtitle("Posterior mixture"))


As you can see, in the example the mixture component that fit the data best got the most weight.
References
Diaconis, P. and Ylvisaker, D. (1985). Quantifying prior opinion. In Bernardo, J. M., DeGroot, M. H., Lindley, D. V., and Smith, A. F. M., editors, Bayesian statistics 2: Proceedings of the second Valencia International Meeting, September 6-10, 1983, pages 133–156, North-Holland, Amsterdam.)
