Don't concatenate/append/stack arrays if you can help it, especially big ones. It's very wasteful of memory and slow.

Assign A = np.empty((m, n+2, n+2)) and then fill it with A[i] = np.r_[S1, np.c_[S2, Sc[i], S2], S1]. Or do it vectorized and get rid of the for loops:

A = np.zeros((m, n+2, n+2))
A[:,1:-1,1:-1] = Sc

or even do it in one line:

A = np.pad(Sc, ((0,0),(1,1),(1,1)), 'constant', constant_values = 0)

Sounds like what you really need is a python list of 3D numpy arrays. Appending to a numpy array is possible with np.append or np.concat, but it's very expensive because it forces the entire array to be remade. Is there any reason you want a 4D array?

Use dstack:

>>> np.dstack((A, B)).shape
(480, 640, 4)

This handles the cases where the arrays have different numbers of dimensions and stacks the arrays along the third axis.

Otherwise, to use append or concatenate, you'll have to make B three dimensional yourself and specify the axis you want to join them on:

>>> np.append(A, np.atleast_3d(B), axis=2).shape
(480, 640, 4)

The most efficient way is to initialize the 3d array first and then sequentially write your data into it. Efficiency here is increased by avoiding unnecessary copies. The gains are however minimal - as long as you are not working with very large arrays (in either of the 3 dimensions).

import numpy as np

arr = np.empty(shape=(2, 100, 100))

arr[0,:,:] = array_2d_1
arr[1,:,:] = array_2d_2

I recommend using numpy for multidimensional arrays. It makes it much more convenient, and much faster. This would look like:

import numpy as np
x = np.zeros((10,20,30)) # Make a 10 by 20 by 30 array
x[0,0,0] = value1

Still, if you don't want to use numpy, or need non-rectangular multi-dimensional arrays, you will need to treat it as a list of lists of lists, and initialize each list:

x = []
x.append([])
x[0].append([])
x[0][0].append(value1)

Edit: Or you could use the compact notation shown in ndpu's answer (x = [[[value1]]]).

You can allocate extra space for your array grid_data, fill it with the NaN, and keep track of the next index to be filled in another array while iterating through and filling with values from data. If you completely fill the third dimension for some lat_idx, lon_idx with non-NaN values, then you just allocate more space. Since appending is expensive with numpy, it's best that this extra space is pretty large so you only do it once or twice (below I allocate twice the original space).

Once the array is filled, you can remove the added space that was unused with numpy.isnan(). This solution does what you want but is very slow (for the example values you gave it took about two minutes), but the slow execution comes from iterating rather than the numpy operations.

Here's the code:

import random
import numpy as np

grid_data = np.ones(shape=(121, 201, 1000))
data = np.random.rand(4800, 4800)

# keep track of next index to fill for all the arrays in axis 2
next_to_fill = np.full(shape=(grid_data.shape[0], grid_data.shape[1]),
                       fill_value=grid_data.shape[2],
                       dtype=np.int32)

# allocate more space
double_shape = (grid_data.shape[0], grid_data.shape[1], grid_data.shape[2] * 2)
extra_space = np.full(shape=double_shape, fill_value=np.nan)
grid_data = np.append(grid_data, extra_space, axis=2)

for row in range(4800):
    for col in range(4800):
        lat_idx = random.randint(0, 120)
        lon_idx = random.randint(0, 200)

        # allocate more space if needed
        if next_to_fill[lat_idx, lon_idx] >= grid_data.shape[2]:
            grid_data = np.append(grid_data, extra_space, axis=2)

        grid_data[lat_idx, lon_idx, next_to_fill[lat_idx, lon_idx]] = data[row,
                                                                           col]
        next_to_fill[lat_idx, lon_idx] += 1

# remove unnecessary nans that were appended
not_all_nan_idxs = ~np.isnan(grid_data).all(axis=(0, 1))
grid_data = grid_data[:, :, not_all_nan_idxs]