Regular Liesegang patterns and precipitation waves in an open system

We investigate the regular and moving Liesegang pattern formation phenomena in an open system.

First, simulations have been performed at fixed coupling between the reactive medium and the reservoir, later this control parameter was varied during the simulations resulting in various phenomena.

We predicted and monitored for the first time various—dynamically changing—precipitation structures and a spatial hysteresis phenomenon, which is beyond the scope of the Turing instability.

The dynamics of the reaction is well detectable using specific quantities: the total amount of precipitate and its center of gravity.


Detailed investigations of homogeneous chemical systems showed that phenomena like bistability, complex and chaotic oscillations may exist only in open chemical systems, far from thermodynamic equilibrium.1,2,3

The appropriate experiments are carried out using a continuously stirred tank reactor (CSTR).

In 1952, Turing predicted that the diffusion of species may destabilize a homogeneous chemical system producing stationary patterns in a spontaneous way.4,5

The pioneering work of De Kepper and his coworkers provided for the first time an experimental evidence for the formation of Turing patterns.6

Up to the last decade, a wide range of spatial or spatiotemporal reaction–diffusion patterns (such as complex Turing patterns, phase waves, autocatalytic chemical fronts, etc). have been designed and investigated also in open systems. 7,8

Formation of the precipitation patterns in the wake of diffusion driven chemical fronts has been extensively studied from the end of the 1800’s.9,10

This is the first chemically designed and controlled spatiotemporal pattern formation.

Although the generation of such patterns (called Liesegang patterns) is quite easy, the research on this topic still has a great relevance.

In the typical experimental setup, one chemical reagent is uniformly distributed in a gelled medium (called inner electrolyte, B), while the other one (called outer electrolyte, A) diffuses from outside.

The initial concentration of A is chosen to be much larger than that of B. This condition ensures the higher diffusion flux of the outer electrolyte into the gel.

The coupling of the diffusion and the precipitation reaction between the two chemical components produces a non-soluble precipitate (AB) in the reaction front.

The distribution of the precipitate is usually non-uniform behind this front.

Appearing precipitation zones (Liesegang bands or rings, depending on the geometry of the experimental setup)11–19 have several empirical regularities (time, 20,21 spacing,22 and width law23–25), in which some well-measurable quantities are connected with each other.

These are Xn, the position of the nth band, measured from the junction point of the electrolytes, tn the appearance time, and wn, the width of the nth band, respectively.

The Matalon–Packter law describes the dependence of the spacing coefficient (P = Xn+1/Xn) on the initial concentration of the inner and outer electrolytes.26–28

Recently, a new and more general law has been proposed, which also includes advective dynamics in the reactive medium.29

For this, two new quantities are introduced: ptot, the total amount of precipitate and Xc, the position of the center of gravity of precipitation system measured from the junction point of the electrolytes.

Besides the theoretical derivation, experimental and numerical evidences confirmed the new law, which states that ptot is linearly proportional to Xc.29

From the beginning of the last decade, more and more complex precipitation pattern formation scenarios have been studied and discussed.30,31,34–36

Several authors reported that the precipitation process followed by complex formation of precipitate may produce moving Liesegang patterns (“moving” bands)34–36 or moving precipitation pulse (precipitation waves).12,30–33

In several papers, deterministic chaotic time oscillation of the total number of bands,35,37 effect of an electric field on pattern structure,12,38 and crossover from the precipitation waves to moving Liesegang patterns39 have been investigated.

These are the recent trends in the study of the precipitation pattern formation phenomena.

Das et al. presented the Liesegang pattern formation in 2D using a gel-ring reactor, in which the two electrolyte solutions were fed from reservoirs13: The outer electrolyte was fed continuously from the center, while the inner one from a gel ring outside the reaction medium.

Several years ago Bhattacharya et al. proposed a model, which included a single autocatalytic reaction step followed by nucleation and growth.40

Based on the fact that a cubic autocatalytic step coupling with CSTR may give rise to Turing patterns,41 they demonstrated formation of spatially periodic precipitation patterns.40

In classical experimental setup, the diffusion front of the outer (invading) electrolyte depletes the inner electrolyte from the reactive medium due to the precipitation reaction.

The aim of this paper is to provide for the first time a systematic investigation of the regular and moving Liesegang pattern formation in an open system, where a continuous feeding of the inner electrolyte is allowed.


To illustrate our idea, we have chosen the simplest precipitation mechanism (1 : 1 type, AB), giving rise to precipitation patterns, according to the reactionA(aq) + B(aq) → AB(s),where A and B yield the outer and inner electrolytes, respectively, while AB is the precipitation product.

In addition, the precipitate reacts with the excess of the outer electrolyte, producing a soluble complexAB(s) + A(aq) → A2B(aq).For simplicity, we supposed that the formed complex A2B is stable, so the reverse step was neglected.

Due to the irreversible step the stable complex does not affect dynamically the pattern formation.

Therefore, it was supposed that neither the complex nor the precipitate diffuses.

A sketch of our system appears in Fig. 1 where for a clear visual interpretation we provide a 2-D description.

In fact, the thickness of the reactive medium is negligible and therefore, the reaction is described using only horizontal coordinates.

The dimensionless governing equations for the reactive medium, coupled with a well-mixed reservoir for the inner electrolyte (open system) in 1-D arewith τ and x being the time and length coordinates, respectively.

Here a, b, c, and p denote the concentrations (or even density for solid phase) of the outer and inner electrolyte, the complex, and the precipitate, respectively.

Da and Db are the diffusion coefficients of the corresponding species, while b0 is the concentration of the inner electrolyte in the reservoir (CSTR).

k is the reaction rate constant for the complex formation (eqn (2)).

κ is related to the coupling between the reactive medium and CSTR, which maintains a continuous mass flux of the inner electrolyte.

κ is the inverse of the residence time of species in the reactor, which can be accurately controlled by the flow rate.

The function Δ(ab, K, L) (will be specified later on) describes the precipitate formation based on Ostwald’s supersaturation theory.

Despite its simplicity and generality this concept provides a perspicuous mechanism in many cases, where the precipitation process plays a partial role.

The efficiency of the above mentioned precipitation mechanism was tested by simulating the effect of an electric field,42,43 stochastic precipitation pattern formation,44 and also moving Liesegang bands and precipitation waves.37

The term kpa in eqns (3), (5) and (6) represents the complex formation rate, whereas κ(b0b) in eqn (4) describes the coupling between the reactive medium and the reservoir.

However, it should be noted that the immobilization of the complex is not a necessary condition for the moving pattern formation.

To describe the precipitation formation presented by the simple mechanism (1) we choose the model proposed by Büki et al.,45,46 which is based on the Ostwald’s ion-product supersaturation theory.47

The precipitation reaction term Δ(ab,K,L) is defined as follows: where k0 is the rate constant of the precipitation reaction (1), K is the nucleation product, L is the solubility product, and Θ is the Heaviside step function.

The coefficient Sp is related to the extent of the supersaturation of the system and influences the intensity of the reaction.

Formally, Sp is calculated using the equation (aSp) (bSp) = L. As far as abL we do not need to define it, since the precipitation reaction term in (7) is zero, anyway.45,46

The precipitation reaction continues at a certain point x until the concentration p in x has reached a threshold pmax.

Numerical simulations were performed using a “method of lines” technique.

Discretization of eqns (3)–(6) was carried out on a 1-D equidistant grid applying a second-order finite difference scheme.

In this procedure, at any fixed time the functions a,b,c and p in (3)–(6) are substituted with vectors [a(τ, x1), a(τ, x2),…,a(τ, xn)] and similarly for b, c and p, respectively, where x1,x2,…,xn are the grid points.

The spatial differential operator turns to be a matrix of dimension n × n and with this, (3)–(6) can be approximated by a system of ordinary differential equations (of size n) with respect to time.

This has been solved using a second-order Runge–Kutta method with the following initial conditions:a(0,x) = a0Θ(−x), b(0,x) = b0Θ(x), p(0,x) = c(0,x) =0.The boundary condition is no flux type both for the inner and the outer electrolyte at the end of the reactive medium.

It is Dirichlet type a(t,0) = a0 for the outer electrolyte at the junction point of the electrolytes.

Although both κ and b0 could be used as control parameters of the system, b0 was fixed in the simulations.

This concentration has been chosen to coincide that of the initial concentration of B in the reactive medium.

The presented model allows to describe the crossover from Liesegang patterns to precipitation waves using only one parameter: in case of k = 0 regular Liesegang pattern formation, and for k ≠ 0 moving patterns are expected.

In all simulations, we used the parameter set Da = Db = 1.0, K = 0.15, L =0.1, k0 = 1.0, and pmax = 50.0.

The grid spacing was Δx = 1.0, and the numerical simulations were performed with the time step Δτ = 0.01.

Results and discussion

Regular and moving precipitation patterns at fixed coupling

First, the dependence of pattern structure on the coupling coefficient (κ) has been investigated in the case of regular Liesegang patterns (k = 0), and precipitation wave (k = 10−3).

Obviously, if κ = 0 (there is no coupling between the reactive medium and reservoir), then regular Liesegang pattern and precipitation wave evolve, and the final pattern will depend on the particular parameter set.

If κ ≫ 0, a single precipitation zone develops, starting from the junction point of the electrolytes instead of the (regular or moving) Liesegang patterns, due to the high flux of the inner electrolyte from the reservoir.

This extremely high coupling reduces and slightly modifies the purely diffusive profile of the invading (outer) electrolyte A. This intensive coupling ensures that the product of the local concentrations is kept above the nucleation product (K), therefore, the precipitate is formed continuously.

Fig. 2a,b shows the final pattern structure for fixed values of κ.

Results of the numerical experiments are in good agreement with the above explanation.

In the case of the precipitation waves, increasing κ results in a thicker precipitation zone, the final position of which is getting more and more close to the junction point.

Increasing κ promotes the precipitation reaction, thus reduced concentration of the outer electrolyte puts back the complex formation.

In contrast to the case of precipitation waves, we predict regular Liesegang patterns for small values of κ, and a single precipitation zone for relatively high values of the control parameter.

The time law for the regular Liesegang patterns (even for 1-D or planar ones) states that the position of each band is linearly proportional to the square root of its formation time.

This is a consequence of the diffusive dynamics of the invading electrolyte.

In several cases, modification of the dynamics (including the presence of external electric field, advection field, curvature effect) results in some deviations compared to the time law.20

Fig. 3a presents the evolution of the pattern for various values of κ.

Relatively high coupling hampers the depletion of B, therefore, the more intensive precipitate formation reduces the concentration of A in the reactive medium, which leads to the slower evolution of the pattern.

Similarly, the center of gravity of the precipitation system (Xc) exhibits slower motion in time as depicted in Fig. 3b.

The total amount of the precipitate ptot and the center of gravity of the precipitation system Xc are defined as follows:where L is the length of the 1-D reactive medium.

At the same time, the enhanced reaction rate increases the total amount of precipitate (ptot), the evolution of which has been plotted in Fig. 3c.

Note that increasing the coupling coefficient affects more significantly ptot than Xc.

In this way increasing κ, dependence of ptot on Xc exhibits an accelerating growing rate (Fig. 3d).

In the regular 1-D case, in all of the subfigures of Fig. 3 one can observe a linear dependence.

Moreover, this linear dependence in Fig. 3d is also valid in the case of an additional advective flux in the system29.

Moving precipitation patterns using time dependent coupling

Up to this time, the control parameter (κ) was fixed at certain values during the numerical experiments.

Now we will focus on the behavior of the system with the coupling rate between the reactive medium and reservoir being varied on a continuous scale: various time dependences of κ have been applied.

Accordingly, variation of the total amount of the precipitate (ptot) and the center of gravity of the precipitation system (Xc) is depicted in Fig. 4 as a function of log κ.

The change of log κ was linear on [nT, (n + 1)T] with integer values of n.

We give the exact time dependence of log κ as follows:which allows to illustrate clearly its effect on the distribution of the precipitate, where T is the time period of one cycle.

All simulations were started at κ = 10−9, which corresponds to a very weak coupling.

Increasing κ results in an increased amount of the precipitate due to the enhanced mass-flux of the inner electrolyte from the reservoir into the reactive medium.

The maximum of ptot is attained at the maximal rate of the coupling (κ = 1 in the simulations).

Similarly, ptot is decreasing when the coupling coefficient κ is decreased (Fig. 4a).

One can observe a slightly different behavior of the center of gravity with varying κ.

The position Xc is an increasing function of κ up to a certain coupling coefficient, afterwards, Xc is decreasing (Fig. 4b).

Decreasing κ from the maximum to its minimum value, Xc monotonically increases.

At the first stage, when the coupling is weak, a precipitation wave travels in the reactive medium.

In the case of higher coupling, the precipitate formation behind the pattern will dominate compared to the precipitate dissolution by complex formation.

At the same time, the precipitate formation is stopped in front of the pattern, because the coupling and precipitation reaction behind it reduces the concentration of the invading electrolyte.

Therefore, Xc is shifted backward, however, the total amount of precipitate increases continuously.

The complex formation of the precipitate has a more significant role behind the pattern, when the coupling rate is decreased.

The same applies to ptot, while Xc increases.

At weak coupling, the precipitation process continues in front of the precipitation wave and the wave evolves further.

The most interesting scenario observed during the simulations is the hysteresis of the pattern.

Accordingly, one can observe the hysteresis of the quantities ptot and Xc related to the pattern.

The time oscillation of ptot is a trivial consequence of the variation of κ in time, but the hysteresis is a unique feature.

Note that similar hysteresis curves (with different maxima) belong to different time dependent functions of κ.

Hysteresis occurs due to the interaction of the reaction rate (for the precipitation and complex formation), diffusion of the electrolytes, and the coupling rate.

When increasing κ, the dominant process is the precipitate formation behind the pattern due to the high coupling, which reduces the diffusion of the invading electrolyte.

When decreasing κ, the precipitate formation process diminishes behind the pattern, and the evolution of the pattern will be driven by complex formation.

This process may be accompanied with the precipitate formation in front of the pattern, if the diffusion profiles of the outer electrolyte are regenerated.

In the simulations, we varied κ periodically over some periods (Fig. 5).

We observed that the values of ptot corresponding to the maximal and minimal values of κ have been increased, compared to those of the first period.

An interesting feature (using a given parameter set) of the observed system is that after the first loop, a periodic change of κ results in the periodic change of total amount of the precipitate.

Analyzing the evolution of Xc while varying the time dependent κ in the same way, a similar behavior of the center of gravity has been found.

Fig. 6 presents the evolution of the precipitation wave over four periods.

Position of the front and back section of the wave and the center of gravity exhibit a forced time oscillation according to the change of κ.

Time dependence of the back section of the wave has different dynamics whenever the variation of log κ (in time) was linear, while decreasing or increasing this control parameter.

Note that the thickness of the precipitation zone is related to the total amount of the precipitate, therefore both of them exhibit the same behavior.


In this study, we have investigated the regular and moving precipitation pattern formation in an open system using a CSTR for the coupling of the inner electrolyte between the reactive medium and reservoir.

Evolution of patterns in both cases, and the distortion of the existing laws (time and “universal” laws) due to the coupling coefficient have been investigated.

The coupling of the inner electrolyte between the CSTR and the reactive medium has a great impact on the pattern structure.

Variation of the coupling coefficient (κ) during the individual simulations resulted in a unique behavior for the moving precipitation wave phenomenon: the pattern structure and the corresponding characterizing quantities, the total amount of the precipitate (ptot) and the center of gravity (Xc) have hysteresis.

This new result is a unique observation on this research field.

Compared to the Turing patterns which are stationary ones, the above investigated moving precipitation wave phenomenon in an open system represents dynamic patterns.

The spatial hysteresis of the Turing patterns has been known for more than one decade.

The novel result of our study is that we observed the hysteresis phenomenon in a precipitation system, which qualitatively differs from activator–inhibitor or cubic autocatalytic systems, where Turing patterns have been reported.

The coming step in our study is to design the experimental setup to validate our new predictions.