Exploring the limits of ensemble forecasting via solutions of the Liouville equation for realistic geophysical models

The atmosphere is an exemplary case of uncertain system. The state of such systems is described by means of probability density functions which encompass uncertainty information. In this regard, the Liouville equation is the theoretical framework to predict the evolution of the state of uncertain systems. This study analyses the morphological characteristics of the time evolution of probability density functions for some low complexity geophysical systems by solving the Liouville equation in order to obtain tractable solutions which are otherwise unfeasible with currently available computational resources. The current and usual modest approach to over-come these obstacles and estimate the probability density function of the system in realistic weather and climate applications is the use of a discrete and small number of samples of the state of the system, evolved individually in a deterministic, perhaps sometimes stochastic, way. We investigate particular solutions of the shallow water equations and the barotropic model that allow to apply the Liouville formalism to explore its topological characteristics and interpret them in terms of the ensemble prediction system approach. We provide quantitative evidences of the high variability that solutions to Liouville equation may present, challenging currently accepted uses and interpretations of ensemble forecasts.


Introduction
Understanding the fundamentals of weather forecasting is one of the most challenging problems the scientific community currently faces, not only for its academic value but also for its potential benefits on multiple socioeconomic assets. Among other less critical applications, weather forecasts are routinely used to guide civil protection agencies in their mission against severe weather events in order to prevent personal and material losses (WMO, 2012;Yano et al., 2018).
From the seminal works of Charney in the early 1950s (Charney, 1951), we recognize that unlike other natural phenomena such as tides, which can be precisely predicted from accurate observations, the atmosphere does not manifest periodicities that enable to predict its state from recurrence. The common alternative forecasting approach consists in solving the governing equations of the dynamical system, initialized from a representation of the actual conditions of the day. However, multiple error sources degrade the forecasting process.
Current numerical weather prediction models have essentially two sources of errors: the description of the initial state of the atmosphere, and the modeling of the physical processes and their interaction with external components (oceans, orography and cryosphere among others). The lack of infinitely precise knowledge about the state of the atmosphere, partly due to an insufficient observational sampling compared to the enormous number of degrees of freedom of the system, induce errors in the description of the system. Despite the advances in variational data assimilation methods currently used to produce the best estimate of the state of the system (Carrassi et al., 2018), initial conditions for numerical forecasts still contain relevant errors, which need to be characterized and adequately sampled.
The second source of uncertainty arises from the incomplete knowledge of some physical processes and their discrete implementation in numerical models. The parametrization of sub-grid scale mechanisms required to account for not fully represented or even absent processes is one of the main causes of model error (Palmer, 2000). There are different approaches to account for model uncertainties, such as multi-model ensembles (Harrison et al., 1999) or stochastic parameterizations, (e.g. Buizza et al., 1999;Berner et al., 2017).
Forecasting the atmosphere involves uncertainty, and the mathematical framework to model uncertainties is probability theory. The degree of uncertainty in the state of a system is expressed by means of a probability density function (PDF). Such multidimensional function is defined over a vector space that represents all possible states of the https://doi.org/10.1016/j.atmosres.2020.105127 Received 19 February 2020; Received in revised form 16 June 2020; Accepted 3 July 2020 system, i.e. the phase space (Lorenz, 1963). While a fully deterministic evolution of the system is represented by a single trajectory in phase space, the perfect-model evolution of an uncertain state (i.e., the state PDF) fulfills the Liouville equation (LE; Ehrendorfer, 2006). A more general framework, which includes stochastic sources of model error is the Fokker-Planck equation (Hasselmann, 1976;Thompson, 1985).
Although the importance of initial conditions was already pointed out by Thompson (1957), the work of Lorenz (1963) was the first that specifically reported that two slightly different initial conditions may evolve into highly different states, eventually rendering a worthless forecast. The notion of deterministic chaos referring to irregular, unpredictable behavior in deterministic, dissipative, and nonlinear dynamical systems such as the atmosphere, as well as the existence of strange attractors, which display high sensitivity to initial conditions is currently well known (Zeng et al., 1993). However, the implications of this chaotic behavior on the topological characteristics of the PDF still remain open. Indeed, even in systems without the aperiodic chaotic behavior, the characteristics of the corresponding PDF may be complex. These aspects are clearly scale dependent. At very small scales, turbulence is the most critical factor in limiting the atmospheric predictability, specially in the convective scale, whereas chaos has only a secondary role. In particular, in the fully turbulent regime of convective scales, the intrinsic characteristics of the flow limit the predictability more critically than errors in the initial conditions, which can be ultimately reduced (Yano et al., 2018).
The true dynamics of the atmosphere is described by the Navier-Stokes equations, which are approximated with ordinary differential equations in a finite dimension space. The high dimensionality of this space severely restrict the applicability of the LE for operational weather and climate applications. In a pioneering work exploring the solutions of the LE, Epstein (1969) computed the first and second moments of the PDF by integrating the stochastic dynamic equations (SDE). These equations describe the time evolution of the PDF moments, and can be derived from the prognostic equations of the dynamical system. However, the artificial requirement for a closure assumption and its computational cost are two intrinsic limitations of this method. Further work on SDE was done by other authors, such as Pitcher (1977), who integrated these equations for a barotropic model using realistic data. Additionally, Pannekoucke et al. (2018) and Pannekoucke and Fablet (2020) offer a recent contribution to the characterization of the PDF moments under nonlinear dynamics in the infinite dimensional functional space by means of a parametric approach in which the error covariance matrix is approximated by evolving parameters in time and space.
On the other hand, Leith (1974) introduced ensemble forecasting with a Monte Carlo technique, which consists in randomly sampling the initial PDF and evolving deterministically each sample with the dynamical model. For real weather and climate applications, ensemble forecasting is the only feasible strategy to estimate the PDF of the atmospheric system as the solution of the LE for a such high dimensional system is unapproachable (Molteni et al., 1996;Toth and Kalnay, 1997). Since 1992 ensemble prediction systems have been produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) and the National Centers for Environmental Prediction (NCEP) among other operational centers. They apply methods based on the identification of fast growing modes (Buizza and Palmer, 1995;Toth and Kalnay, 1993) and sampling model error usually by stochastic techniques (Leutbecher et al., 2017).
In the weather and climate ensemble prediction framework, most methods used to derive PDFs, explicitly or implicitly, assume the underlying unknown forecast PDF has a locally connected support (i.e. the set of points where the function is not 0), is smooth, correlated among neighboring states and presents weak gradients of probability. Nevertheless, the real topological characteristics of the PDF are unknown, so this general assumptions may not be adequate. For instance, the actual PDF may present strong gradients of probability with high probability states surrounded by a low probability region in phase space (a schematic example of this topology was described by Penrose, 1989, andalso by Ehrendorfer, 2006). In that case, computing probabilities under the a connected PDF support assumption would render misleading probabilistic forecasts for either high and low probability states.
Regarding the full solution of the LE, Ehrendorfer (1994aEhrendorfer ( , 1994b found solutions for some low-dimensional systems. In addition, Ehrendorfer (2006) introduces solutions of the LE which display nonlinear features as time progresses. In more recent works, Yano and Ouchtar (2017) solved the LE for a simple and low dimensional convective system and Garret (2019) applies the LE to obtain the size distribution of precipitating particles.
This study is framed in the context of geophysical systems in which a numerical model containing the dynamical equations and the physics parametrization is used and only ensemble forecasting techniques can be used to estimate the PDF of the system due to the extremely high number of degrees of freedom of the system. The purpose of this work is to explore the topological characteristics of Liouville solutions and report on topological features present in low complexity geophysical systems which challenge current uses and interpretations of ensemble prediction systems. Special emphasis is put onto solutions of the LE presenting complex support topologies due to its implications for ensemble forecasting. This study extends the analysis of Ehrendorfer (1994bEhrendorfer ( , 2006, exploring the implications for ensemble prediction systems of solutions of the LE obtained for simplified models.
The rest of the paper is organized as follows. Section 2 presents the general solution of the LE, in Sections 3 and 4 an analytical solution of the LE is obtained for simplified cases of the shallow water and barotropic models. The main conclusions are summarized in Section 5.

General solution of the Liouville equation
Given a dynamical system with N degrees of freedom = X Ẋ Φ( ), where X is the system vector state, the Liouville equation for the probability density function ρ has the following expression: is the divergence of the trajectories of the system in phase space.
LE is a continuity equation for probability in phase space, equivalent to mass continuity equation in physical space. The PDF ρ must also satisfy the following conditions: The Liouville solutions remain normalized as long as the initial ρ 0 is normalized (e.g. Ehrendorfer, 1994a).
A general solution of the LE can be found by using the method of characteristics (Zwillinger, 1989). See Appendix A for details of the derivation of the following general solution for the LE: where α is the state vector at t = 0 that takes the system to state X at time t. ρ 0 is the initial PDF evaluated at state α. In order to evaluate ρ 0 , α must be expressed in terms of X (i.e. α(X, t)). In other words, the trajectory connecting states α and X must be computed. The trajectory can be calculated analytically when both the analytical solution of the dynamical system = X Ẋ Φ( ) (i.e. X(α, t)) can be derived, and this solution can be inverted. To put it another way, the computation of ρ 0 (α) involves the computation of backward trajectories which connect state X at time t with the corresponding state α at t = 0.
The integral in Eq.
(3) depends on the divergence in phase space χ. When no divergence occurs in phase space, the integral equals 0 and the solution for ρ reduces to the first term. However, in general the integral of the second term must usually be computed numerically, although for some simple systems an analytical expression exists.
Therefore, the possibility to solve the LE analytically, even having a general expression for ρ is restricted to simplified dynamical systems Fig. 1. Cross-sections of the PDF ρ in the U 1 -H plane for the value of A with maximum probability at: a) t = 0, b) t = π/12, c) t = π/6 d) t = π/4, e) t = 3π/4, f) t = π, g) t = 5π/4 and h) t = 7π/4. Note the different scales for each time step. A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 which have analytical solution. Nevertheless, Eq.
(3) provides an additional method to solve the LE numerically by computing a set of forward trajectories in phase space. In particular, the integration of the dynamical system connects α and X. According to the general solution, ρ at time t equals ρ 0 (α) at time 0 in a non-divergent system in phase space.

Solution of the Liouville equation for a simplified shallow water equations system
The general solution of the LE is applied to a special case of the shallow water equations with the effect of the Earth rotation over an arbitrary domain: where u = ∂ t x and v = ∂ t y are the x, y velocity components, h is the depth, f is the Coriolis parameter, g is gravity and D/ Dt = ∂ t + u∂ x + v∂ y . A solution quadratic in x and y for h and linear for u and v satisfies the system (4). The analytical solution of the LE for this system after some simplifications is the following: where ρ 0 is the probability density function at t = 0, which is assumed to be known. The second term of the solution corresponds to the integral of the divergence in phase space χ, which can be computed analytically for this simple system. In order to obtain an explicit solution for ρ at time t and any phase space point U 1 , A, H, the solution of the system at time t depending on the initial conditions U 1 0 , A 0 , H 0 must be inverted (see Appendix B). In this case, this procedure can be done analytically, which yields an expression for the initial condition that takes the system to state (U 1 , A, H) at time t: cos ( ) (cos( ) 2 sin( )) , cos( ) cos( ) 2 sin( ) .
As an illustrative example, we take a three dimensional Gaussian distribution as initial PDF. The distribution is centered at point (μ U 1 = 0.25, μ A = − 0.01, μ H = 1) and standard deviations are:

Results
The evolution of the system, solved analytically, is analyzed by representing the probability at each gridbox, calculated from ρ(t), and the corresponding ΔV in phase space (ΔV = ΔU 1 ΔAΔH). At t = 0, ρ is a three dimensional Gaussian distribution centered at (μ U 1 , μ A , μ H ) (Fig. 1a). As time passes, the distribution moves through phase space maintaining the initial connected nature of the PDF and changing its shape and support due to the action of phase space divergence (Fig. 1). Since U 1 , A and H are periodic and the periodicity is the same for all states in the initial condition, ρ is also periodic, so ρ(2πn) = ρ 0 for any integer n. With regard to the predictability of this simple dynamical system, besides its periodicity, the high probable states remain connected at all times, its predictability will be large, and conventional discrete ensemble forecasts would be adequate to cope with its uncertainty. Additionally, improved predictability occurs at the beginning of the evolution due to a convergence phase. Indeed, uncertainty is not always a growing function in time. This feature was also pointed out by Yano and Ouchtar (2017).
The evolution of this system in physical space shows an oscillation between a wide shallow vortex and a small deep vortex. The maximum depth of the vortex is represented by H 0 and the radius R by −H A / . The oscillation in the expected value of the depth and radius, as well as their uncertainty is shown in Fig. 2. The uncertainty of coefficients A (not shown) and H is larger in the last part of the period (Fig. 2a), which is consistent with the evolution of ρ (Fig. 1). Larger uncertainty in the state variables is associated with divergent phase space trajectories. However, the uncertainty of quantities such as R derived from state variables (Fig. 2b) does not satisfy the same relation with divergence in phase space as it is greater during the first part of the period, when trajectories in phase space are convergent.

Solution of the equation
An additional system which admits an analytical solution of the LE and has some degree of geophysical realism is a special solution of the barotropic vorticity equation: where ψ is the streamfunction and β = df/dy with f the Coriolis parameter.
We consider an extension of the solution proposed by Lorenz (1960), introducing a linear variation of the Coriolis parameter with latitude (e.g. a β-plane solution), also used by Paegle and Robl (1977), which consists in a highly truncated Fourier series expansion for the Laplacian of the streamfunction with three time dependent amplitudes and phase angles that depend linearly on time: Substituting (9) into (8) and equating the terms with the same trigonometric function yields: The equations in system (10) are dimensionless and the time unit is arbitrary as the topology of the Liouville solution is not sensitive to these dimensions. The evolution in phase space under the dynamics of this system can be inspected before solving the differential equations. The three equations of the system can be expressed as Taking the first equality yields which is integrable. With an analogous procedure for the other two equalities and integrating we obtain the following solutions: where K 12 , K 13 and K 23 are constants of integration. Absolute value is used to make the sign of the coefficients C 1 , C 2 , C 3 explicit. In the X 1 -X 2 and X 2 -X 3 planes the orbits are ellipses, whereas in the X 1 -X 3 plane the curve is a hyperbola (Fig. 3). The result of this development is obtained assuming that X 1 , X 2 and X 3 are not 0. When X i , i=1, 2 or 3 is 0, the trajectory changes its direction in a plane perpendicular to X i . This can be checked in the dynamical system, when X i changes its sign, Xj, j ≠ i also does and thus the direction of the trajectory is reversed.
The dynamical system (10) has analytical solution given by: , .
The solutions of X 1 , X 2 and X 3 are periodic with periods 2κ/h for X 1 and 4κ/h for X 2 and X 3 , with Given that k 0 2 depends on the initial conditions, the value of κ and thus the period of the solutions of the dynamical system (10) has a dependence on the initial conditions as well. Since Xk, k = 1,2,3 are independent of X k , the χ term in Eq. (1) is zero, and thus the non-divergent solution of the LE is obtained (i.e. ρ(X, t) = ρ(α)): = ρ X X X t ρ X X X X t X X X X t X X X X t ( , , , ) ( ( , , , ), ( , , , ), ( , , , )).
An explicit analytical solution for the LE requires to have an explicit A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 expression for X k0 , k = 1,2,3, which is the initial state (α in Eq. (3)) that takes the system to state X at time t. Due to the time invariance of system (10), the initial condition can be expressed in terms of the state at time t by using the functions in (13) and changing X k0 by X k and t byt (Ehrendorfer, 1994b). Note that, in strict sense, this solution is not purely analytical as the Jacobi elliptic functions must are computed numerically by using the Python SciPy package.
As for the previous case, a three dimensional Gaussian distribution is taken as initial ρ 0 : Fig. 3. Projection of individual trajectories on a), b) X 1 , X 2 , c), d) X 1 , X 3 , e) f) X 2 , X 3 planes for two different initial conditions: X 1 = 0.12, X 2 = 0.24, X 3 = 0.1 (a, c, e) and X 1 = 0.5, X 2 = 0.7, X 3 = 0.7 (b, d, f) with δ = 2. A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 where K is a normalization constant, μ k are the expected values in the X k directions and σ is the standard deviation, taken equal for all three state variables.
A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 spanning in the 3 dimensions, surrounded by low probability states. This generates increasingly strong gradients of probability in smaller and smaller regions. Despite the fact that each X k is periodic, ρ(X, t) is not because periods are different for each initial condition as they depend on parameters h and k 0 2 (Eq. (14)), which in turn depend on X 0 . Given a point in phase space, X p , ρ(X p ) is periodic with period 4κ(X p )/h(X p ). A different point X q has a different period, so the full probability density function is not periodic. After some time, trajectories initiated from neighboring points evolve into distant points. Conversely, initially distant states can evolve towards the same phase space region. Nevertheless, the distribution of periods in phase space (Fig. 5) is continuous, so that the continuity between high probability states is not broken with time. In other words, high probability states which are initially connected remain connected as the system evolves, but the distance between these states increases with time, while the distance between some high and low probability states can decrease and become neighbors at specific times (Fig. 6). A segment of high probability states stretches during his evolution under system (10) increasing the distance between these high probability states. Simultaneously, a segment of low probable states evolves analogously, so that filaments of high and low probability are mixed (Fig. 7). This stretching effect produced by this linear system resembles the stretching and folding of chaotic systems. The folding process does not occur in this system as the periodicity of the solutions already bounds the support to a limited region. The X 1 direction presents the highest rate of change of the period (Fig. 5b) and also the highest gradients of probability (Fig. 4e,f,h,i,k,l), which enhances the previously described effect for small variations of the initial condition in this direction. This proximity between high and low probability states generated by the evolution of dynamical system produces strong gradients in the resulting PDF, rendering the highly variable distribution observed in Fig. 4. The filamentary topology of the solution is better illustrated considering a plane perpendicular to trajectories such as the plane X 2 =0 (Fig. 8). For short times, only a thick filament crosses the plane. As time evolves, the filament stretches and due to the different periodicities of each state, regions of appreciable probability density have low probability states in between. Bear in mind that all initial states which belong to a curve with high probability density states on its extremes remain connected under the system evolution, but the distance between extremes increases due to the stretching effect, which renders on the one hand, a reduction of the curve section and on the other hand, the possibility that the same curve crosses a perpendicular plane more than once, as the set of possible are bounded by the system dynamics and thus, the longitude of the curve can only be increased if the filament makes various revolutions. The presence of these high gradients of probability reveals a serious predictability challenge for this system: usual probability forecasts derived from ensemble prediction systems would artificially assign large values of probability density to marginally probable states. As an illustrative example of this situation applied to the previously described system, we consider the probability p of specific events, characterized by a volume Ω in phase space. For different times, the volume in which p is computed is chosen so that it is centered around the maximum value of ρ. The probability p is computed as follows: The probability is also computed for a smoothed PDF, computed applying a Gaussian smoothing filter with σ = 2 to assess the impact of neglecting the filamentary topology of the actual PDF. The ratio between the probability computed from the smoothed PDF, p sm and from the real PDF, p as a function of the volume considered for different times shows underestimations of probabilities due to an incorrect assignment of lower probabilities to states in Ω, which must be compensated by higher probabilities outside the region considered (Fig. 9). For larger volumes, p sm converges to p, due to the filtering out effect of the small scale variability involved in the averaging of Eq. (16). The exact same effect is observed for large times: given a volume size in phase space and a tolerable accuracy error for the probabilistic forecast, there exist a time after which the stretching of the solution renders p sm acceptable. Before that time, p sm is strongly affected by the scale of the PDF variability. Therefore, the error in the estimation of probability depends on the relation of the scale of the variability of the PDF (i.e. the volume of the high and low probability filaments) and the user's tolerance range of the event whose probability is computed (i.e. the volume Ω considered). If the scale of the phase space volume is of the same order than the scale of the filaments, probability computations are strongly affected. Conversely, if the scale of interest is larger than the scale of the variability of the PDF and the distribution of high and low probable states is uniform, the presence of probability gradients in the PDF has no influence on the resulting probability forecasts. In this system, the scale of the filaments reduces with time, rendering more accurate probabilities for small volumes and large times. The translation of the previous example into a real meteorological application would imply that a specific temperature range (e.g. T = [10,11]°C), could include negligible probability states in the forecast. However, if  10) in a) X 2 -X 3 plane, X 1 =0.12 and b) X 1 -X 2 plane, X 3 =0.1. Dashed contour encloses the region where initial probability densities are higher than 1% of the maximum probability density. A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 two ensemble members predicted values of e.g. 10 and 11°C, typical postprocessing tools would assign high probabilities to all states in the range T ∈ [10, 11]°C. Hence, ignoring the potential highly variable character of the actual PDF would render incorrect probabilistic forecasts of T ∈ [10, 11]°C, assigning artificially high probabilities to this range. Eventually, in a filamentary PDF topology with filament scale decreasing with time, this effect would not significantly affect the calculation of probability of T ∈ [10, 11]°C as underestimations and overestimations of probability density could be compensated as the scale of high and low probability density filaments decrease with time. An alternative initial condition for ρ (i.e. ρ 2 ) with μ 1 = 0.5, μ 2 = 0.7, μ 3 = 0.7 and the same value of σ renders a similar evolution for ρ 2 (Fig. 10) to that depicted in Fig. 4 for ρ 1 . The initial probability density function is formed by spherical shells of constant probability density (Fig. 10a) that evolve into an ellipse in the X 2 -X 3 plane centered close to the point (μ 1 ,0,0) displaying high gradients of probability (Fig. 10d). Nevertheless, with this initial condition the spread along X 2 and X 3 is larger and the blend of high and low probable states arises at earlier times so that it can be seen even at t=50 (Fig. 10b). The difference between ρ 2 and ρ 1 in terms of Eq. (12) is the value of the constants K 12 , K 13 and K 23 , which are associated to each individual trajectory in phase space.
These solutions depend on an additional parameter δ, which is the ratio between the wavenumber in the two spatial directions (Eq. (11)).
To study its effect on the solution for the probability density function, the value of δ is changed to 50 (C 1 = −8⋅10 −6 , C 2 = 49.98, C 3 = −24.99). For this value of δ (Fig. 11), in the X 2 -X 3 plane the effect of the different periods of the solution, which produces a highly variable PDF (Fig. 4), is reduced but it is increased in the X 1 direction producing extremely high gradients of the PDF in this direction (Fig. 11). Hence, initially neighboring states in the X 1 direction separate much faster than for δ = 2. In this case, ≈ Ẋ0 1 , so that state evolution is mainly twodimensional and different for each X 2 − X 3 plane. The change in δ modifies constants C 1 , C 2 and C 3 (Eq. (11)), rendering an evolution with X 1 approximately constant and the elliptical trajectory in the X 2 -X 3 plane.

Impact of the beta effect
When β is not zero, the argument of the Jacobi elliptic functions in (13) is periodic in time. Unlike the solution with no β effect, not only X k , but also ρ is periodic due to the sinus term, not present in the β = 0 solutions. A solution with high gradients of probability is still obtained for high values in the argument of the Jacobi elliptic functions (Eq. (13)). The β effect produces periodicity in the argument of the Jacobi elliptic functions, conditioning the high variability of the PDF to the period of τ (Eq. (13)), which is given by 2π/γ. For small values of β of O (10 −11 ), the impact of the β effect on the probability density function is Fig. 6. States with probabilities > 80% of probability maximum (blue) and states with probability < 10 −3 of the probability maximum (red) in a region surrounding the probability maximum at a) t = 100, b) t = 500 and c) t = 5000. Projections along the three two-dimensional plains are provided as well.
noted for times of the order of the period of τ. For short times compared to the period of τ, a PDF with high gradients of probability is obtained. However, because of the periodicity of the argument of the Jacobi elliptic functions, at t = nπ/γ, the sinus term vanishes and ρ returns to its initial form ρ 0 (not shown).
For higher values of β, the period of ρ is shorter, and for large enough values, such as β = 2.5 ⋅ 10 −7 , the spiral obtained for the β = 0 solution does not form and ρ is rather smooth for all times (Fig. 12). High variability is not produced in this case because the period of ρ is shorter than the time required to form the distribution with high gradients of probability.
The different choices of parameters and initial conditions reveal that the characteristic of the Liouville solution for this simplified system are fairly sensitive to the model parameters. The direction of the highest Fig. 7. Time evolution of a segment of high probable states (X 1 = [0.1 − 0.14], X 2 = 0.24, X 3 = 0.1; in blue) and low probable states (X 1 = [0.1 − 0.14], X 2 = − 0.24, X 3 = 0.1; in red) and projections along the three two-dimensional plains: a) t = 0, b) t = 100, c) t = 250, d) t = 500, e) t = 1000, f) t = 5000. A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 variability of probability density changes with the ratio between wavenumber in x and y directions, δ, obtaining higher gradients of probability in the X 1 direction for higher values of δ. The β-effect controls the period of the full PDF obtaining a periodic solution when this effect is considered. On the other hand, the sensitivity of the displayed characteristics to the initial condition is low.

Solution in physical space
For the sake of a better understanding of the implications of the solution in physical space, the characteristics in physical space of states with high probability are investigated in detail for the β = 0, δ = 2 case. States with X 3 close to zero present streamfunctions which are symmetric with respect to maxima and minima, while increasing X 2 Fig. 8. Normalized values of probability density (ρ max =1) on plane X 2 =0 for t = 85 (a, b), t = 200 (c, d), t = 500 (e, f) and t = 2000 (g,h). Panels b, d, f and g display an extended range of X 1 . A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 generates streamfunctions with sharper and narrower maxima and minima. For a given X 2 , the effect of changing X 3 is the rotation of the symmetry axis. For X 3 = 0, this axis is parallel to y axis. X 1 controls wave amplitude. An approximation to the mean and standard deviation of the coefficients is computed for different times to associate the shape of the probability density function to uncertainties in physical space. At t = 0 (Fig. 13), with the same spread for all three coefficients, the spread of X 1 translates onto larger spread of the streamfunction field, whereas the spread in X 2 and X 3 causes only minor changes on the waves at maxima and minima, as well as in the space in between, respectively. However, for t > 0, the spread in X 1 remains approximately constant while the spread in X 2 and X 3 increases. This spread of X 2 and X 3 displays fields which are very different from the mean. Remarkably, the mean is a state roughly at the center of the spiral (Fig. 4) with a very low value of ρ.
The reflect of uncertainty in X 1 , X 2 and X 3 in physical space is assessed by plotting various states with high probability at t = 100 (Fig. 14). At the time the ring is forming (Fig. 4b), uncertainties in X 2 and X 3 are large because points with high values of ρ are found for very different values of X 2 and X 3 and consequently, the corresponding streamfunctions are also very different. The main differences between them are the position of maxima and minima and the streamline distribution. Points with the highest values of X 2 and with X 3 close to zero have closed streamlines in the central part of the depicted domain in the y direction and lines of constant streamfunction lay symmetric with respect to the maxima and minima (Fig. 14b). The value of X 3 modifies the orientation of the closed streamlines. For larger values of |X 3 |, these closed streamlines are more elongated and the axis of symmetry does not cross maxima and minima (Fig. 14a, c, d). For small |X 2 |, no closed streamlines (Fig. 14e) emerge.
For t = 5000, with the filamented distribution well formed (Fig. 4d), we investigate in detail the state of the system for large and small values of X 2 and X 3 (Fig. 15). Fields with X 3 close to zero and varying X 2 ( Fig. 15a and b) or vice versa ( Fig. 15c and d) have different locations for maxima and minima. The difference between the two couples of fields is the shape of streamlines with the emergence of closed streamlines in the first couple.
Unlike X 2 and X 3 , the most probable values of X 1 are close to the initial X 1 and the spread is small. Despite the remarkable differences in physical space between states with high probability, the evolution of system (10) only gives rise to a few number of streamfunctions with specific combinations of position, amplitudes and orientations. Although states with very different characteristics have high probability, some other regions of phase space are unreachable for this combination of system dynamics and initial condition.

Summary and discussion
The Liouville equation has been solved analytically for some low complexity systems by using the general solution obtained with the method of characteristics. The analytical solution describes the evolution of the PDF of a dynamical system generated by an imperfect knowledge of the initial conditions. The different examples considered in this study exhibit interesting characteristics.
The PDF evolution for the shallow water solution illustrates the effect of divergence (convergence) in phase space trajectories, which leads to a decrease (increase) in predictability. As one might intuitively expect, high probable states for the shallow water system remain compact for all times, so that a discrete sampling of the distribution function in phase space (i.e. ensemble of evolutions) will render a satisfactorily precise representation of the PDF evolution.
More striking solutions are obtained for the barotropic model. The Liouville solution for this system displays a PDF with high gradients of probability. Although individual trajectories are periodic in phase space, the full evolution of the PDF is more complex, rendering distributions with highly probable states nearby low probability states. This would provide misleading interpretation of finite-size ensemble predictions as high probability regions are separated by states with very low probability. In the context of ensemble prediction systems, and under the ever-present constrain of limited availability of computational resources, the sampling strategy of the initial PDF targets highly probable states only. However, standard postprocessing and interpretations of ensembles would assume a smooth underlying distribution in which states in regions between ensemble members are also probable, which would be misleading if the actual PDF presented high gradients of probability. This effect worsens when the high variability is not uniform across the PDF, so that its support combines smooth and connected regions of high probability with others with high gradients of probability. The impact of this effect is critical when the scale of the event of interest is of the same order than the regions in phase space with high and low probability. On the contrary, the presence of low probability regions embedded in large volumes of interest is not problematic because underestimations of high probabilities are compensated by overestimations of low probabilities. The cases which the scale of an event of interest is of the same order of the variability of the PDF would pose fundamental challenges on current ensemble prediction systems as their common interpretations are based on the hypothesis of smooth forecast PDF. In this context, this misconception would alter probabilistic forecasts as probabilities of unlikely states would be overestimated at the expense of surrounding highly probable states.
A further important feature of the solution for the analytical barotropic case is its dependence on structural parameters δ and β. The former, which is the ratio between wave numbers in x and y, affects the gradients of the solution, producing a solution with higher variability in one of the phase space directions for higher values of δ. Conversely, β, which accounts for a linear variation of the Coriolis parameter, generates a periodic PDF solution when is different from zero, and a connected support PDF for high β values. This evidences that PDF evolutions can display multiple topological characteristics depending on the governing dynamical systems.
The extension of this study to explore the features of the support for more realistic systems with larger number of degrees of freedom (e.g. operational or research weather models) is an open question as it requires unforeseeable amounts of computational resources. Even for lower complexity models, the system dimension is excessively high to apply the Liouville formalism without making significant Fig. 9. Ratio between probabilities computed from a smoothed PDF (p sm ) and from the actual PDF (p) over a range of volumes in phase space centered at the maximum value of ρ for t = 50 (solid), t = 100 (dashed), t = 250 (dotted) and t = 500 (dash-dotted).
A. Hermoso, et al. Atmospheric Research 246 (2020) 105127 approximations to reduce the dimensionality of the system. Admittedly, properties identified in simplified systems may not be present in real world applications. Indeed, diffusive processes in more complex systems may produce smoother PDFs. In any case, current limitations to solve the LE prevent these issues to be thoroughly addressed. Furthermore, the formalism described in Section 2 does not apply to the infinite partial differential equations problems, but to the discretized finite dimension version system, which is the approximation to the true dynamics considered by weather and climate prediction models. Regardless of the above comments, the first quantitative evidence of high variability in the solutions of the LE, described and characterized in this study, poses a fundamental challenge on current probabilistic geophysical forecasting systems as it compromises the hypothesis of compactness and smoothness assumed by most current ensemble interpretation and postprocessing tools. The identification of these characteristics in probability distribution functions of realistic geophysical systems paves the ground towards more precise, accurate and valuable forecasts.

Declaration of Competing Interest
None. The initial conditions in s necessary to obtain this solution come from the initial condition of Eq. (A.1), which is given by the general form f(x, u) = 0. This surface is taken as s = 0. If x and u depend on {s, t 1 , t 2 , …t n }, variables {t 1 , t 2 , …t n } can be used to parametrize the initial conditions: The solution of (A.6) is t = s, the set of Eqs. (A.7) are the equations of the dynamical system, whose solution is required to obtain an explicit solution of Eq. (A.8). This is the result presented in Section 2:  Fig. 4 for β = 2.5 ⋅ 10 −7 at: a) t = 2nπγ −1 , b) t = 2nπγ −1 +50, c) t = 2nπγ −1 +100, and d) t = 2nπγ −1 +200 with n an integer. A. Hermoso, et al. Atmospheric Research 246 (2020 where = − f γ 1 2 2 and γ, ϕ, W and H 0 are parameters to be defined by the initial conditions U 1 0 , V 1 0 , A 0 and H 0 . These parameters have limited ranges of application: −1 < γ < 1, 0 < W ≤ f 2 , H 0 > 0. ϕ is a phase shift in the solution determined by the initial condition time. For simplicity, we take ϕ = 0 at t 0 = 0. Furthermore, γ, W and H 0 are determined by the initial conditions: Fig. 13. Representation of mean field (black lines) and mean field plus (dashed grey lines) or minus (dotted grey lines) two standard deviations in physical space for: a) X1, b) X2 and c) X3 at t = 0 for ρ 1 .