Denoted by $\mu_X$ and $\mu_Y$ respectively. The probability density function of an exponential distribution with mean $\mu$ is given by:

$$f(x) = \frac{1}{\mu} e^{-x/\mu} \quad \text{for } x \geq 0$$

The cumulative distribution function of an exponential distribution with mean $\mu$ is given by:

$$F(x) = 1 - e^{-x/\mu} \quad \text{for } x \geq 0$$

Using these formulas, we can find the mean and variance of X and Y:

$$\text{Mean of X:} \quad E[X] = \mu_X$$

$$\text{Variance of X:} \quad Var[X] = \mu_X^2$$

$$\text{Mean of Y:} \quad E[Y] = \mu_Y$$

$$\text{Variance of Y:} \quad Var[Y] = \mu_Y^2$$

To find the distribution of Z = X + Y, we can use the convolution formula:

$$f_Z(z) = \int_{-\infty}^{\infty} f_X(x) f_Y(z-x) dx$$

Since X and Y are independent, their joint probability density function is given by:

$$f_{X,Y}(x,y) = f_X(x) f_Y(y) = \frac{1}{\mu_X \mu_Y} e^{-x/\mu_X} e^{-y/\mu_Y} \quad \text{for } x, y \geq 0$$

Using this joint probability density function, we can find the distribution of Z as follows:

$$f_Z(z) = \int_{0}^{z} \frac{1}{\mu_X \mu_Y} e^{-x/\mu_X} e^{-(z-x)/\mu_Y} dx$$

$$= \frac{1}{\mu_X \mu_Y} e^{-z/\mu_Y} \int_{0}^{z} e^{x/\mu_X} e^{-(z-x)/\mu_Y} dx$$

$$= \frac{1}{\mu_X \mu_Y} e^{-z/\mu_Y} \int_{0}^{z} e^{-(\mu_Y/\mu_X)x} e^{-z/\mu_Y} dx$$

$$= \frac{1}{\mu_X \mu_Y} e^{-z/\mu_Y} \left(\frac{1-e^{-z(\mu_Y/\mu_X+\mu_Y)}}{\mu_Y/\mu_X+\mu_Y}\right)$$

Simplifying this expression, we get:

$$f_Z(z) = \frac{1}{\mu_X \mu_Y} e^{-z/\mu_X} \left(\frac{1-e^{-z(\mu_X/\mu_Y+\mu_X)}}{\mu_X/\mu_Y+\mu_X}\right)$$

This is the probability density function of Z, which is a Gamma distribution with shape parameter 2 and scale parameter $\mu_X/\mu_Y+\mu_X$.

Therefore, the mean and variance of Z are given by:

$$\text{Mean of Z:} \quad E[Z] = 2(\mu_X + \mu_Y)$$

$$\text{Variance