Custom Matrix product

You can also achieve your ultimate (or ulterior) goal with ImageResize:

imgdata = Array[0.2 #1 + 0.1 #2 - 0.2 + 0.01 #3 &, {2, 2, 3}]
(* {{{0.11, 0.12, 0.13}, {0.21, 0.22, 0.23}}, {{0.31, 0.32, 0.33}, {0.41, 0.42, 0.43}}} *)

img = Image @ imgdata;
img2 = ImageResize[img, Scaled[2]];

ImageData[img2] // MatrixForm

P.S. There are various Resampling algorithms built into Mathematica.

One of the most direct ways hasn't been shown yet, which is to expand each element at level 2 with e.g. ConstantArray and then ArrayFlatten the entire result. Edit: Actually Nasser did this manually, without Map and using the less efficient Table, but the idea is identical.

ArrayFlatten @ Map[ConstantArray[#, {2, 4}] &, list1, {2}]

Slightly more complicated but significantly faster is to expand the entire array and then Flatten as necessary:

ConstantArray[#, {2, 4}] ~Flatten~ {{3, 1}, {4, 2}} & @ list1

On my system this tests faster than any other code I have tried, including ImageResize (see below).

---

Timings

Timings for all methods posted so far, in version 7.

Edit: Michael's ouput does not match the question or other answers. The code should be ImageResize[img, {Scaled[4], Scaled[2]}] which I will use below.

SetAttributes[timeAvg, HoldFirst]
timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]

list1 = RandomInteger[99, {500, 500, 3}];

ArrayFlatten @ Map[ConstantArray[#, {2, 4}] &, list1, {2}]    // timeAvg
ConstantArray[#, {2, 4}] ~Flatten~ {{3, 1}, {4, 2}} & @ list1 // timeAvg

ArrayFlatten[
  Transpose[Outer[Times, {{1, 1, 1, 1}, {1, 1, 1, 1}}, list1], {3, 4, 1, 2}]] // timeAvg

ImageData@ImageResize[Image@list1, {Scaled[4], Scaled[2]}] // timeAvg

Module[{f},
  f[x_] := ArrayFlatten@Table[x, {2}, {4}];
  ArrayFlatten@Map[f, list1, {2}]
] // timeAvg

With[{
   grid = ArrayFlatten[Table[ConstantArray[Slot@n, {2, 4}], {n, 4}]~Partition~2]
   },
  grid & @@ Flatten[list1, 1]
] // timeAvg

With[{exp = {2, 4}},
  Array[
   Extract[list1, Ceiling[{#1, #2}/exp]] &,
   exp Most@Dimensions[list1]
  ]
] // timeAvg

0.2714

0.03684

1.373

0.0512

1.529

0.0686

6.209

My second function takes first place, Michael's takes second, and rm-rf's takes third. Note that Michael's is less general, applying only to data that is handled by Image.

Observing the pattern in the desired output, you can construct something like this (extensible to larger grids):

With[{grid = ArrayFlatten[Table[ConstantArray[Slot@n, {2, 4}], {n, 4}] ~Partition~ 2]},
    grid & @@ Flatten[list, 1] 
] // MatrixForm

grid here is a pure function (well, only the slots) that looks like this, mimicking the structure of your output:

to which we pass the individual sublists as arguments.