diff --git a/models/bubblecube.obj b/models/bubblecube.obj
new file mode 100644
index 0000000000000000000000000000000000000000..3de663a1c664bcfe796dfce377b5aaaa562c56f5
--- /dev/null
+++ b/models/bubblecube.obj
@@ -0,0 +1,38 @@
+# Blender 4.3.2
+# www.blender.org
+o Cube.001
+v 1.486580 0.021086 1.427890
+v 1.486580 0.687774 1.427890
+v 1.486580 0.021086 -1.427890
+v 1.486580 0.687774 -1.427890
+v 3.486580 0.021086 1.427890
+v 3.486580 0.687774 1.427890
+v 3.486580 0.021086 -1.427890
+v 3.486580 0.687774 -1.427890
+vn -1.0000 -0.0000 -0.0000
+vn -0.0000 -0.0000 -1.0000
+vn 1.0000 -0.0000 -0.0000
+vn -0.0000 -0.0000 1.0000
+vn -0.0000 -1.0000 -0.0000
+vn -0.0000 1.0000 -0.0000
+vt 0.375000 0.000000
+vt 0.625000 0.000000
+vt 0.625000 0.250000
+vt 0.375000 0.250000
+vt 0.625000 0.500000
+vt 0.375000 0.500000
+vt 0.625000 0.750000
+vt 0.375000 0.750000
+vt 0.625000 1.000000
+vt 0.375000 1.000000
+vt 0.125000 0.500000
+vt 0.125000 0.750000
+vt 0.875000 0.500000
+vt 0.875000 0.750000
+s 0
+f 1/1/1 2/2/1 4/3/1 3/4/1
+f 3/4/2 4/3/2 8/5/2 7/6/2
+f 7/6/3 8/5/3 6/7/3 5/8/3
+f 5/8/4 6/7/4 2/9/4 1/10/4
+f 3/11/5 7/6/5 5/8/5 1/12/5
+f 8/5/6 4/13/6 2/14/6 6/7/6
diff --git a/shaders/bubble.frag b/shaders/bubble.frag
new file mode 100644
index 0000000000000000000000000000000000000000..5b6b505a12643dac903b56cc7dbe77f3822ce998
--- /dev/null
+++ b/shaders/bubble.frag
@@ -0,0 +1,121 @@
+#version 150
+
+const float step_size = 0.1;
+const int NUMBER_OF_STEPS = 48;
+const int NUMBER_OF_BALL_STEPS = 16;
+const float START_STEP = 0.1/30.0; // Half the size of the smallest ball radius
+const float MINIMUM_HIT_DISTANCE = 0.002;
+const float MAXIMUM_TRACE_DISTANCE = 100.0;
+
+// Constants related to reflection/refraction
+const float R0 = 0.04;
+const float WATER_REFRACTION = 0.75;
+const float PI = 3.1415926535897f;
+const int REFRACTION_ITER = 3;
+
+in vec3 world_pos;
+
+uniform sampler2D sky;
+uniform int num_balls;
+uniform vec3 camera_pos;
+
+out vec4 out_Color;
+
+struct Ball {
+  vec4 pos;
+  float radius;
+};
+
+layout(std140) uniform BubbleBuffer {
+  Ball balls[350];  // TODO: Need to manually update size
+};
+
+float distance_from_sphere(vec3 pos, Ball ball) {
+  return length(pos - ball.pos.xyz) - ball.radius;
+}
+
+// Function taken from: https://iquilezles.org/articles/distfunctions/
+float opSmoothUnion(float d1, float d2, float k) {
+  float h = clamp(0.5 + 0.5*(d2-d1)/k, 0.0, 1.0);
+  return mix(d2, d1, h) - k*h*(1.0-h);
+}
+
+float SDF(vec3 position) {
+  float min_dist = distance_from_sphere(position, balls[0]);
+
+  // TODO: iterate over uniform num_balls
+  for (int i = 0; i < 200; ++i) {
+    float dist = distance_from_sphere(position, balls[i]);
+
+    min_dist = opSmoothUnion(min_dist, dist, 0.08);
+  }
+
+  return min_dist;
+}
+
+vec3 approximate_normal(vec3 position) {
+  const vec3 small_step = vec3(0.0001, 0, 0);
+  return normalize(
+    vec3(
+      SDF(position + small_step.xyy) - SDF(position - small_step.xyy),
+      SDF(position + small_step.yxy) - SDF(position - small_step.yxy),
+      SDF(position + small_step.yyx) - SDF(position - small_step.yyx)
+    )
+  );
+}
+
+float fresnel(vec3 normal, vec3 view) {
+  float cos_theta = dot(normal, -view);
+  float x = 1.0 - cos_theta;
+  return mix(R0, 1.0, x * x * x * x * x);
+}
+
+vec2 sphere_uv(vec3 direction)
+{
+  float u = 0.5 + atan(direction.z, direction.x) / (2 * PI);
+  float v = 0.5 + asin(direction.y) / PI;
+
+  return vec2(u, v);
+}
+
+vec3 ray_march(vec3 ro, vec3 rd) {
+  float total_distance_traveled = 0.0;
+
+  for (int i = 0; i < NUMBER_OF_STEPS; ++i) {
+    vec3 current_position = ro + total_distance_traveled * rd;
+
+    float min_dist = SDF(current_position);
+
+    // hit
+    if (min_dist < MINIMUM_HIT_DISTANCE) {
+      const vec3 color = vec3(0.2, 0.3, 0.8);
+
+      const vec3 light_position = vec3(0, 10, 0);
+      vec3 normal = normalize(approximate_normal(current_position));
+      vec3 direction_to_light = normalize(light_position - current_position);
+      float diffuse_intensity = max(0.0, dot(normal, direction_to_light));
+      
+      if (diffuse_intensity < 0.95) {
+        discard;
+      }
+
+      //return color + diffuse_intensity * vec3(1, 1, 1) * 0.5;
+      return vec3(diffuse_intensity, diffuse_intensity, diffuse_intensity);
+    }
+
+    if (total_distance_traveled > MAXIMUM_TRACE_DISTANCE) {
+      discard;
+    }
+
+    total_distance_traveled += min_dist;
+  }
+  return vec3(0.9, 0.9, 0.9);
+}
+
+
+void main(void)
+{
+  vec3 dir = normalize(world_pos - camera_pos);
+  vec3 shaded_color = ray_march(world_pos, dir);
+	out_Color = vec4(shaded_color, 1.0);
+}
diff --git a/shaders/bubble.vert b/shaders/bubble.vert
new file mode 100644
index 0000000000000000000000000000000000000000..655b91d54b661f3108c5595e0a7267a7ceb835f8
--- /dev/null
+++ b/shaders/bubble.vert
@@ -0,0 +1,21 @@
+#version 150
+
+in  vec3  in_Position;
+in  vec3  in_Normal;
+in  vec2  in_TexCoord;
+
+uniform mat4 projectionMatrix;
+uniform mat4 modelToWorldToView;
+uniform mat4 invModelToWorld;
+uniform int num_balls;
+uniform float screenWidth;
+uniform float screenHeight;
+uniform vec3 camera_pos;
+
+out vec3 world_pos;
+
+void main(void)
+{
+	gl_Position = projectionMatrix * modelToWorldToView * vec4(in_Position, 1.0);
+	world_pos = in_Position;
+}
diff --git a/shaders/waterfall.frag b/shaders/waterfall.frag
index 03b855d20ab53e3e354eaa13b759196c5a46cbc8..292d11f2b73c461576f0321363f3d1423acfd374 100644
--- a/shaders/waterfall.frag
+++ b/shaders/waterfall.frag
@@ -27,7 +27,7 @@ struct Ball {
 };
 
 layout(std140) uniform BallBuffer {
-  Ball balls[100];  // TODO: Need to manually update size
+  Ball balls[250];  // TODO: Need to manually update size
 };
 
 float distance_from_sphere(vec3 pos, Ball ball) {
@@ -44,7 +44,7 @@ float SDF(vec3 position) {
   float min_dist = distance_from_sphere(position, balls[0]);
 
   // TODO: iterate over uniform num_balls
-  for (int i = 1; i < 30; ++i) {
+  for (int i = 0; i < 250; ++i) {
     float dist = distance_from_sphere(position, balls[i]);
 
     min_dist = opSmoothUnion(min_dist, dist, 0.2);
@@ -86,7 +86,6 @@ vec3 reflection(vec3 pos, vec3 dir) {
   return sky_color;
 }
 
-// TODO: implement this
 vec3 refraction(vec3 world_pos, vec3 dir) {
 
   float total_distance_traveled = START_STEP;
diff --git a/src/bubble.h b/src/bubble.h
new file mode 100644
index 0000000000000000000000000000000000000000..c8fba31db6f4132b6099788eb1c6efc806fc54ea
--- /dev/null
+++ b/src/bubble.h
@@ -0,0 +1,41 @@
+#pragma once
+
+#include "VectorUtils4.h"
+#include "object.h"
+
+struct Bubble : Object {
+
+  Bubble() = default;
+  Bubble(Model *model, GLuint program);
+  void move_bubble();
+  void gen_bubble();
+  void draw() const;
+
+  static GLuint const BINDING_POINT{2};
+  char const *UNIFORM_BLOCK = "BubbleBuffer";
+  static int const NUM_BALLS{350};
+
+  // Dimensions of the bubblearea
+  // on the form <min, max>
+  std::pair<float, float> x;
+  std::pair<float, float> y;
+  std::pair<float, float> z;
+
+  // Properties of each ball
+  struct Ball {
+    vec4 pos;
+    float radius;
+
+    // Glsls padding to make struct a multiple of vec4 (cringe)
+    vec3 padding;
+
+    Ball(vec3 pos, float radius) : pos{pos}, radius{radius} {}
+    Ball() : pos{0, 0, 0, 0}, radius{0} {}
+  };
+
+  Ball balls[NUM_BALLS];
+  vec4 ball_velocities[NUM_BALLS];
+
+  GLuint ubo;
+  GLuint block_index;
+};
diff --git a/src/bubbles.cpp b/src/bubbles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..49722ac45e383242ef90f8e5f37f63620018eb61
--- /dev/null
+++ b/src/bubbles.cpp
@@ -0,0 +1,142 @@
+#include "bubble.h"
+#include "GL_utilities.h"
+#include "object.h"
+#include <cmath>
+#include <iostream>
+#include <ostream>
+
+Bubble::Bubble(Model *model, GLuint program) : Object{model, program} {
+
+  float min_x = std::numeric_limits<float>::max();
+  float max_x = std::numeric_limits<float>::min();
+  float min_y = std::numeric_limits<float>::max();
+  float max_y = std::numeric_limits<float>::min();
+  float min_z = std::numeric_limits<float>::max();
+  float max_z = std::numeric_limits<float>::min();
+
+  for (int i = 0; i < model->numVertices; ++i) {
+    min_x = std::min(min_x, model->vertexArray[i].x);
+    max_x = std::max(max_x, model->vertexArray[i].x);
+    min_y = std::min(min_y, model->vertexArray[i].y);
+    max_y = std::max(max_y, model->vertexArray[i].y);
+    min_z = std::min(min_z, model->vertexArray[i].z);
+    max_z = std::max(max_z, model->vertexArray[i].z);
+  }
+
+  x = {min_x, max_x};
+  y = {min_y, max_y};
+  z = {min_z, max_z};
+
+  gen_bubble();
+  use();
+
+  // Create a uniform buffer object, used to send waves.
+  glGenBuffers(2, &ubo);
+  glBindBuffer(GL_UNIFORM_BUFFER, ubo);
+
+  // Allocate memory in the buffer.
+  glBufferData(GL_UNIFORM_BUFFER, sizeof(balls), nullptr, GL_DYNAMIC_DRAW);
+
+  block_index = glGetUniformBlockIndex(program, UNIFORM_BLOCK);
+
+  if (block_index == GL_INVALID_INDEX) {
+    std::cout << "Bubbles: Block name could not be found." << std::endl;
+  }
+
+  // Bind the buffer to a binding point.
+  glBindBufferBase(GL_UNIFORM_BUFFER, BINDING_POINT, ubo);
+
+  // Link the UBO to the shader’s uniform block
+  glUniformBlockBinding(program, block_index, BINDING_POINT);
+}
+
+void Bubble::draw() const {
+
+  // Bind the buffer before updating it.
+  glBindBuffer(GL_UNIFORM_BUFFER, ubo);
+  // Assign values to the `balls` array and update the UBO data
+  glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(balls), balls);
+
+  // Upload the number of balls.
+  glUniform1i(glGetUniformLocation(program, "num_balls"), NUM_BALLS);
+
+  DrawModel(model, program, "in_Position", "in_Normal", "in_TexCoord");
+  printError("Bubble draw");
+}
+
+void Bubble::move_bubble() {
+  // Moves the balls one distance in their velocity. If they fall
+  // below 0 (end of the waterfall) their height is reset of waterfall
+  // height.
+
+  for (int i{0}; i < NUM_BALLS; ++i) {
+    balls[i].pos += ball_velocities[i];
+    ball_velocities[i].y -= 0.0001;
+
+    if (ball_velocities[i].y < -0.01) {
+      ball_velocities[i].y = 0.01;
+    }
+
+    if (balls[i].pos.y > y.second - balls[i].radius) {
+      balls[i].pos.y = 0 + balls[i].radius * 3.5;
+    }
+
+    if (balls[i].pos.x < x.first) {
+      balls[i].pos.x = x.second - balls[i].radius * 3.5;
+    }
+
+    if (balls[i].pos.x > x.second) {
+      balls[i].pos.x = x.first + balls[i].radius * 3.5;
+    }
+
+    if (balls[i].pos.z < z.first) {
+      balls[i].pos.z = z.second - balls[i].radius * 3.5;
+    }
+
+    if (balls[i].pos.z > z.second) {
+      balls[i].pos.z = z.first + balls[i].radius * 3.5;
+    }
+  }
+}
+
+float lerp(float a, float b, float t);
+
+void Bubble::gen_bubble() {
+  for (int i{0}; i < NUM_BALLS; ++i) {
+
+    float radius = ((rand() / (float)RAND_MAX) * 0.3 + 1.4) / 50.0;
+
+    // pseudo random value between min-radius and max-radius
+    radius *= 3.5;
+    float pos_x =
+        lerp(x.first + radius, x.second - radius, (rand() / (float)RAND_MAX));
+    float pos_y =
+        lerp(y.first + radius, y.second - radius, (rand() / (float)RAND_MAX));
+    float pos_z =
+        lerp(z.first + radius, z.second - radius, (rand() / (float)RAND_MAX));
+    vec4 pos = vec4(pos_x, pos_y, pos_z, 0);
+    radius /= 3.5;
+
+    float x_target;
+    if (abs(x.first - pos.x) < abs(x.second - pos.x)) {
+      x_target = (pos_x - x.first)/abs(pos_x - x.first);
+    } else {
+      x_target = (pos_x - x.second)/abs(pos_x - x.first);
+    }
+    float z_target;
+    if (abs(z.first - pos.z) < abs(z.second - pos.z)) {
+      z_target = (pos_z - z.first)/abs(pos_z - z.first);
+    } else {
+      z_target = (pos_z - z.second)/abs(pos_z - z.first);
+    }
+
+    // pseudo random value between -0.001 and -0.01
+    x_target *= ((rand() / (float)RAND_MAX) * 0.5 + 0.4) / 100.0;
+    z_target *= ((rand() / (float)RAND_MAX) * 0.5 + 0.4) / 100.0;
+    float velocity_y = ((rand() / (float)RAND_MAX) * 0.5 + 0.3) / 40.0;
+    vec4 velocity = vec4(x_target, velocity_y, z_target, 0);
+
+    balls[i] = Bubble::Ball(pos, radius);
+    ball_velocities[i] = velocity;
+  }
+}
diff --git a/src/main.cpp b/src/main.cpp
index e75b340b02539cf092048730d96cebcbf5e18423..b934e443d51ae7b10ff4eacb3693ef61a911dfba 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,7 @@
 // Revised 2019 with a bit better variable names.
 // Experimental C++ version 2022. Almost no code changes.
 
+#include "bubble.h"
 #include "object.h"
 #include "surface.h"
 #include "waterfall.h"
@@ -22,6 +23,7 @@ struct Scene {
   Object ground;
   Surface surface;
   Waterfall waterfall;
+  Bubble bubble;
   Object skybox;
   vec3 sun = vec3(0, 100, 0);
 
@@ -43,6 +45,7 @@ struct Scene {
     Model *ground_model = LoadModel("models/ground.obj");
     Model *surface_model = LoadModel("models/surface.obj");
     Model *waterfall_model = LoadModel("models/wide_cube.obj");
+    Model *bubble_model = LoadModel("models/bubblecube.obj");
     Model *skybox_model = LoadModel("models/skybox.obj");
 
     // This is not a ground logic program
@@ -52,12 +55,15 @@ struct Scene {
         loadShaders("shaders/surface.vert", "shaders/surface.frag");
     GLuint waterfall_program =
         loadShaders("shaders/waterfall.vert", "shaders/waterfall.frag");
+    GLuint bubble_program =
+        loadShaders("shaders/bubble.vert", "shaders/bubble.frag");
     GLuint skybox_program =
         loadShaders("shaders/skybox.vert", "shaders/skybox.frag");
     printError("compiled shaders");
     ground = Object{ground_model, ground_program};
     surface = Surface{surface_model, surface_program};
     waterfall = Waterfall{waterfall_model, waterfall_program};
+    bubble = Bubble{bubble_model, bubble_program};
     skybox = Object{skybox_model, skybox_program};
 
     ground_fbo = initFBO2(1024, 1024, GL_LINEAR, true);
@@ -180,6 +186,26 @@ struct Scene {
     glEnable(GL_CULL_FACE);
   }
 
+  void draw_bubble() {
+    bubble.use();
+    GLuint program = bubble.program;
+    glUniformMatrix4fv(glGetUniformLocation(program, "projectionMatrix"), 1,
+                       GL_TRUE, proj_matrix.m);
+    glUniformMatrix4fv(glGetUniformLocation(program, "worldToView"), 1, GL_TRUE,
+                       view_matrix.m);
+    glUniformMatrix4fv(glGetUniformLocation(program, "modelToWorldToView"), 1,
+                       GL_TRUE, view_matrix.m);
+    glUniform3f(glGetUniformLocation(program, "camera_pos"), pos.x, pos.y,
+                pos.z);
+
+    // Skybox
+    glActiveTexture(GL_TEXTURE0);
+    glBindTexture(GL_TEXTURE_2D, skybox_tex);
+
+    bubble.draw();
+    bubble.move_bubble();
+  }
+
   void draw_waterfall() {
     waterfall.use();
     GLuint program = waterfall.program;
@@ -228,6 +254,7 @@ struct Scene {
     draw_surface();
     draw_ground();
     draw_waterfall();
+    draw_bubble();
   }
 };
 
diff --git a/src/waterfall.cpp b/src/waterfall.cpp
index 3f2e5c3fa114ad1de67e037c2fa47f133eaf57cf..f8be22d6b66f96087f1f19c174fb679f58d1cfe4 100644
--- a/src/waterfall.cpp
+++ b/src/waterfall.cpp
@@ -72,10 +72,8 @@ void Waterfall::move_waterfall_balls() {
   for (int i{0}; i < NUM_BALLS; ++i) {
     balls[i].pos += ball_velocities[i];
 
-    // TODO: check if balls are outside in x and z plane aswell
-    // maybe just generate a new ball if it falls outside
     if (balls[i].pos.y < 0) {
-      balls[i].pos.y = y.second;
+      balls[i].pos.y = y.second - balls[i].radius * 3.5;
     }
   }
 }
@@ -89,19 +87,21 @@ void Waterfall::gen_balls() {
     float radius = ((rand() / (float)RAND_MAX) * 0.9 + 0.1) / 10.0;
 
     // pseudo random value between min-radius and max-radius
-    float pos_x = lerp(x.first, x.second, (rand() / (float)RAND_MAX));
-    float pos_y = lerp(y.first, y.second, (rand() / (float)RAND_MAX));
-    float pos_z = lerp(z.first, z.second, (rand() / (float)RAND_MAX));
+    radius *= 3.5;
+    float pos_x =
+        lerp(x.first + radius, x.second - radius, (rand() / (float)RAND_MAX));
+    float pos_y =
+        lerp(y.first + radius, y.second - radius, (rand() / (float)RAND_MAX));
+    float pos_z =
+        lerp(z.first + radius, z.second - radius, (rand() / (float)RAND_MAX));
     vec4 pos = vec4(pos_x, pos_y, pos_z, 0);
+    radius /= 3.5;
 
     // pseudo random value between -0.001 and -0.01
-    float velocity_y = -((rand() / (float)RAND_MAX) * 0.9 + 0.1) / 100.0;
+    float velocity_y = -((rand() / (float)RAND_MAX) * 0.5 + 0.3) / 20.0;
     vec4 velocity = vec4(0, velocity_y, 0, 0);
 
     balls[i] = Waterfall::Ball(pos, radius);
     ball_velocities[i] = velocity;
   }
-  // balls[0].pos = vec4(3.2, 1.0, 0, 0);
-  // balls[0].radius = 0.55;
-  // ball_velocities[0] = vec4(0);
 }
diff --git a/src/waterfall.h b/src/waterfall.h
index 194c4303f6ff28da8301a9b416cd48690bd82759..9052606affff7e5ef7eb1a24be73cb738406bb66 100644
--- a/src/waterfall.h
+++ b/src/waterfall.h
@@ -13,7 +13,7 @@ struct Waterfall : Object {
 
   static GLuint const BINDING_POINT{0};
   char const *UNIFORM_BLOCK = "BallBuffer";
-  static int const NUM_BALLS{100};
+  static int const NUM_BALLS{250};
 
   // Dimensions of the waterfall
   // on the form <min, max>