What is the difference between isoviscous and layered viscosity models




















For viscosity contrasts greater than the critical value a sluggish-lid solution exists and decreasing the LVL viscosity decreases the plate velocity. Results for case StPlate. With , energy needs to be supplied to overcome the net dissipation in the plate. This energy can be provided in three different ways, each corresponding to one of the three boundary layer solutions. For low viscosity contrasts, the mantle pressure gradient provides the additional energy.

However, for a very low LVL viscosity, the LVL cannot support a large pressure gradient and the contribution from the pressure term becomes too small. Basal tractions can provide the required energy to drive the plate. But in order to generate a negative shear stress the plate must be moving slower than flow in the mantle and a sluggish-lid solution is obtained. Finally, the plate can be driven by the local rate of change of potential energy in the lithosphere.

This term is small and again results in a slow moving sluggish-lid solution. For very low LVL viscosities, the traction and pressure terms are too small and plates are driven by the local rate of change of potential energy, leading to slow moving sluggish lid solutions. This example demonstrates that while a LVL can promote plate tectonics, it can also inhibit plate motions by reducing the available energy from mantle coupling.

Thus, a sluggish-lid solution does not necessarily imply a slow moving and thick plate, as was the case for the isoviscous solutions considered in Section 3. For very large viscosity contrasts, and thus very low , the model predicts either a return flow in the LVL with fast moving plates when or a channelized flow in the LVL with slow moving plates when.

Thus, the sign of is critical in determining the behaviour of the system. The results of Figs 6 and 7 demonstrate that the difference in dynamics may only be noticeable for large viscosity contrasts.

It is also worth noting that most of the terms in the lithospheric energy balance are of comparable size. The behaviour is complex and is best captured by solving the full energy balance equations. Simplified scalings that balance only two terms will not capture these dynamics.

While the intermediate and lower branches of solutions are theoretically and conceptually viable steady-state solutions they satisfy the energy and fluid flow equations , the stability of solutions on these branches has not yet been determined.

As such, it is reasonable to ask whether there is any evidence of these forms of solutions in nature or in the large body of literature on numerical mantle convection simulations.

All outer boundaries were free slip. Their model had a simple depth-dependent viscosity structure with a high viscosity lithosphere upper 10 per cent , a low viscosity channel for the asthenosphere 10 per cent , and an intermediate viscosity lower mantle lower 80 per cent.

Their numerical simulations focussed on the effects of changing aspect ratio. Their results demonstrated a channelization of flow in the asthenosphere at low aspect ratios and a rapid increase in plate velocity with increasing aspect ratio. They noted that two scaling regimes exist, one at low wavelength and one at larger wavelengths.

They proposed that the negative buoyancy increases with cell wavelength and at the regime transition it overcomes the bending resistance. Comparison to numerical simulations. Inset figures show the rms average of the horizontal velocity profiles through the convective cell. Dashed lines indicate negative values in the energy balance work done by the lithosphere while solid lines indicate positive values work done on the lithosphere.

As in their study, we use dimensionless parameters , , , , , , , ,. For simplicity we assume that internal heating is the main mode of heating and we use a dimensionless heat flux of Again, the depth averaged advective heat flow through both the convective cell and the lithosphere are assumed to be. In all previous calculations the thermal and mechanical thickness of the lithosphere was assumed to be the same. For this case they are independent and we define the thickness of the lithosphere, , using the fixed mechanical thickness.

Their simulations did not include lateral viscosity variations. As such, the slab has little strength in the asthenosphere and we set the slab-pull normal stress, , to zero. A consequence of the absence of slab-pull is that deformation in the lithosphere is not focused at the ends of the cell and can be diffuse [see surface velocity in fig. As a result, the plate does not deform like a bending viscous beam and we treat dissipation in the lithosphere by using the constant stress model given in eq.

The characteristic stress for the lithosphere is , where is the mantle velocity a characteristic velocity for the problem , is the effective viscosity of the lithosphere and is the thickness of the lithosphere. From their simulations, we use for the characteristic mantle velocity. Figs 8 b and c show flow profiles from their numerical simulations black lines for convective cells with aspect ratios of 4 and 10, respectively [these are the same flow profiles shown in fig.

A channelized flow occurs in the asthenosphere for small wavelengths. For larger wavelengths the asthenosphere acts as a lubricating layer for the plate and the flow in the asthenosphere is Couette-like.

The inset plots in b and c show the RMS horizontal velocity averaged over the convective cell for the numerical simulations black line and calculated from our model red line. Given the simplicity of the our analytic model, the agreement is quite good. The analytic model slightly overestimates mantle flow rates for the small aspect ratio cell. Nonetheless, the model does a reasonable job at getting the plate velocity and characteristic flow profiles correct.

For small aspect ratios, the short and low viscosity asthenosphere cannot support a large pressure drop and work done by the mantle pressure gradient is small.

Both the magnitude of the shear stress on the base of the plate and the surface area over which the stresses can be applied are small for a low viscosity asthenosphere and a small aspect ratio convective cell. Thus, energy provided by basal tractions is also very small.

All that remains to drive the motion of the plate at small aspect ratios is the local rate of change of potential energy in the lithosphere. This term is small and results in a solution that is a slow moving self-driven plate. This is another example of how the presence of an asthenosphere can inhibit plate motions by effectively decoupling the plate and mantle.

As the aspect ratio of the convective cell increases, so does the length of the asthenosphere. A longer channel can, for similar flow rates, support a larger pressure drop across the mantle and work done by the mantle pressure gradient is larger for large aspect ratio cells.

Similarly, a longer channel provides more surface area for tractions to act on the lithosphere and the contribution of tractions to the lithospheric energy balance is larger for large aspect ratios. The boundary terms drive the plate at large aspect ratios and result in a solution that is a fast moving mantle-driven plate.

Our analytic model predicts that tractions help to drive the plate at low aspect ratios but can provide drag that inhibits the motion of the plate at larger aspect ratios. At this aspect ratio the flow in the asthenosphere changes from mainly Poiseuille flow to mainly Couette flow. Furthermore, our energy balance results suggest that the transition between the two different scaling regimes is not the result of an increase in plate buoyancy relative to plate resistance which decreases , but is instead due to the transition from a sluggish-lid solution at small aspect ratios driven by the local buoyancy of a slow and thick lithosphere to a mobile-lid solution at large aspect ratios driven by mantle coupling via the mantle pressure gradient.

The Earth is a complex fluid dynamical system that has many plates, three dimensional geometry, and time dependence. Neighboring plates interact and transfer energy to each other through boundary forces. By comparison, our analytic model represents a much simpler system that has a single plate and only two dimensions. As such, it would be tenuous at best to use our model to try to predict plate velocities for specific plates on the Earth.

Nonetheless, if the first order physics are correct, then our model can be expected to predict average values for plate velocity, heat flow, plate thickness, etc.

Increasing the effective lithospheric viscosity decreases the plate velocity. Increasing the amount of slab-pull increases the plate velocity. This is all consistent with observations and our current understanding of present day plate tectonics on the Earth. The magnitude of the plate velocity in this sluggish lid regime would likely have been similar to the magnitude of present day plate tectonics.

This would have slowed the cooling of the Earth and regulated plate tectonic rates through time. Eventually, due to a slow cooling of the mantle and gradual increase in viscosity, the mantle viscosity would become large enough that mantle flow could more efficiently couple to the plates. This would cause plate tectonics to transition from the sluggish-lid regime to the mobile-lid regime, increasing heat flow. This has significant implications for the thermal evolution of the Earth and for the present day Urey ratio ratio of internal heating to heat loss from the mantle.

It also has significant implications for the geochemical evolution of the Earth, as chemical differentiation of the mantle at ridges and the creation of large scale heterogeneity as plates are subducted at trenches are both modulated by plate tectonic rates. A detailed study of the thermal evolution of the Earth, with a heat flow scaling that allows for the full stagnant-lid convective regime, will be the focus of a separate study.

The sluggish-lid solutions presented in Section 3. The thickness of these plates, through the dependence in eq. In fact, solutions on the lower branch plotted in Fig.

For plates this thick, it is likely that other physical processes, such as small scale convection Korenaga , would take over and might limit the thickness of the plate. However, Section 3. Moreover, it would be premature to dismiss these solutions, regardless of their plate thickness, as thick plate solutions might be possible under different situations, such as early Earth or on other planets.

While our analytic model makes no claims with regard to the episodicity or any time dependence of convection, it does predict that a region of multiple solutions can exist see for example Fig.

Episodic convection could occur in this region as the system transitions between the different branches of solutions due to either perturbations in the system irregular or a natural hysteresis periodic. Such behaviour might be explained by jumps between the various solution branches. Transitions between branches could be smooth, such as the transitions in Figs 6 and 7 , or could be abrupt, as suggested in Figs 4 and 5. This behaviour will largely depend on the stability of the different solution branches and this will be the subject of future work.

It is worth noting, however, that the sluggish-lid solutions shown in Fig. This suggests that at least some of the sluggish-lid solutions predicted by our model are stable.

Such a numerical counterpart to our analytic solution provides some vicarious verification. Multiple solutions also allow for the possibility that two planets with similar properties, such as the Earth and Venus, could be on different solution branches and therefore exhibit very different convective states.

The particular state of the system, in a region of parameter space with multiple solutions, could depend on the evolution of the system leading to that region of parameter space. Such path dependency may be fundamental to our understanding of planetary evolution. This is demonstrated in Fig. Furthermore, there are regions of parameter space where only the sluggish-lid boundary layer solution exists. Thus, while the classic boundary layer scaling may capture the dynamics of present day plate tectonics on the Earth, it should not be blithely applied to the early Earth, where conditions may have been very different, or to other terrestrial planets with different properties, sizes e.

Our analytic model provides a tool that may be used to explore and map out the parameter space for the different solution branches. We have kept the model as simple as possible in order to keep the physics clear and we recognize that it can be improved and made more accurate. Our aim has been to determine the nature of the forces driving plates and to be able to extract simple analytic expressions that capture the model behaviour and reveal the effects of various parameters and features of the model.

We have explored solutions where the lithospheric thickness is controlled by the thermal thickness. Other models for predicting lithospheric thickness have been proposed. Our model allows for the simple substitution of other parametrizations for plate thickening and dissipation and provides a general framework for studying a simple convection cell with a finite-strength plate. We have used simple representations of the flow in the model.

This could be made more accurate by using eigenfunctions of solutions of the Navier Stokes equations, as in Busse et al. We have also limited our analysis to a conductively defined lithosphere with a thickness that can increase until it subducts. Our model is steady state, and does not address the formation of subduction zones and ridges that are necessary for plate tectonics, which is an important problem in its own right.

The depth of the frictional resistance at a subduction zone is taken as the plate thickness, although other choices could be made from considerations of the mechanics and state of material subducted. These may change the results somewhat, but as our analysis shows eq.

Little attention was given to the effect of slab-pull in this study. Slab-pull is an important driving force for plates Wu et al.

However, we demonstrated in Section 2. The presence of slab-pull simply reduces the effective yield stress or bending stress of the plate. Thus, explicitly including slab-pull in the model will not change the behaviour of the system.

Including the lithospheric energy balance is important for properly predicting the plate velocity. The plate velocity modulates the heat flow out of the system eq. As the rate of change of potential energy scales with the heat flow eq. Thus, even though the magnitude of the lithospheric dissipation term may be small in the global energy balance, the total energy available for driving mantle convection is determined by the plate velocity and thus by the lithospheric energy balance.

As the strength of the plate is increased the model deviates from the classic parametrization and convective solutions in the sluggish-lid regime become possible. The strength of the lithosphere is determined by the effective lithospheric viscosity and radius of curvature for the plate bending model, or the plate yield stress for the shearing model.

The two different models for calculating lithospheric dissipation are equivalent and predict constant stress when the thickness of the lithosphere is determined by its thermal boundary layer thickness for fixed aspect ratio.

Thus, all of the results presented are independent of the choice of model for calculating lithospheric dissipation. The contrast between the strength of the lithosphere and the strength of the mantle determines the convective regime of the solutions. Solutions in the mobile-lid convective regime occur when either the plate is weak or the mantle viscosity is high and the contrast in strength is low.

Solutions in the sluggish-lid convective regime occur when the contrast in strength is high and either the plate is strong or the mantle viscosity is low. Multiple solutions for the plate velocity are present for intermediate contrasts in strength between the lithosphere and the mantle.

The upper branch of plate velocity solutions corresponds to the classic form of convective cell in which a simple cellular flow carries a weak plate on the surface. This mode of convection is controlled by the material properties of the mantle.

The intermediate branch of plate velocity solutions represents a convective cell in the sluggish-lid regime with a plate velocity that is less than the mantle velocity. Basal tractions are the dominant driving force on the plate and the dynamics for this branch depend on both the mantle and lithosphere material properties.

The lower branch of plate velocity solutions represents a mode of convection in which the lithosphere controls the dynamics of the system. The local negative buoyancy of a thick and strong lithosphere drives the motion of the plate. Work done by basal tractions and any mantle pressure gradient are negligible compared to the rate of change of potential energy and local dissipation in the lithosphere.

This mode of convection occurs for strong plates and can exist with little slab-pull. Furthermore, it is controlled by the strength of the lithosphere and is independent of the mantle properties. The lithosphere modulates the heat flow out of the convective cell, regulating the amount of energy available to drive mantle convection, and thus also controlling flow rates in the mantle. The lower branch of plate velocity solutions exists only for plates with a positive net resisting stress and when the contrast in strength between the lithosphere and mantle is high.

The introduction of an asthenosphere LVL can significantly change the dynamics by altering the plate-mantle coupling. Lowering the asthenospheric viscosity, relative to the lower mantle, can promote plate motion by providing a lubricating layer. However, a very low viscosity asthenosphere beneath a strong plate can inhibit plate motion and produce a slow moving plate with a channelized flow in the asthenosphere.

Our model is able to capture and explain the transition in plate driving forces responsible for this behaviour. Thus, an asthenosphere does not always promote plate tectonics. The distinct multiple solutions that occur are all energetically equivalent.

Thus, the state of the system may depend on its history [this idea of transitioning between different modes of convection is discussed by Sleep ]. It is also possible that the system may transition between the various states due to natural perturbations in the velocity and density fields that occur in high Rayleigh number thermal convection flows. Examining the stability of the various solution branches will be the subject of future work.

We wish to thank the editor Bruce Buffett as well as Clint Conrad and another anonymous reviewer for their positive and constructive comments. Booker J. Thermal-convection with strongly temperature-dependent viscosity , J.

Fluid Mech. Google Scholar. Buffett B. Plate force due to bending at subduction zones , J. Busse F. Richards M. Lenardic A. A simple model of high Prandtl and high Rayleigh number convection bounded by thin low-viscosity layers , Geophys. Capitanio F. Morra G. Goes S. Dynamics of plate bending at the trench and slab-plate coupling , Geochem. Conrad C. Hager B. Effects of plate bending and fault strength at subduction zones on plate dynamics , J.

Mantle convection with strong subduction zones , Geophys. Lithgow-Bertelloni C. How mantle slabs drive plate tectonics , Science , , — The thermal evolution of an Earth with strong subduction zones , Geophys.

Davies G. Thermal histories of convective Earth models and constraints on radiogenic heat production in the Earth , J. Effect of plate bending on the Urey ratio and the thermal evolution of the mantle , Earth planet.

Elder J. Google Preview. Elsasser W. Runcorn S. Forsyth D. Uyeda S. On the relative importance of the driving forces of plate motion , Geophys.

Golitsyn G. Simple theoretical and experimental study of convection with some geophysical applications and analogies , J. Labrosse S. Tackley P. Convective heat transfer as a function of wavelength: implications for the cooling of the Earth , J. Oceanic plate motions driven by lithospheric thickening and subducted slabs , Nature , , — A simple global model of plate dynamics and mantle convection , J. Hirth G.

Kohlstedt D. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere , Earth planet. Long wavelength convection, Poiseuille-Couette flow in the low-viscosity asthenosphere and the strength of plate margins , Geophys. Howard L. Convection at high Rayleigh number ,in Proceedings of the 11th Int. Munich , Springer , Berlin. Korenaga J. Energetics of mantle convection and the fate of fossil heat , Geophys.

How does small-scale convection manifest in surface heat flux? Depth-dependent rheology and the horizontal length scale of mantle convection , J.

McKenzie D. Weiss N. Speculations on the thermal and tectonic history of the Earth , Geophys. Roberts J. Some remarks on heat flow and gravity anomalies , J. Mitrovica J. Forte A. A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data , Earth planet.

Moresi L. Gurnis M. Constraints on the lateral strength of slabs from three-dimensional dynamic flow models , Earth planet. Morris S. Viscosity stratification and the horizontal scale of end-driven cellular flow , Phys.

Fluids , 20 , Dziewonski A. Boschi E. Schubert G. Turcotte D. Olson P. Stevenson D. Cassen P. Whole planet cooling and the radiogenic heat-source contents of the Earth and moon , J. Sharpe H. Peltier W. A thermal history model for the Earth with parameterized convection , Geophys. Silver P. Behn M. Intermittent plate tectonics?

Sleep N. Evolution of the mode of convection within terrestrial planets , J. Torrance K. Thermal convection with large viscosity variations , J. Oxburgh E. Finite amplitude convective cells and continental drift , J.

Geodynamics , Cambridge University Press , Cambridge. Valencia D. Convection scaling and subduction on Earth and super-Earths , Earth planet. Heuret A. Lallemand S. Reconciling strong slab pull and weak plate bending: the plate motion constraint on the strength of mantle slabs , Earth planet. Zhong S. Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature , J.

The purpose of this section is to evaluate the boundary term in eq. This section will describe how the energy balance equations are solved in practice. Flow and dissipation in the mantle may be calculated for an arbitrary number of layers. To keep things simple we will solve the equations for a plate above a two layer system with an upper mantle or an asthenosphere and a lower mantle.

Eqs A7 , A8 , A9 , A10 , A12 , A13 and A14 provide seven equations to solve for the seven unknown variables , , , , , A and b in terms of the model parameters and the plate velocity. With , , , , A and b now known, eq. A4 may be used to calculate the velocity in any of the layers. Oxford University Press is a department of the University of Oxford.

It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide. Sign In or Create an Account. Sign In. Advanced Search. Search Menu. Article Navigation. Close mobile search navigation Article Navigation. Volume Article Contents Summary. An analytic model of convection in a system with layered viscosity and plates. Crowley , John W. E-mail: jcrowley fas. Oxford Academic.

Richard J. Cite Cite John W. Select Format Select format. Permissions Icon Permissions. Summary We present an analytic boundary layer model for thermal convection with a finite-strength plate and depth-dependent viscosity.

Planetary tectonics , Dynamics of lithosphere and mantle , Heat generation and transport , Mechanics, theory, and modelling. The first model was a classic boundary layer model with a conductive thermal boundary layer at the surface overlying a deeper convective layer that was essentially isothermal.

The horizontal speed of the boundary layer matched the horizontal speed of the top of the convecting layer. The boundary layer cooled by conduction as it moved from ridge to trench on the surface, and the sinking boundary layer provided the energy to balance viscous dissipation in the convecting region.

The average surface heat flow, per unit length in the third dimension, for a plate of length L with constant uniform speed is. The boundary layer model gives a dimensionless heat flow the Nusselt number and dimensionless plate velocity of for fixed aspect ratio. This model was based on the growth of a static conductive boundary layer on the heated surface that became unstable to convection when its thickness d grew to the critical value for convection to occur, that is,.

Open in new tab Download slide. We can separate dissipation in the lithosphere from dissipation in the asthenosphere and mantle. The third term in eq. The body force will be directed vertically, and with the vertical velocity component v , the rate of change of potential energy in the system is then given by.

Only the heat flow is required and for several simple cases the calculation of is straightforward. For a high Rayleigh number system at steady state and with uniform internal heating, eq. Using eqs 9 , 13 and 18 , the integral form of the mechanical energy equation for the convective cell, given by eq.

The dissipation is. At subduction zones, an infinitesimally small block is sheared from the plate into the mantle. The two dissipation models for the lithosphere represent end member models for how the plate would behave if it deformed either entirely viscously eq. The behaviour of a real plate is expected to lie somewhere between these two end members. If the thickness of the plate is defined by its thermal thickness eq.

The dissipation for the plastically deforming plate scales as and the models scale the same way with the plate velocity. In general, both models give the same result for a fixed aspect ratio when.

Thus, either deformation model captures the physics of both deformation mechanisms. We can define an effective stress for the plate equal to. Dissipation occurs at the interface between the subducting and overriding plates. This dissipation may be parametrized using an effective shear stress on the fault zone,. We assume that the surface area of the interface will be proportional to the thickness of the overriding plate other criteria could be used.

We consider first the interior horizontal flow. Flow in the interior of the cell is driven by the motion of the plate on the surface and by pressure gradients across the cell that arise from density gradients. The equation to be solved is the momentum equation given by eq. The horizontal velocity is then governed by. We approximate the pressure gradient in the mantle with a simple linear dependence on the depth. Substituting eq. Using constraints 29 - 32 provides the depth-dependent horizontal velocity in the convective cell as a function of the viscosity structure, plate velocity , and the basal shear stress.

As such, A and b are completely determined by the boundary conditions and the mass conservation constraint. They are not free parameters in the model. The analytic solution for U y is a third-order polynomial in depth y. The coefficients of U y are a complicated function of the layer thicknesses and viscosities and do not provide any physical insight see Appendix for details.

Thus, we will continue to express the horizontal velocity field as U y to avoid unnecessarily large equations. Then and the dissipation from lateral flow in the asthenosphere and mantle, using eq.

Dissipation in the corner flow of the cell is difficult to calculate without a more complex flow model that accounts for the vertical flow in the ends of the cell. Peltier, The kinematics and dynamics of poloidal-toroidal coupling in the mantle flow: The importance of surface plates and lateral viscosity variations, Adv. Woodward, Global 3D mantle structure and vertical mass and heat transfer across the mantle from joint inversions of seismic and geodynamic data, J.

Woodward, and A. Dziewonski, Joint inversions of seismic and geodynamic data for models of three-dimensional mantle heterogeneity, J. Dziewonski, and R. Earth Planet. Goldberg, D. Hager, B. Clayton, Constraints on the structure of mantle convection using seismic observations, flow models and the geoid, in Mantle Convection , edited by W.

Peltier, pp. Jordan, T. Karato, S. B, 61—66, Kido, M. Cadek, Inferences of viscosity from the oceanic geoid: Indication of a low viscosity zone below the km discontinuity, Earth Planet.

Yuen, O. Cadek, and T. Nakakuki, Mantle viscosity derived by genetic algorithm using oceanic geoid and seismic tomography for whole-mantle versus blocked-flow situations, Phys.

Google Scholar. King, S. Masters, An inversion for radial viscosity structure using seismic tomography, Geophys. Li, X. Romanowicz, Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, J. Mitrovica, J. Forte, Radial profile of mantle viscosity: Results from the joint inversion of convection and postglacial rebound observables, J. Panasyuk, S. Hager, and A. Forte, Understanding the effects of mantle compressibility on geoid kernels, Geophys.

Ribe, N. Ricard, Y. Froidevaux, and L. Fleitout, Global plate motion and the geoid: a physical model, Geophys. Vigny, and C. Doglioni, and R. Sabadini, Differential rotation between litho-sphere and mantle: A consequence of lateral mantle viscosity variation, J. Richards, C. Lithgow-Bertelloni, and Y. Le Stunff, A geodynamic model of mantle density heterogeneity, J. Nataf, and J. Richards, M. Hager, Geoid anomalies in a dynamic earth, J. Hager, Effects of lateral viscosity variation on long-wavelength geoid anomalies and topography, J.

Sen, M. Stoffa, Rapid sampling of model space using genetic algorithms: examples from seismic waveform inversion, Geophys. Tarantola, A. Valette, Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Space Phys. Thoraval, C. Machetel, and A. Cazenave, Influence of mantle compressibility and ocean warping on dynamical models of the geoid, Geophys. Wen, L. Anderson, Layered mantle convection: A model for geoid and topography, Earth Planet.

Zhang, S. Christensen, The effect of lateral viscosity variations on geoid, topography and plate motions induced by density anomalies in the mantle, Geophys. Download references. You can also search for this author in PubMed Google Scholar. Correspondence to Motoyuki Kido. Reprints and Permissions. Synthetic tests of geoid-viscosity inversion: A layered viscosity case. Earth Planet Sp 50, — Download citation.

Received : 02 April Revised : 02 November



0コメント

  • 1000 / 1000