*This is a mammoth post which is a product of a week of obsessing over a question that came to me from nowhere. I’ve written it for a general audience but the technical details are in red – they are ‘extra reading’ for anyone wanting to know a little more.*

# Introduction

One of the scary thoughts that keeps me awake at night is the ‘Kessler Syndrome’. Any readers who have seen the film ‘Gravity’ are already familiar with it. The concept was first proposed by Donald J Kessler back in the 1970’s and goes like this:

- As time goes along we will launch more and more satellites into orbit.
- As the number of satellites increases, it is inevitable that they will eventually start colliding.
- A satellite collision causes both satellites to fragment into thousands of pieces which remain in orbit (this is an understatement – see below).
- As fragments build up, the chance of these fragments colliding with other satellites becomes higher and higher.
- Eventually the resulting debris will build up so much that we will be unable to launch new satellites to orbit.

So, it’s a classic case of a ‘cascade failure’, where a single failure is amplified repeatedly, eventually leading to the failure of the entire system.

Several events in recent years have brought this scenario into the public eye:

- The 2007 test of a Chinese anti satellite missile, which caused a single satellite to fragment into about 2000 trackable fragments.
- The subsequent US deployment of a similar system (which was actually of no consequence in this context, as explained later).
- The 2009 collision between an iridium and Kosmos satellite.

Is it possible that events like these will continue and eventually prevent us from reaching space ever again?

In short, no. The Kessler syndrome hypothesis is about a build up of debris *in orbit.* Even in a worst case scenario, we would still be able to leave Earth orbit and ‘explore the universe’.

## So we can’t have astronauts floating around in orbit, what’s the big deal?

Turns out, losing the capability to deploy and maintain satellites in LEO (Low Earth Orbit) is a huge deal. Let’s look at just one example – Satellite Navigation.

Imagine if the world lost its Satellite Navigation capability. Doesn’t sound too bad does it? You woldn’t be able to use your maps app on your phone, what’s the big deal? Off the top of my head:

- Pretty much all commercial and military airlines rely on GPS. Sure, they have alternative systems (which are vastly inferior). But the savings introduced by using GPS no doubt add up to billions.
- Since we started using it, emergency and disaster response services have saved thousand of lives which they wouldn’t otherwise without satellite navigation capability.
- Science/research etc. So many studies rely on reliable and accurate positional tracking.
- Self driving cars – they will die a death with no sat nav.
- GPS time signal. The GPS (and GLONASS) systems rely on atomic clocks which broadcast an incredibly precise time signal. I don’t even know where to begin with the uses for this – one example is power generation and distribution. If we lost GPS, we’d also lose the ability to effectively distribute electricity until some backup option could be implemented.

When I talk about huge costs, it’s not just about some folk not being as rich as before. When I talk about additional money being required for aircraft or self driving cars, that is money that otherwise may have been used for building schools or hospitals. So the loss of the satellite navigation network would be a real disaster; and that is just one of the many things you rely on satellites for every day.

# Modelling Satellite Collisions

It is easy to start worrying about the Kessler Syndrome happening, especially with some of the literature available online. But is it really that likely? I decided to try and come up with a model to determine exactly that. Here was my initial train of thought:

## A deterministic model

It is clearly an exercise in probability so one way to predict would be to determine the factors which make a collision possible. Some common sense and reading around, and I decided the main factors which influence satellite collisions are:

- Satellite altitude. If there are more satellites at a given altitude they are more likely to collide with each other.
- Orbital inclination. This is a measure of how inclined the orbit is. A satellite orbiting the equator has zero inclination. One orbiting around the poles has 90 degrees of inclination.
- Satellite surface area. Obviously a larger satellite is more likely to be involved in a collision.
- Atmospheric drag. Orbiting objects gradually decay and fall back to Earth, this serves to clear objects from orbit.
- The number of satellites launched per year.
- The number objects which spontaneously explode per year – turns out objects in orbit explode surprisingly regularly which adds more debris.

To create a model based on the above, we would have to devise a time dependent function for each and then integrate over the time we are interested in. The output for each timestep would become the input for the next. This would be a nightmare.

## A probabilistic model

Make a simulation. We could take the orbital parameters of satellites currently in orbit and model the orbit of each one over time. New objects could be added to simulate launches/explosions/collisions etc. A collision could only occur if 2 orbits intersect. In that case, we would determine the probability both objects were are the intersection point and simulate them colliding in only a select few cases corresponding to that probability.

This approach is incredibly nightmarish. Simulating the orbits for each object is fine on a modern computer. But there are about 13000 trackable objects in Low Earth Orbit right now. We would have to compare each orbit with all others to determine whether they intersect.

This may still be possible – there are some shortcuts which we can take to exclude certain orbits. But then, the calculation to determine whether an intersection occurs is actually surprisingly difficult. My back of envelope calculations put the time to simulate 100 years of this is on the order of years!

But, it has actually been done before. Nasa have actually run numerour runs for the above simulation for the next 100 years on their supercomputers – amazing!

My quad core AMD A8 is going to need something a bit more simple.

## A probabilistic model which doesn’t need a supercomputer

We just go back to the first approach and let the computer simulate the product of all the functions for each timestep. That sentence doesn’t really make much sense by itself, so let’s look at what my model does step by step:

We start with a table of every object in orbit. This table includes the object’s mass, cross sectional area and orbital parameters

For day 1:

- From all the altitudes, determine the ‘density’ of satellites at each altitude.
- From all the inclinations, determine the ‘density of satellites at each inclination.
- Calculate the solar activity for the day – this is time dependent, so only uses the day number as an input.
- Determine if a launch occurs on that day. About 1 satellite is launched every 5 days at present.
- Determine if an object in orbit explodes. If so, fragment a random object.

Then for each object in orbit:

- Use the input parameters (available online) to calculate the altitude of each object in orbit.
- Calculate the loss of altitude for the day – this depends on the solar flux, the current altitude, the mass and the cross sectional area.
- Take the density of objects at the current altitude, cross sectional area, density at the given inclination, and a calibration factor. Multiply these together to get the probability of the given object colliding with another on that day. Call this probability P.
- Generate a random number between 1 and (1/P).
- If that number is equal to 1 – a collision occurs – fragment the satellite into pieces, along with another mass to be chosen at random.

We then have a new set of objects and respective orbits. Repeat all the above for the next day. Continue repeating for 36500 days (approx 100 years).

There is tons to go through here so the next section is a bit technical.

# How to determine the functions and parameters

Let’s go through each aspect of the simulation:

## The raw data

The initial conditions of the simulation were the current orbital parameters of all trackable objects in obit. These are tracked with NORAD radar and the information is freely available online in the form of TLE (two line element) data. The data is current on 19/02/16.

I decided to only include objects below 2000km as this is considered low Earth orbit and is, by far, the most crowded ‘band’ of Earth orbit.

NORAD is capable of tracking objects down to about 10cm in diameter. Luckily, this seems to be around the minimum size capable of causing a violent fragmentation due to collision, so we will only concern ourselves with objects above this size.

## Satellite Altitudes

The TLE data includes the inclination and period of each satellite. The periapsis (lowest point of orbit) can be calculated from these (see wikipedia page on orbital calculations). At present my model uses the periapsis as an ‘effective altitude’ as I found it difficult to model effects of eccentricity, for both collision risk and orbital decay. This probably does not have too much effect as most LEO orbits have a low inclination.

## Density of objects at each altitude

This was easily determined using the ‘density’ function in R and fitting a function using ‘Approxfun’. For each altitude, we could then determine the object density. I’ll call this value f(d)

Here is the density of objects by altitude at the time of writing. This was the ‘initial state of the model:

## Change in periapsis over a day

Satellite orbits decay depending on the mass,XSA, current periapsis and current solar activity. The function for change in orbital period over a day was taken from the paper in the sources, it is:

**dP/dt = -3πaρ (Ae/m) **

**Where:**

**dP/dt = rate of change of orbital period**

**a = semimajor axis (calculated from TLE data)**

**ρ = atmospheric density (see orbital decay paper in references for calculation method)**

**Ae = effective cross sectional area. Assumed to be 2 X XSA.**

**m = object mass**

This serves to remove objects from our ‘pool’. However, the formula is only valid up t oaround 1200km. I have assumed objects above this altitude do not decay for now. I may go back and improve this model in the future.

## Inclination enhancement factor

The number of satellites at different inclinations varies. The majority of satellites are in near-polar orbits. Satellites sharing inclinations spend a longer time intersecting orbits than those which differ in inclination.

Coming up with a function to model this was difficult. However, someone over at NASA has already done this. So I used the following enhancement factors:

## Solar activity

A bit of research led me to conclude that solar activity follows an 11 year cycle and varies between about 70 and 150 sfu in the 10.7cm band which is apparently the standard indicator for this.

## Satellite mass

All the objects were assigned a random mass from a Weibull distribution. Some research indicated that the total mass of objects in orbit above the 10cm size is around 550000kg, so I used a size parameter of 55000 and a shape factor of 1. This gives a total of around 5500T of objects.

I chose a Weibull distribution based on my ‘engineering instinct’. This is generally the distribution you’d expect to see for different sized objects, where less larger ones and more smaller ones are expected.

A check of the maximum mass from this distribution returned a value of 8600kg which is consistent with large satellites such as Envisat. (The ISS is far larger but that has debris evasion capability so is not part of the simulation).

## Satellite cross sectional area (XSA)

A little reading around gave me values of mass and cross sectional area for various satellites. I decided to only include the surface of the main body as a collision with a solar panel is unlikely to produce a full fragmentation (as evidenced by a Cubesat collision in 2013).

We’d expect the cube root of the mass squared to be proportional to the surface area (ie, the square cube law). The following relationship held up remarkably well, so I used it to determine the XSA of all objects in orbit:

**A=√(m^(1/3) )/16**

where:

A = cross sectional area of object (excluding solar panels)

m=object mass

A 10cm diameter object would have a mass of 0.044kg in this case, so that is the smallest mass we will model. All masses below 0.044kg are automatically removed from the simulation. These objects are generally incapable of causing meaningful fragmentations as discussed in the ‘modelling a collision’ section.

## Does a launch occur?

The first means by which objects are added to orbit is by launches. Looking at historical data, the chance of a satellite being launched on a given day is 0.222. So, the simulation generates a random number between 1 and 1000. If this number is less than 223, a new object is added to the simulation. The size of this object is randomly chosen from the Weibull distribution which determined the masses. The periapsis is chosen randomly from the initial distribution of periapses.

## Does an object explode on a given day?

Historical data indicates that an average of about 5 objects in orbit explode every year. This may be due to residual fuel in discarded stages, undocumented collisions etc. This translates to a rate of 0.015 per day. So, R generates a random number between 1 and 1000 every ‘day’ of the simulation. If this number is less than or equal to 15, a random object is assigned to explode.

The explosion is modelled as a fragmentation into 50 objects, following a Weibull distribution with a shape factor of 1. The inclinations and orbital period of the new objects are assumed to vary from the original inclination and orbital period by a normal distribution with standard deviations of 0.1 degrees and 0.1 seconds respectively.

These standard deviations are based on eyeballing data from previous explosions (see the NASA paper in the links). I may refine these values at a later date.

## Does a collision occur?

As stated in the summary above, each satellite has a ‘collision flux’ value – that is the probability of it colliding with another object per dat. This is calculated as follows:

**For an object with:**

**Periapsis = Pe**

**Inclination = i**

**cross sectional area =A**

**flux callibration factor = k**

**Collision flux [/day] = density of objects at Pe X enhancement factor for i X A X k**

The value of k was determined empirically. Historically only one satellite collision resulting in a full scale fragmentation has occured. But it seems that one collision every 10 years is a reasonable starting point. So I summed the collision flux for all objects for one day and adjusted k until the total flux represented one collision every 10 years.

For every object, for each day of the simulation, a random number was generated between 1 and (1/collision flux). If this number was equal to 1, the object was determined to collide. The collision mas was chosen at random from the existing pool of object masses in orbit. Likewise for the inclination of the ‘collider’.

# Modelling a collision

A collision is much (much much) more energetic than an explosion. For reference, one 1000kg object at orbital velocity has 30GJ of kinetic energy. An explosion of, say 1kg of hydrazine propellant (the kind of thing that would cause an in orbit explosion) is 1.6MJ. This is why collisions are so much worse – they completely fragment the bodies.

## Mass of the fragmented bodies

The 2007 Chinese ASAT test fragmented a satellite of mass approx 2000kg into about 2000 trackable pieces. Likewise, the iridium satellite which collided in 2009 had a mass of around 1000kg and fragmented into about 1000 trackable pieces.

So, empirically it’s fair to say a collision results in M trackable fragments for an object of mass M.

Again, we’ll assume a Weibull distribution with shape factor 1.

For objects below 50kg, I assumed 50 trackable pieces were created, as around this size, many of the fragments become smaller than the trackable size. This isn’t necesserily the strongest assumption but it isn’t of much consequence since either way, a small object fragmenting does not produce many trackable objects compared to a large one.

Incidentally, the minimum size (0.044kg as dictated by the minimum trackable diameter of 10cm) is about the minimum size of a fragment which can cause a meaningful collision. The NASA paper which discusses their deterministic simulation states that the minimum size of an object which can fragment a larger object is 0.01% of the larger objects mass. So, bodies smaller than 0.044kg will not fragment those above 50kg. This is another reason for the model changing at the 50kg limit and for the minimum size of 10cm.

## Periapsis (orbital period) of the fragmented masses

I assumed the fragmented masses are distributed normally around the orbit of the original body with a standard distribution of 0.1s. This is a bit of a guess but seems to match up visually pretty well with the 2009 and 2007 collisions. I may revisit this assumption in the future.

This gives a significant ‘step up’ in the density of objects at the given periapsis following an explosion. Such step ups are clearly visible in historical NORAD data immediately after the 2007 ASAT test and the 2009 collision.

## Inclination of orbital masses

I made a similar assumption to that of the periapsis. I assumed the fragmented masses have inclinations distributed normally to the original body with a standard deviation of 1 degree.

**All the above are calculated for each object in orbit (where applicable) for each timestep in the simulation (1 timestep = 1 day). Following each timestep, the new orbital periods of each object are calculated and the above process is repeated for 365000 days.**

# The Results

## Simulation time

I ran the above simulation 10 times, using a different seed number for random number generation in each case (basically a Monte Carlo simulation). Ideally I’d run it far more times but it took 11 hours and I need my laptop for everyday use too. I may eventually re-run it on my upcoming Raspberry Pi cluster (see future post!). I may also re-write it in Python. This would be a little more difficult to code but would run faster…probably.

This is a rare case for my hobby projects – I have actually tried to make my code efficient as possible and commented throughout. Anyone wishing to download it may do so at the following link:

https://www.dropbox.com/s/h1k2t7xqd8czwke/final%20working%20simulation.R?dl=0

And here is the TLE data:

https://www.dropbox.com/s/fumi6z2x3n0lfuz/tle%20text.txt?dl=0

## First Observations

I think there is some good news and some bad news in the results.

Reading anything in the media related to the Kessler Syndrome gives the impression that one day everything is going to go completely haywire and every object in orbit is going to violently collide and spray debris everywhere (like in Gravity). This does not happen in any of the 10 cases.

However, let’s remember these cases assume the future launch rate for satellites does not increase (not a strong assumption based on rapidly reducing launch costs at the time of writing). Even with this assumption, the amount of debris increases over time significantly. More worryingly, the most commonly used altitude band 700-900km looks really quite unstable.

## Variation of number of objects at given altitudes with time

For the ten simulation runs, here is the evolution of the number of objects in orbit over time:

Here is a video of the density of objects at each altitude, changing over time (ie for the whole 100 year simulation run – the day is given at the top of the histogram). This video is for the case which generated the most debris, however, it is not actually *much* worse than any of the other cases as can be seen in the plot below:

There are all sorts of things to note here:

- Collisions make a complete mess. At first everything looks pretty stable – a few explosions can be seen but their effect isn’t really too major. When we see a sudden jump in the amount of debris at a given periapsis, that’s a collision.
- Anything below 400km in altitude decays
*quickly.*Just look at the dropoff. - A single unfortuntely placed collision can make a given band of altitudes really unstable. Look how quickly the number of objects at the most common altitudes grows towards the end.
- Even though my model doesn’t account for orbital decay above 1200km, there doesn’t look like much risk of orbits this high becoming unstable due to debris – there just isn’t enough there to pose a major risk.

## Effect of collider mass on evolution of the model

Let’s look at the masses of objects which collided in this run of the simulation:

Now let’s look at the masses of objects which collided in the simulation run which produced the least number of overall objects at the end of the simulation period:

So, it looks like the number of collisions was similar but the masses were quite different. This explains the difference in the number of objects in orbit at the end of the two runs. It also shows that a few large objects colliding can quickly cause a whole band of orbits to become unstable.

A real life concern here would be the now disused satellite Envisat. This weight 8000kg and orbits at about 800km. The largest collider object in my simulation was 4000kg and look at the debris cloud it created. Now imagine if Envisat was involved in a similar collision today. Based on my model, it would increase the number of trackable objects in orbit by about another 50%.

Another ‘what if’ scenario would be the ISS being involved in a fragmenting collision. The ISS weighs about 400000kg. Even a partial fragmentation would be a complete and utter mess on paper. In reality it may not be so bad since it orbits at just 400km, so the fragments would quickly decay. Fortunately this isn’t really a major concern anyway – the ISS has thoroughly tried and tested debris evasion capability.

## Probability of a collision over the run life

Here is the total probability of any satellite colliding on a given day for the 10 simulation runs:

So, the probability of an impact increases significantly over the run life. It is not surprising that it almost perfectly mirrors the number of objects in orbit (see plot above). At the start of the simulation runs, the probability of an impact is 0.022% on a given day. By the end of the simulation period, the probability increases to between 0.07% and 0.19%. This is a significant increase but is still a low figure overall and would probably fall within the acceptable risk for unmanned satellites.

And just for confirmation, here is the solar flux over the simulation period. It is somewhat simplified here – just following a sin wave with an 11 year period. In reality, we’d expect peaks up to a maximum value of 400.

# Conclusion

Although it is inevitable that collisions between orbiting objects will accelerate over the next few decades, my model does not suggest the rate will be enough the prohibit the use of LEO. This does, however depend on a few assumptions:

- The launch rate does not significantly increase.
- No more government idiots decide to destroy satellites in orbit deliberately. Reading this post makes it pretty evident that the 2007 test was probably signed off by some idiotic bureacrat who did not understand the ramifications of his/her (probably ‘his’ – only men get involved in stupid ‘show of force’ decisions like this) actions.
- Eccentricity effects are ignored.
- An exceptionally large object in a crowded orbit does not collide in the near future
- Exceptional solar activity does not significantly increase the collision rate.

So ~~un~~fortunately, the happenings of *Gravity *probably won’t be going on anytime soon. However, there is a significant chance we will see a lot more losses of active satellites in the most densely populated orbits. In a worst case scenario, this could certainly lead to real knock on effects down here on Earth.

As much as I don’t condone silly ‘everyone look how scary our military technology is’ events by governments, we can also see that the American ASAT test in 2008 had effectivly zero consequence in space-debris terms. The satellite in question was due to de-orbit naturally in the next few days. It would have been at the extreme end of the plots in the above video, so the debris created would have naturally fallen to Earth within hours. This did not stop the world’s media completley misreporting the event and comparing it to the Chinese test (the result of which is still the largest single contributor to orbital debris ever).

Besides, being able to shoot down a car sized object travelling at 10 times the speed of a rifle bullet on the first attempt is pretty damn impressive.

**Spot any mistakes? Thoughts? Disagree with any of the above? Just let me know in the comments.**

# References and further reading

http://webpages.charter.net/dkessler/files/Collision%20Frequency.pdf The original paper by Donald Kessler.

http://www.sws.bom.gov.au/Category/Educational/Space%20Weather/Space%20Weather%20Effects/SatelliteOrbitalDecayCalculations.pdf Orbital decay calculation

http://www.amostech.com/TechnicalPapers/2012/Orbital_Debris/NIKOLAEV.pdf The probabilistic model of orbital debris discussed at the start of this article. This required significant time to run on a supercomputer network. I am frankly amazed my model gave similar results!

http://orbitaldebris.jsc.nasa.gov/library/SatelliteFragHistory/TM-2008-214779.pdf NASA data on on-orbit fragmentations. I based the explosion modelling and aspects of the collision model on this.

http://www.duncansteel.com/archives/1475 A similar study with some extra detail included

http://satellitedebris.net/Database/LaunchHistoryView.php?hname=LaunchHistoryYearViewDetailEdit1LaunchHistoryView_handler&fk0=US Source of launches per year

Pingback: How to build a supercomputer – part 1 | brain -> blueprint -> build