Here is the model-free explanation: statistical mechanics concerns systems of many (N>>1) interacting elements. Observables can be computed from the partition function $Z$, possibly with added fields. In principle this requires some microscopic probability distribution over dynamical elements, although maximum entropy approaches are also widely used. If we have universality, then not all details of this distribution will be important.
More precisely, observables are obtained by derivatives of $\log Z$. In the thermodynamic limit, the finite observables are of the form
$\mathcal{O} = \lim_{N \to \infty} \frac{1}{N}\frac{\partial \log Z}{\partial h}$
where $h$ is some field-like parameter.
In the simplest class of systems, those we call equilibrium, all but a finite number of degrees of freedom are dynamical: we average over them to obtain expected values of observables. In these systems we compute $Z$, possibly under some approximation, and compute observables as above.
However, in the real world many systems have a large number of degrees of freedom that are not dynamical on the timescales of interest. For example, a spin glass has impurities that don't move; in a glass, the molecules are not free to move arbitrarily, but maintain their neighbours (roughly speaking); in an economy, we may be interested in dynamics over a background of fixed technology and resources. Such systems are called disordered, where the `disorder' is the set of fixed DOF. Call this $J$.
To treat disordered systems, ideally we would like to compute $\mathcal{O}$ at a fixed value of disorder, i.e. $\mathcal{O}[J]$. This is generally impossible because the disorder creates too many parameters to keep track of. So we assign some probability distribution to the disorder $J$, and then we would be happy to compute
$\overline{\mathcal{O}[J]}$
where $\overline{ \cdot }$ denotes the average over the distribution of disorder. Now the problem is that to average $\log Z[J]$, we seem to need to compute $Z[J]$, which is the impossible task. The replica method is based on the trivial identity
$\log z = \left.\frac{\partial z^n}{\partial n}\right|_{n=0}$
which allows us to replace computation of $\log Z[J]$ by $Z[J]^n$. If we can compute this for general $n$, then we can compute $\overline{\mathcal{O}[J]}$ as desired.
Now if $n$ is an integer, then $Z[J]^n$ is the partition function for a system of $n$ replicas of the original system, where the disorder is the same in all the replicas. Generally one considers models where it is easy to average $Z[J]^n$ over the disorder, so we are back to a pure system. But, in this new pure system, the replicas are coupled. These new interactions have the following physical interpretation: in our original system, we are exploring dynamical configurations over a landscape of fixed disorder. In general the disorder restricts the configurations that are possible (or that have a non-negligible probability). (This is what is called `frustration'). Thus, for a given disorder, the different dynamical configurations are correlated. In the replica method, once we have averaged over the disorder, what remains is this inter-replica coupling.
What happens typically is that if frustration is weak, the inter-replica coupling is not important, $\overline{Z^n} \approx \overline{Z}^n$, and the replicas drop out of the problem. In this regime it is equivalent to average $Z$ over the disorder, so physically the disorder may as well have been an extra set of dynamical variables.
Instead if frustration is strong, the inter-replica coupling is important, and the disorder-averaged partition function does not factorize over replicas. In mean-field models, the latter can be reduced to depend only on an $n\times n$ matrix $Q$, called the overlap matrix.
Originally there was a permutation symmetry over the replicas. But, like in any statistical mechanical problem, this symmetry can be broken. Giorgio Parisi understood what is the correct pattern of replica symmetry breaking that applies to the overlap matrix.
It turns out that there is a family of matrices, called hierarchical, in terms of which one can take the $n \to 0$ limit required above. It is important that the matrix size $n$ is under analytical control. Amazingly, in the limit $n \to 0$ there is a new phase of matter, the Parisi-Gardner phase${}^*$, in which the energy landscape is fractal: basins have sub-basins have sub-sub-basins, and so on, to eternity. This is a new phase of matter in which disorder has a very strong effect on the physics. For examples, response functions are generally jerky, leading to avalanches, etc.
Giorgio and collaborators worked out many of the physical consequences of this phase, and applied these constructions to many problems, among them spin glasses, structural glasses, neural networks, random manifolds, and constraint satisfaction problems. Others have applied it to a variety of problems (economies, ecosystems, languages, quantum gravity).
${}^*$ Elisabeth Gardner showed in 1985 that one can have a phase transition from a simple glass to the more exotic Parisi-Gardner phase.